FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_mpc_prepost.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5 
6 module hecmw_mpc_prepost
7  use hecmw_util
8  use m_hecmw_comm_f
9  use hecmw_matrix_misc
10  use hecmw_matrix_ass
11  use hecmw_local_matrix
12  use hecmw_solver_las
13  implicit none
14 
15  private
16  public :: hecmw_mpc_mat_init
17  public :: hecmw_mpc_mat_init_explicit
18  public :: hecmw_mpc_mat_finalize
19  public :: hecmw_mpc_mat_finalize_explicit
20  public :: hecmw_mpc_mat_ass
21  public :: hecmw_mpc_trans_rhs
22  public :: hecmw_mpc_tback_sol
23  public :: hecmw_mpc_trans_mass
24  public :: hecmw_mpc_tback_eigvec
25  public :: hecmw_mpc_mark_slave
26  public :: hecmw_mpc_check_solution
27  public :: hecmw_mpc_check_equation
28 
29  integer, parameter :: DEBUG = 0
30  logical, parameter :: DEBUG_VECTOR = .false.
31 
32 contains
33 
34  !C
35  !C***
36  !C*** hecmw_mpc_mat_init
37  !C***
38  !C
39  subroutine hecmw_mpc_mat_init(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMAT, conMATmpc)
40  implicit none
41  type (hecmwst_local_mesh), intent(inout), target :: hecmesh
42  type (hecmwst_matrix), intent(in), target :: hecmat
43  type (hecmwst_local_mesh), pointer :: hecmeshmpc
44  type (hecmwst_matrix), pointer :: hecmatmpc
45  type (hecmwst_matrix), intent(in), target, optional :: conmat
46  type (hecmwst_matrix), pointer, optional :: conmatmpc
47  integer(kind=kint) :: totalmpc, mpc_method, solver_type
48 
49  totalmpc = hecmesh%mpc%n_mpc
50  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
51 
52  if (totalmpc == 0) then
53  hecmeshmpc => hecmesh
54  hecmatmpc => hecmat
55  if (present(conmat).and.present(conmatmpc)) conmatmpc => conmat
56  return
57  endif
58 
59  call hecmw_mpc_scale(hecmesh)
60 
61  mpc_method = hecmw_mat_get_mpc_method(hecmat)
62  if (mpc_method < 1 .or. 3 < mpc_method) then
63  solver_type = hecmw_mat_get_solver_type(hecmat)
64  if (solver_type > 1) then ! DIRECT SOLVER
65  mpc_method = 1 ! default: penalty
66  else ! ITERATIVE SOLVER
67  mpc_method = 3 ! default: elimination
68  endif
69  call hecmw_mat_set_mpc_method(hecmat, mpc_method)
70  endif
71 
72  if (mpc_method == 2) then
73  write(*,*) 'WARNING: MPCMETHOD=2 (MPCCG) is deprecated; may not work correctly'
74  ! MPC_METHOD = 3
75  ! call hecmw_mat_set_mpc_method(hecMAT, MPC_METHOD)
76  endif
77 
78  select case (mpc_method)
79  case (1) ! penalty
80  hecmeshmpc => hecmesh
81  hecmatmpc => hecmat
82  if (present(conmat).and.present(conmatmpc)) conmatmpc => conmat
83  case (2) ! MPCCG
84  hecmeshmpc => hecmesh
85  hecmatmpc => hecmat
86  if (present(conmat).and.present(conmatmpc)) conmatmpc => conmat
87  case (3) ! elimination
88  allocate(hecmeshmpc)
89  call hecmw_mpc_mesh_copy(hecmesh, hecmeshmpc)
90  allocate(hecmatmpc)
91  call hecmw_mat_init(hecmatmpc)
92  if (present(conmat).and.present(conmatmpc)) then
93  allocate(conmatmpc)
94  call hecmw_mat_init(conmatmpc)
95  endif
96  end select
97 
98  end subroutine hecmw_mpc_mat_init
99 
100  !C
101  !C***
102  !C*** hecmw_mpc_mat_init_explicit
103  !C***
104  !C
105  subroutine hecmw_mpc_mat_init_explicit(hecMESH, hecMAT, hecMATmpc)
106  implicit none
107  type (hecmwst_local_mesh), intent(inout), target :: hecmesh
108  type (hecmwst_matrix), intent(in), target :: hecmat
109  type (hecmwst_matrix), pointer :: hecmatmpc
110  integer(kind=kint) :: totalmpc, mpc_method
111 
112  totalmpc = hecmesh%mpc%n_mpc
113  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
114 
115  if (totalmpc == 0) then
116  hecmatmpc => hecmat
117  return
118  endif
119 
120  call hecmw_mpc_scale(hecmesh)
121 
122  ! Force MPC_METHOD=3
123  mpc_method = 3
124  call hecmw_mat_set_mpc_method(hecmat, mpc_method)
125 
126  allocate(hecmatmpc)
127  call hecmw_mat_init(hecmatmpc)
128 
129  hecmatmpc%N = hecmat%N
130  hecmatmpc%NP = hecmat%NP
131  hecmatmpc%NDOF = hecmat%NDOF
132  allocate(hecmatmpc%B(size(hecmat%B)))
133  allocate(hecmatmpc%X(size(hecmat%X)))
134  end subroutine hecmw_mpc_mat_init_explicit
135 
136  !C
137  !C***
138  !C*** hecmw_mpc_mat_finalize
139  !C***
140  !C
141  subroutine hecmw_mpc_mat_finalize(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMATmpc)
142  implicit none
143  type (hecmwst_local_mesh), intent(in) :: hecmesh
144  type (hecmwst_matrix), intent(in) :: hecmat
145  type (hecmwst_local_mesh), pointer :: hecmeshmpc
146  type (hecmwst_matrix), pointer :: hecmatmpc
147  type (hecmwst_matrix), pointer, optional :: conmatmpc
148  integer(kind=kint) :: totalmpc, mpc_method
149 
150  totalmpc = hecmesh%mpc%n_mpc
151  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
152 
153  if (totalmpc == 0) then
154  nullify(hecmeshmpc)
155  nullify(hecmatmpc)
156  if (present(conmatmpc)) nullify(conmatmpc)
157  return
158  endif
159 
160  mpc_method = hecmw_mat_get_mpc_method(hecmat)
161 
162  select case (mpc_method)
163  case (1) ! penalty
164  nullify(hecmeshmpc)
165  nullify(hecmatmpc)
166  if (present(conmatmpc)) nullify(conmatmpc)
167  case (2) ! MPCCG
168  nullify(hecmeshmpc)
169  nullify(hecmatmpc)
170  if (present(conmatmpc)) nullify(conmatmpc)
171  case (3) ! elimination
172  call hecmw_mpc_mesh_free(hecmeshmpc)
173  deallocate(hecmeshmpc)
174  nullify(hecmeshmpc)
175  call hecmw_mat_finalize(hecmatmpc)
176  deallocate(hecmatmpc)
177  nullify(hecmatmpc)
178  if (present(conmatmpc)) then
179  call hecmw_mat_finalize(conmatmpc)
180  deallocate(conmatmpc)
181  nullify(conmatmpc)
182  endif
183  end select
184 
185  end subroutine hecmw_mpc_mat_finalize
186 
187  !C
188  !C***
189  !C*** hecmw_mpc_mat_finalize_explicit
190  !C***
191  !C
192  subroutine hecmw_mpc_mat_finalize_explicit(hecMESH, hecMAT, hecMATmpc)
193  implicit none
194  type (hecmwst_local_mesh), intent(in) :: hecmesh
195  type (hecmwst_matrix), intent(in) :: hecmat
196  type (hecmwst_matrix), pointer :: hecmatmpc
197  integer(kind=kint) :: totalmpc, mpc_method
198 
199  totalmpc = hecmesh%mpc%n_mpc
200  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
201 
202  if (totalmpc == 0) then
203  nullify(hecmatmpc)
204  return
205  endif
206 
207  mpc_method = hecmw_mat_get_mpc_method(hecmat)
208 
209  select case (mpc_method)
210  case (1) ! penalty
211  nullify(hecmatmpc)
212  case (2) ! MPCCG
213  nullify(hecmatmpc)
214  case (3) ! elimination
215  call hecmw_mat_finalize(hecmatmpc)
216  deallocate(hecmatmpc)
217  nullify(hecmatmpc)
218  end select
219 
220  end subroutine hecmw_mpc_mat_finalize_explicit
221 
222  !C
223  !C***
224  !C*** hecmw_mpc_mat_ass
225  !C***
226  !C
227  subroutine hecmw_mpc_mat_ass(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMAT, conMATmpc, hecLagMAT)
228  implicit none
229  type (hecmwst_local_mesh), intent(in) :: hecmesh
230  type (hecmwst_matrix), intent(inout) :: hecmat
231  type (hecmwst_local_mesh), pointer :: hecmeshmpc
232  type (hecmwst_matrix), pointer :: hecmatmpc
233  type (hecmwst_matrix), intent(in), optional :: conmat
234  type (hecmwst_matrix), pointer, optional :: conmatmpc
235  type (hecmwst_matrix_lagrange), intent(inout), optional :: heclagmat
236  integer(kind=kint) :: totalmpc, mpc_method
237 
238  totalmpc = hecmesh%mpc%n_mpc
239  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
240 
241  if (totalmpc == 0) return
242 
243  mpc_method = hecmw_mat_get_mpc_method(hecmat)
244 
245  select case (mpc_method)
246  case (1) ! penalty
247  !if (hecMESH%my_rank.eq.0) write(0,*) "MPC Method: Penalty"
248  call hecmw_mat_ass_equation ( hecmesh, hecmat )
249  case (2) ! MPCCG
250  !if (hecMESH%my_rank.eq.0) write(0,*) "MPC Method: MPC-CG"
251  case (3) ! elimination
252  !if (hecMESH%my_rank.eq.0) write(0,*) "MPC Method: Elimination"
253  call hecmw_trimatmul_ttkt_mpc(hecmeshmpc, hecmat, hecmatmpc)
254  if (present(conmat).and.present(conmatmpc).and.present(heclagmat)) then
255  call hecmw_trimatmul_ttkt_mpc(hecmeshmpc, conmat, conmatmpc)
256  call resize_heclagmat(conmat%NP, conmatmpc%NP, conmat%NDOF, heclagmat)
257  endif
258  end select
259 
260  end subroutine hecmw_mpc_mat_ass
261 
262 
263  subroutine resize_heclagmat(NP_orig, NP_new, ndof, hecLagMAT)
264  integer(kind=kint), intent(in) :: np_orig, np_new, ndof
265  type (hecmwst_matrix_lagrange), intent(inout) :: heclagmat
266  integer(kind=kint), pointer :: itemp(:)
267 
268  if (heclagmat%num_lagrange == 0) return
269 
270  allocate(itemp(0:np_new))
271  itemp(0:np_orig) = heclagmat%indexU_lagrange(0:np_orig)
272  itemp(np_orig+1:np_new) = heclagmat%indexU_lagrange(np_orig)
273 
274  deallocate(heclagmat%indexU_lagrange)
275  heclagmat%indexU_lagrange => itemp
276 
277  end subroutine resize_heclagmat
278 
279  !C
280  !C***
281  !C*** hecmw_mpc_trans_rhs
282  !C***
283  !C
284  subroutine hecmw_mpc_trans_rhs(hecMESH, hecMAT, hecMATmpc)
285  implicit none
286  type (hecmwst_local_mesh), intent(inout) :: hecmesh
287  type (hecmwst_matrix), intent(inout) :: hecmat
288  type (hecmwst_matrix), pointer :: hecmatmpc
289  real(kind=kreal), allocatable :: btmp(:)
290  real(kind=kreal) :: time_dumm
291  integer(kind=kint) :: totalmpc, mpc_method, i
292 
293  totalmpc = hecmesh%mpc%n_mpc
294  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
295 
296  if (totalmpc == 0) return
297 
298  mpc_method = hecmw_mat_get_mpc_method(hecmat)
299 
300  select case (mpc_method)
301  case (1) ! penalty
302  call hecmw_mat_ass_equation_rhs ( hecmesh, hecmatmpc )
303  case (2) ! MPCCG
304  allocate(btmp(hecmat%NP*hecmat%NDOF))
305  do i = 1, hecmat%NP*hecmat%NDOF
306  btmp(i) = hecmat%B(i)
307  enddo
308  call hecmw_trans_b(hecmesh, hecmat, btmp, hecmatmpc%B, time_dumm)
309  deallocate(btmp)
310  case (3) ! elimination
311  call hecmw_trans_b(hecmesh, hecmat, hecmat%B, hecmatmpc%B, time_dumm)
312  hecmatmpc%Iarray=hecmat%Iarray
313  hecmatmpc%Rarray=hecmat%Rarray
314  end select
315 
316  end subroutine hecmw_mpc_trans_rhs
317 
318  !C
319  !C***
320  !C*** hecmw_mpc_tback_sol
321  !C***
322  !C
323  subroutine hecmw_mpc_tback_sol(hecMESH, hecMAT, hecMATmpc)
324  implicit none
325  type (hecmwst_local_mesh), intent(in) :: hecmesh
326  type (hecmwst_matrix), intent(inout) :: hecmat
327  type (hecmwst_matrix), pointer :: hecmatmpc
328  real(kind=kreal) :: time_dumm
329  integer(kind=kint) :: totalmpc, mpc_method, i
330  integer(kind=kint) :: npndof, npndof_mpc, num_lagrange
331 
332  totalmpc = hecmesh%mpc%n_mpc
333  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
334 
335  if (totalmpc == 0) return
336 
337  mpc_method = hecmw_mat_get_mpc_method(hecmat)
338 
339  select case (mpc_method)
340  case (1) ! penalty
341  ! do nothing
342  case (2) ! MPCCG
343  call hecmw_tback_x(hecmesh, hecmat%NDOF, hecmat%X, time_dumm)
344  case (3) ! elimination
345  npndof = hecmat%NP * hecmat%NDOF
346  do i = 1, npndof
347  hecmat%X(i) = hecmatmpc%X(i)
348  enddo
349  call hecmw_tback_x(hecmesh, hecmat%NDOF, hecmat%X, time_dumm)
350  num_lagrange = size(hecmat%X) - npndof
351  npndof_mpc = hecmatmpc%NP * hecmatmpc%NDOF
352  do i = 1, num_lagrange
353  hecmat%X(npndof+i) = hecmatmpc%X(npndof_mpc+i)
354  enddo
355  hecmat%Iarray=hecmatmpc%Iarray
356  hecmat%Rarray=hecmatmpc%Rarray
357  end select
358  end subroutine hecmw_mpc_tback_sol
359 
360  !C
361  !C***
362  !C*** hecmw_mpc_trans_mass
363  !C***
364  !C
365  subroutine hecmw_mpc_trans_mass(hecMESH, hecMAT, mass)
366  implicit none
367  type (hecmwst_local_mesh), intent(inout) :: hecmesh
368  type (hecmwst_matrix), intent(inout) :: hecmat
369  real(kind=kreal), intent(inout) :: mass(:)
370 
371  real(kind=kreal), allocatable :: mtmp(:)
372  real(kind=kreal) :: time_dumm
373  integer(kind=kint) :: totalmpc, mpc_method, i
374 
375  totalmpc = hecmesh%mpc%n_mpc
376  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
377 
378  if (totalmpc == 0) return
379 
380  mpc_method = hecmw_mat_get_mpc_method(hecmat)
381 
382  select case (mpc_method)
383  case (1) ! penalty
384  ! do nothing
385  case (2,3) ! MPCCG or elimination
386  allocate(mtmp(hecmat%NP*hecmat%NDOF))
387  !C-- {Mt} = [T'] {w}
388  call hecmw_ttvec(hecmesh, hecmat%NDOF, mass, mtmp, time_dumm)
389  do i = 1, hecmat%NP*hecmat%NDOF
390  mass(i) = mtmp(i)
391  enddo
392  deallocate(mtmp)
393  end select
394 
395  end subroutine hecmw_mpc_trans_mass
396 
397  !C
398  !C***
399  !C*** hecmw_mpc_tback_eigvec
400  !C***
401  !C
402  subroutine hecmw_mpc_tback_eigvec(hecMESH, hecMAT, neig, eigvec)
403  implicit none
404  type (hecmwst_local_mesh), intent(in) :: hecmesh
405  type (hecmwst_matrix), intent(inout) :: hecmat
406  integer(kind=kint), intent(in) :: neig
407  real(kind=kreal), intent(inout) :: eigvec(:,:)
408 
409  real(kind=kreal) :: time_dumm
410  integer(kind=kint) :: totalmpc, mpc_method, i
411 
412  totalmpc = hecmesh%mpc%n_mpc
413  call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
414 
415  if (totalmpc == 0) return
416 
417  mpc_method = hecmw_mat_get_mpc_method(hecmat)
418 
419  select case (mpc_method)
420  case (1) ! penalty
421  ! do nothing
422  case (2,3) ! MPCCG or elimination
423  do i = 1, neig
424  call hecmw_tback_x(hecmesh, hecmat%NDOF, eigvec(:,i), time_dumm)
425  !!! need normalization???
426  enddo
427  end select
428  end subroutine hecmw_mpc_tback_eigvec
429 
430  !C
431  !C***
432  !C*** hecmw_mpc_mark_slave
433  !C***
434  !C
435  subroutine hecmw_mpc_mark_slave(hecMESH, hecMAT, mark)
436  implicit none
437  type (hecmwst_local_mesh), intent(in) :: hecmesh
438  type (hecmwst_matrix), intent(inout) :: hecmat
439  integer(kind=kint), intent(out) :: mark(:)
440 
441  integer(kind=kint) :: ndof, i, j, k, kk
442 
443  ndof = hecmat%NDOF
444  mark(:) = 0
445  outer: do i = 1, hecmesh%mpc%n_mpc
446  do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
447  if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
448  enddo
449  k = hecmesh%mpc%mpc_index(i-1)+1
450  kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
451  mark(kk) = 1
452  enddo outer
453  end subroutine hecmw_mpc_mark_slave
454 
455  !C
456  !C***
457  !C*** hecmw_mpc_scale
458  !C***
459  !C
460  subroutine hecmw_mpc_scale(hecMESH)
461  implicit none
462  type (hecmwst_local_mesh), intent(inout) :: hecmesh
463  integer(kind=kint) :: i, j, k
464  real(kind=kreal) :: wval
465 
466  !$omp parallel default(none),private(i,j,k,WVAL),shared(hecMESH)
467  !$omp do
468  do i = 1, hecmesh%mpc%n_mpc
469  k = hecmesh%mpc%mpc_index(i-1)+1
470  wval = 1.d0 / hecmesh%mpc%mpc_val(k)
471  hecmesh%mpc%mpc_val(k) = 1.d0
472  do j = hecmesh%mpc%mpc_index(i-1)+2, hecmesh%mpc%mpc_index(i)
473  hecmesh%mpc%mpc_val(j) = hecmesh%mpc%mpc_val(j) * wval
474  enddo
475  hecmesh%mpc%mpc_const(i) = hecmesh%mpc%mpc_const(i) * wval
476  enddo
477  !$omp end do
478  !$omp end parallel
479 
480  end subroutine hecmw_mpc_scale
481 
482 
483  !C
484  !C***
485  !C*** hecmw_trans_b
486  !C***
487  !C
488  subroutine hecmw_trans_b(hecMESH, hecMAT, B, BT, COMMtime)
489  implicit none
490  type (hecmwst_local_mesh), intent(in) :: hecmesh
491  type (hecmwst_matrix), intent(in) :: hecmat
492  real(kind=kreal), intent(in) :: b(:)
493  real(kind=kreal), intent(out), target :: bt(:)
494  real(kind=kreal), intent(inout) :: commtime
495 
496  real(kind=kreal), allocatable :: w(:)
497  real(kind=kreal), pointer :: xg(:)
498  integer(kind=kint) :: ndof, i, j, k, kk, flg_bak
499 
500  ndof = hecmat%NDOF
501 
502  call debug_write_vector(b, 'original RHS', 'B', ndof, hecmat%N, hecmat%NP, .true.)
503 
504  allocate(w(hecmesh%n_node * ndof))
505 
506  !C===
507  !C +---------------------------+
508  !C | {bt}= [T']({b} - [A]{xg}) |
509  !C +---------------------------+
510  !C===
511  xg => bt
512  do i = 1, hecmat%N * ndof
513  xg(i) = 0.d0
514  enddo
515 
516  !C-- Generate {xg} from mpc_const
517  !$omp parallel default(none),private(i,k,kk),shared(hecMESH,XG),firstprivate(ndof)
518  !$omp do
519  outer: do i = 1, hecmesh%mpc%n_mpc
520  do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
521  if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
522  enddo
523  k = hecmesh%mpc%mpc_index(i-1) + 1
524  kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
525  xg(kk) = hecmesh%mpc%mpc_const(i)
526  enddo outer
527  !$omp end do
528  !$omp end parallel
529 
530  !C-- {w} = {b} - [A]{xg}
531  flg_bak = hecmw_mat_get_flag_mpcmatvec(hecmat)
532  call hecmw_mat_set_flag_mpcmatvec(hecmat, 0)
533  call hecmw_matresid(hecmesh, hecmat, xg, b, w, commtime)
534  call hecmw_mat_set_flag_mpcmatvec(hecmat, flg_bak)
535 
536  !C-- {bt} = [T'] {w}
537  call hecmw_ttvec(hecmesh, ndof, w, bt, commtime)
538 
539  deallocate(w)
540 
541  call debug_write_vector(bt, 'transformed RHS', 'BT', ndof, hecmat%N, hecmat%NP, .true.)
542  end subroutine hecmw_trans_b
543 
544 
545  !C
546  !C***
547  !C*** hecmw_tback_x
548  !C***
549  !C
550  subroutine hecmw_tback_x(hecMESH, ndof, X, COMMtime)
551  implicit none
552  type (hecmwst_local_mesh), intent(in) :: hecmesh
553  integer(kind=kint), intent(in) :: ndof
554  real(kind=kreal), intent(inout) :: x(:)
555  real(kind=kreal), intent(inout) :: commtime
556 
557  real(kind=kreal), allocatable :: w(:)
558  integer(kind=kint) :: i, j, k, kk
559 
560  call debug_write_vector(x, 'solution for transformed eqn', 'X', ndof, hecmesh%nn_internal, &
561  hecmesh%n_node, .true.)
562 
563  allocate(w(hecmesh%n_node * ndof))
564 
565  !C-- {tx} = [T]{x}
566  call hecmw_tvec(hecmesh, ndof, x, w, commtime)
567 
568  !C-- {x} = {tx} + {xg}
569  !$omp parallel default(none),private(i,k,kk),shared(hecMESH,X,W),firstprivate(ndof)
570  !$omp do
571  do i= 1, hecmesh%nn_internal * ndof
572  x(i)= w(i)
573  enddo
574  !$omp end do
575 
576  !$omp do
577  outer: do i = 1, hecmesh%mpc%n_mpc
578  do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
579  if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
580  enddo
581  k = hecmesh%mpc%mpc_index(i-1) + 1
582  kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
583  x(kk) = x(kk) + hecmesh%mpc%mpc_const(i)
584  enddo outer
585  !$omp end do
586  !$omp end parallel
587 
588  deallocate(w)
589 
590  call hecmw_update_r(hecmesh, x, hecmesh%n_node, ndof)
591  call debug_write_vector(x, 'recovered solution', 'X', ndof, hecmesh%nn_internal, &
592  hecmesh%n_node, .true.)
593  end subroutine hecmw_tback_x
594 
595  subroutine hecmw_mpc_check_solution(hecMESH, hecMAT)
596  use hecmw_solver_misc
597  implicit none
598  type (hecmwst_local_mesh), intent(in) :: hecmesh
599  type (hecmwst_matrix), intent(in) :: hecmat
600 
601  real(kind=kreal), allocatable :: kx(:), ttkx(:), ttb(:), resid(:)
602  integer(kind=kint) :: ndof, nndof, npndof, i
603  real(kind=kreal) :: bnrm, rnrm, commtime
604 
605  ndof = hecmat%NDOF
606  nndof = hecmat%N * ndof
607  npndof = hecmat%NP * ndof
608  allocate(kx(npndof), ttkx(npndof), ttb(npndof), resid(npndof))
609  call hecmw_matvec(hecmesh, hecmat, hecmat%X, kx)
610  call hecmw_ttvec(hecmesh, ndof, kx, ttkx, commtime)
611  call hecmw_ttvec(hecmesh, ndof, hecmat%B, ttb, commtime)
612  do i = 1, nndof
613  resid(i) = ttb(i) - ttkx(i)
614  enddo
615  call hecmw_innerproduct_r(hecmesh, ndof, ttb, ttb, bnrm)
616  call hecmw_innerproduct_r(hecmesh, ndof, resid, resid, rnrm)
617  bnrm = sqrt(bnrm)
618  rnrm = sqrt(rnrm)
619  write(0,*) 'DEBUG: hecmw_mpc_check_solution: resid(abs), resid(rel):',rnrm, rnrm/bnrm
620  end subroutine hecmw_mpc_check_solution
621 
622  subroutine hecmw_mpc_check_equation(hecMESH, hecMAT)
623  implicit none
624  type (hecmwst_local_mesh), intent(in) :: hecmesh
625  type (hecmwst_matrix), intent(in) :: hecmat
626 
627  integer(kind=kint) :: ndof
628  integer(kind=kint) :: i, j, kk
629  real(kind=kreal) :: xmax, sum, resid, resid_max
630 
631  ndof = hecmat%NDOF
632  xmax = 0.d0
633  do i = 1, hecmesh%nn_internal * ndof
634  if (xmax < abs(hecmat%X(i))) xmax = abs(hecmat%X(i))
635  enddo
636  call hecmw_allreduce_r1 (hecmesh, xmax, hecmw_max)
637 
638  resid_max = 0.d0
639  outer: do i = 1, hecmesh%mpc%n_mpc
640  do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
641  if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
642  enddo
643  sum = 0.d0
644  do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
645  kk = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
646  sum = sum + hecmat%X(kk) * hecmesh%mpc%mpc_val(j)
647  enddo
648  resid = abs(hecmesh%mpc%mpc_const(i) - sum)
649  if (resid_max < resid) resid_max = resid
650  enddo outer
651  call hecmw_allreduce_r1 (hecmesh, resid_max, hecmw_max)
652  write(0,*) 'DEBUG: hecmw_mpc_check_equation: resid(abs), resid(rel):', resid_max, resid_max/xmax
653  end subroutine hecmw_mpc_check_equation
654 
655  subroutine hecmw_mpc_mesh_copy(src, dst)
656  implicit none
657  type (hecmwst_local_mesh), intent(in) :: src
658  type (hecmwst_local_mesh), intent(out) :: dst
659  dst%zero = src%zero
660  dst%MPI_COMM = src%MPI_COMM
661  dst%PETOT = src%PETOT
662  dst%PEsmpTOT = src%PEsmpTOT
663  dst%my_rank = src%my_rank
664  dst%n_subdomain = src%n_subdomain
665  dst%n_node = src%n_node
666  dst%nn_internal = src%nn_internal
667  dst%n_elem = src%n_elem
668  dst%ne_internal = src%ne_internal
669  dst%n_elem_type = src%n_elem_type
670  dst%n_dof = src%n_dof
671  dst%n_neighbor_pe = src%n_neighbor_pe
672  if (src%n_neighbor_pe > 0) then
673  allocate(dst%neighbor_pe(dst%n_neighbor_pe))
674  dst%neighbor_pe(:) = src%neighbor_pe(:)
675  allocate(dst%import_index(0:dst%n_neighbor_pe))
676  dst%import_index(:)= src%import_index(:)
677  allocate(dst%export_index(0:dst%n_neighbor_pe))
678  dst%export_index(:)= src%export_index(:)
679  allocate(dst%import_item(dst%import_index(dst%n_neighbor_pe)))
680  dst%import_item(:) = src%import_item(:)
681  allocate(dst%export_item(dst%export_index(dst%n_neighbor_pe)))
682  dst%export_item(:) = src%export_item(:)
683  endif
684  allocate(dst%global_node_ID(dst%n_node))
685  dst%global_node_ID(1:dst%n_node) = src%global_node_ID(1:dst%n_node)
686  allocate(dst%node_ID(2*dst%n_node))
687  dst%node_ID(1:2*dst%n_node) = src%node_ID(1:2*dst%n_node)
688  allocate(dst%elem_type_item(dst%n_elem_type))
689  dst%elem_type_item(:) = src%elem_type_item(:)
690  !
691  dst%mpc%n_mpc = src%mpc%n_mpc
692  dst%mpc%mpc_index => src%mpc%mpc_index
693  dst%mpc%mpc_item => src%mpc%mpc_item
694  dst%mpc%mpc_dof => src%mpc%mpc_dof
695  dst%mpc%mpc_val => src%mpc%mpc_val
696  dst%mpc%mpc_const => src%mpc%mpc_const
697  !
698  dst%node_group%n_grp = src%node_group%n_grp
699  dst%node_group%n_bc = src%node_group%n_bc
700  dst%node_group%grp_name => src%node_group%grp_name
701  dst%node_group%grp_index => src%node_group%grp_index
702  dst%node_group%grp_item => src%node_group%grp_item
703  dst%node_group%bc_grp_ID => src%node_group%bc_grp_ID
704  dst%node_group%bc_grp_type => src%node_group%bc_grp_type
705  dst%node_group%bc_grp_index => src%node_group%bc_grp_index
706  dst%node_group%bc_grp_dof => src%node_group%bc_grp_dof
707  dst%node_group%bc_grp_val => src%node_group%bc_grp_val
708  !
709  dst%node => src%node
710  end subroutine hecmw_mpc_mesh_copy
711 
712  subroutine hecmw_mpc_mesh_free(hecMESH)
713  implicit none
714  type (hecmwst_local_mesh), intent(inout) :: hecmesh
715  if (hecmesh%n_neighbor_pe > 1) then
716  deallocate(hecmesh%neighbor_pe)
717  deallocate(hecmesh%import_index)
718  deallocate(hecmesh%export_index)
719  deallocate(hecmesh%import_item)
720  deallocate(hecmesh%export_item)
721  endif
722  deallocate(hecmesh%global_node_ID)
723  deallocate(hecmesh%node_ID)
724  deallocate(hecmesh%elem_type_item)
725  end subroutine hecmw_mpc_mesh_free
726 
729  subroutine debug_write_vector(Vec, label, name, ndof, N, &
730  NP, write_ext, slaves)
731  real(kind=kreal), intent(in) :: vec(:)
732  character(len=*), intent(in) :: label
733  character(len=*), intent(in) :: name
734  integer(kind=kint), intent(in) :: ndof
735  integer(kind=kint), intent(in) :: n
736  integer(kind=kint), intent(in), optional :: np
737  logical, intent(in), optional :: write_ext
738  integer(kind=kint), intent(in), optional :: slaves(:)
739  !
740  integer(kind=kint) :: iunit
741  character(len=128) :: fmt
742 
743  if (.not. debug_vector) return
744 
745  write(fmt,'(a,i0,a)') '(',ndof,'f12.3)'
746 
747  iunit = 1000 + hecmw_comm_get_rank()
748  write(iunit,*) trim(label),'------------------------------------------------------------'
749  write(iunit,*) 'size of ',trim(name),size(vec)
750  write(iunit,*) trim(name),': 1-',n*ndof
751  write(iunit,fmt) vec(1:n*ndof)
752  if (present(write_ext) .and. present(np)) then
753  if (write_ext) then
754  write(iunit,*) trim(name),'(external): ',n*ndof+1,'-',np*ndof
755  write(iunit,fmt) vec(n*ndof+1:np*ndof)
756  endif
757  endif
758  if (present(slaves)) then
759  if (size(slaves) > 0) then
760  write(iunit,*) trim(name),'(slave):',slaves(:)
761  write(iunit,fmt) vec(slaves(:))
762  endif
763  endif
764  end subroutine debug_write_vector
765 end module hecmw_mpc_prepost
hecmw_local_matrix
Definition: hecmw_local_matrix.f90:6
hecmw_mpc_prepost
Definition: hecmw_mpc_prepost.f90:6
hecmw_util::hecmw_max
integer(kind=kint), parameter hecmw_max
Definition: hecmw_util_f.F90:25
hecmw_matrix_misc::hecmw_mat_init
subroutine, public hecmw_mat_init(hecMAT)
Definition: hecmw_matrix_misc.f90:162
hecmw_mpc_prepost::hecmw_mpc_mat_finalize
subroutine, public hecmw_mpc_mat_finalize(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMATmpc)
Definition: hecmw_mpc_prepost.f90:142
hecmw_matrix_misc::hecmw_mat_get_solver_type
integer(kind=kint) function, public hecmw_mat_get_solver_type(hecMAT)
Definition: hecmw_matrix_misc.f90:601
hecmw_solver_las::hecmw_matvec
subroutine, public hecmw_matvec(hecMESH, hecMAT, X, Y, COMMtime)
Definition: hecmw_solver_las.f90:43
hecmw_matrix_misc::hecmw_mat_get_flag_mpcmatvec
integer(kind=kint) function, public hecmw_mat_get_flag_mpcmatvec(hecMAT)
Definition: hecmw_matrix_misc.f90:643
hecmw_mpc_prepost::hecmw_mpc_mark_slave
subroutine, public hecmw_mpc_mark_slave(hecMESH, hecMAT, mark)
Definition: hecmw_mpc_prepost.f90:436
hecmw_util::hecmw_sum
integer(kind=kint), parameter hecmw_sum
Definition: hecmw_util_f.F90:23
hecmw_mpc_prepost::hecmw_mpc_mat_ass
subroutine, public hecmw_mpc_mat_ass(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMAT, conMATmpc, hecLagMAT)
Definition: hecmw_mpc_prepost.f90:228
hecmw_solver_las
Definition: hecmw_solver_las.f90:6
hecmw_matrix_misc::hecmw_mat_set_mpc_method
subroutine, public hecmw_mat_set_mpc_method(hecMAT, mpc_method)
Definition: hecmw_matrix_misc.f90:429
hecmw_mpc_prepost::hecmw_mpc_mat_init
subroutine, public hecmw_mpc_mat_init(hecMESH, hecMAT, hecMESHmpc, hecMATmpc, conMAT, conMATmpc)
Definition: hecmw_mpc_prepost.f90:40
hecmw_solver_las::hecmw_matresid
subroutine, public hecmw_matresid(hecMESH, hecMAT, X, B, R, COMMtime)
Definition: hecmw_solver_las.f90:95
hecmw_solver_misc::hecmw_innerproduct_r
subroutine hecmw_innerproduct_r(hecMESH, ndof, X, Y, sum, COMMtime)
Definition: hecmw_solver_misc.f90:47
hecmw_mpc_prepost::hecmw_mpc_trans_rhs
subroutine, public hecmw_mpc_trans_rhs(hecMESH, hecMAT, hecMATmpc)
Definition: hecmw_mpc_prepost.f90:285
hecmw_mpc_prepost::hecmw_mpc_check_equation
subroutine, public hecmw_mpc_check_equation(hecMESH, hecMAT)
Definition: hecmw_mpc_prepost.f90:623
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::hecmwst_local_mesh
Definition: hecmw_util_f.F90:234
m_hecmw_comm_f::hecmw_allreduce_r1
subroutine hecmw_allreduce_r1(hecMESH, s, ntag)
Definition: hecmw_comm_f.F90:403
hecmw_mpc_prepost::hecmw_mpc_tback_eigvec
subroutine, public hecmw_mpc_tback_eigvec(hecMESH, hecMAT, neig, eigvec)
Definition: hecmw_mpc_prepost.f90:403
m_hecmw_comm_f
Definition: hecmw_comm_f.F90:6
hecmw_mpc_prepost::hecmw_mpc_mat_init_explicit
subroutine, public hecmw_mpc_mat_init_explicit(hecMESH, hecMAT, hecMATmpc)
Definition: hecmw_mpc_prepost.f90:106
hecmw_solver_las::hecmw_ttvec
subroutine, public hecmw_ttvec(hecMESH, ndof, X, Y, COMMtime)
Definition: hecmw_solver_las.f90:180
hecmw_mpc_prepost::hecmw_mpc_mat_finalize_explicit
subroutine, public hecmw_mpc_mat_finalize_explicit(hecMESH, hecMAT, hecMATmpc)
Definition: hecmw_mpc_prepost.f90:193
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_hecmw_comm_f::hecmw_allreduce_i1
subroutine hecmw_allreduce_i1(hecMESH, s, ntag)
Definition: hecmw_comm_f.F90:458
hecmw_matrix_ass
Definition: hecmw_mat_ass.f90:6
hecmw_mpc_prepost::hecmw_mpc_trans_mass
subroutine, public hecmw_mpc_trans_mass(hecMESH, hecMAT, mass)
Definition: hecmw_mpc_prepost.f90:366
hecmw_mpc_prepost::hecmw_mpc_tback_sol
subroutine, public hecmw_mpc_tback_sol(hecMESH, hecMAT, hecMATmpc)
Definition: hecmw_mpc_prepost.f90:324
hecmw_solver_las::hecmw_tvec
subroutine, public hecmw_tvec(hecMESH, ndof, X, Y, COMMtime)
Definition: hecmw_solver_las.f90:154
hecmw_matrix_misc::hecmw_mat_set_flag_mpcmatvec
subroutine, public hecmw_mat_set_flag_mpcmatvec(hecMAT, flag_mpcmatvec)
Definition: hecmw_matrix_misc.f90:637
hecmw_matrix_misc
Definition: hecmw_matrix_misc.f90:6
hecmw_mpc_prepost::hecmw_mpc_check_solution
subroutine, public hecmw_mpc_check_solution(hecMESH, hecMAT)
Definition: hecmw_mpc_prepost.f90:596
hecmw_matrix_misc::hecmw_mat_get_mpc_method
integer(kind=kint) function, public hecmw_mat_get_mpc_method(hecMAT)
Definition: hecmw_matrix_misc.f90:436
hecmw_solver_misc
Definition: hecmw_solver_misc.f90:6
hecmw_matrix_misc::hecmw_mat_finalize
subroutine, public hecmw_mat_finalize(hecMAT)
Definition: hecmw_matrix_misc.f90:204
hecmw_local_matrix::hecmw_trimatmul_ttkt_mpc
subroutine, public hecmw_trimatmul_ttkt_mpc(hecMESH, hecMAT, hecTKT)
Definition: hecmw_local_matrix.f90:842
hecmw_matrix_ass::hecmw_mat_ass_equation
subroutine, public hecmw_mat_ass_equation(hecMESH, hecMAT)
Definition: hecmw_mat_ass.f90:176
hecmw_util::hecmwst_matrix_lagrange
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...
Definition: hecmw_util_f.F90:427
hecmw_util::hecmw_comm_get_rank
integer(kind=kint) function hecmw_comm_get_rank()
Definition: hecmw_util_f.F90:582
m_hecmw_comm_f::hecmw_update_r
subroutine hecmw_update_r(hecMESH, val, n, m)
Definition: hecmw_comm_f.F90:683
hecmw_matrix_ass::hecmw_mat_ass_equation_rhs
subroutine, public hecmw_mat_ass_equation_rhs(hecMESH, hecMAT)
Definition: hecmw_mat_ass.f90:234
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444