MFIX  2016-1
leqsol_mod.f
Go to the documentation of this file.
1 MODULE leqsol
2 
3  use param, only: dim_eqs
4  use exit, only: mfix_exit
5 
6 ! Automatic adjustment of leq parameters possible (set in iterate after
7 ! the completion of first iteration).
8  LOGICAL :: leq_adjust
9 
10 ! Maximum number of linear equation solver iterations
11  INTEGER :: leq_it(dim_eqs)
12 
13 ! Linear equation solver method
14  INTEGER :: leq_method(dim_eqs)
15 
16 ! Total Iterations
17  INTEGER :: iter_tot(dim_eqs) = 0
18 
19 ! Linear equation solver sweep direction
20  CHARACTER(LEN=4) :: leq_sweep(dim_eqs)
21 
22 ! Linear equation solver tolerance
23  DOUBLE PRECISION :: leq_tol(dim_eqs)
24 
25 ! Preconditioner option
26  CHARACTER(LEN=4) :: leq_pc(dim_eqs)
27 
28 ! Option to minimize dot products
29  LOGICAL :: minimize_dotproducts
30 
31 ! Option to transpose A_m
32  LOGICAL :: do_transpose
33 
34 ! Frequency of convergence check in BiCGStab
35  INTEGER :: icheck_bicgs
36 
37 ! Optimize for massively parallel machine
38  LOGICAL :: opt_parallel
39 
40 ! Linear and non-linear solver statistics
41  LOGICAL :: solver_statistics
42 
43 CONTAINS
44 
45 !----------------------------------------------------------------------!
46 ! !
47 ! !
48 ! !
49 !----------------------------------------------------------------------!
50  SUBROUTINE report_solver_stats(TNIT, STEPS)
51 
52  use error_manager
53 
54  IMPLICIT NONE
55 
56  INTEGER, INTENT(IN) :: TNIT, STEPS
57 
58  INTEGER :: LC
59 
60  WRITE(err_msg,1100) ival(tnit), ival(tnit/steps)
61  CALL flush_err_msg(header=.false., footer=.false.)
62 
63 1100 FORMAT(/2x,'Total number of non-linear iterations: ', a,/2x,&
64  'Average number per time-step: ',a)
65 
66  WRITE(err_msg,1200)
67  CALL flush_err_msg(header=.false., footer=.false.)
68 
69 1200 FORMAT(2x,'|',10('-'),'|',13('-'),'|',14('-'),'|',/&
70  2x,'| Equation | Number of | Avg Solves |',/&
71  2x,'| Number | Solves | for NIT |',/&
72  2x,'|',10('-'),'|',13('-'),'|',14('-'),'|')
73 
74  DO lc = 1, dim_eqs
75  WRITE(err_msg,1201) lc, iter_tot(lc), iter_tot(lc)/tnit
76  CALL flush_err_msg(header=.false., footer=.false.)
77  ENDDO
78 
79 1201 FORMAT(2x,'|',3x,i3,4x,'|',2x,i9,2x,'|',2x,i10,2x,'|',/ &
80  2x,'|',10('-'),'|',13('-'),'|',14('-'),'|')
81 
82 
83  RETURN
84  END SUBROUTINE report_solver_stats
85 
86 
87 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
88 ! C
89 ! Subroutine: LEQ_MATVEC C
90 ! Purpose: Compute matrix vector multiplication C
91 ! (for linear equation Ax=b compute Ax C
92 ! C
93 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
94 ! Reviewer: Date: C
95 ! C
96 ! Literature/Document References: C
97 ! Variables referenced: C
98 ! Variables modified: C
99 ! Local variables: C
100 ! C
101 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
102 
103  SUBROUTINE leq_matvec(VNAME, VAR, A_M, Avar)
104 
105 !-----------------------------------------------
106 ! Modules
107 !-----------------------------------------------
108  USE compar, ONLY: istart, iend, jstart, jend, kstart, kend, ijkstart3, ijkend3, nlayers_bicgs, c0, c1, c2, mype
109  USE cutcell, ONLY: re_indexing, cartesian_grid
110  USE geometry, ONLY: do_k, use_corecell_loop, core_istart, core_iend, core_jstart, core_jend, core_kstart, core_kend
111  USE indices
112  USE param, ONLY: dimension_3
113  USE sendrecv, ONLY: send_recv
114  IMPLICIT NONE
115 !-----------------------------------------------
116 ! Dummy arguments
117 !-----------------------------------------------
118 ! Variable name
119  CHARACTER(LEN=*), INTENT(IN) :: Vname
120 ! Variable
121 ! DOUBLE PRECISION, INTENT(IN) :: Var(ijkstart3:ijkend3)
122  DOUBLE PRECISION, INTENT(IN) :: Var(dimension_3)
123 ! Septadiagonal matrix A_m
124 ! DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
125  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3)
126 ! Vector AVar
127 ! DOUBLE PRECISION, INTENT(OUT) :: AVar(ijkstart3:ijkend3)
128  DOUBLE PRECISION, INTENT(OUT) :: AVar(dimension_3)
129 !-----------------------------------------------
130 ! Local variables
131 !-----------------------------------------------
132 ! Variable
133  INTEGER :: I, J, K, IJK
134  integer :: im1jk, ip1jk, ijm1k, ijp1k, ijkm1, ijkp1
135  integer :: class, interval
136  integer :: j_start(2), j_end(2)
137 !-----------------------------------------------
138 
139  IF(re_indexing) THEN
140 
141  DO ijk = ijkstart3, ijkend3 ! Loop only over active cells
142 
143  im1jk = im_of(ijk)
144  ip1jk = ip_of(ijk)
145  ijm1k = jm_of(ijk)
146  ijp1k = jp_of(ijk)
147 
148  avar(ijk) = a_m(ijk,-2) * var(ijm1k) &
149  + a_m(ijk,-1) * var(im1jk) &
150  + a_m(ijk, 0) * var(ijk) &
151  + a_m(ijk, 1) * var(ip1jk) &
152  + a_m(ijk, 2) * var(ijp1k)
153 
154  if (do_k) then
155  ijkm1 = km_of(ijk)
156  ijkp1 = kp_of(ijk)
157 
158 
159  avar(ijk) = avar(ijk) + a_m(ijk,-3) * var(ijkm1) &
160  + a_m(ijk, 3) * var(ijkp1)
161 
162  endif
163 
164  enddo
165 
166  ELSE
167 
168  core_istart = istart+2
169  core_iend = iend-2
170 
171  core_jstart = jstart+2
172  core_jend = jend-2
173 
174  if (do_k) then
175  core_kstart = kstart+2
176  core_kend = kend-2
177  else
178  core_kstart = 1
179  core_kend = 1
180  kstart = 1
181  kend = 1
182  endif
183 
184  if (use_corecell_loop) then
185 
186  class = cell_class(funijk(core_istart,core_jstart,core_kstart))
187 
188 !$omp parallel do default(none) shared(c0,c1,c2,avar,a_m,var,do_k,increment_for_mp,istart,jstart,kstart,iend,jend,kend,cell_class,core_istart,core_jstart,core_kstart,core_iend,core_jend,core_kend,use_corecell_loop,class) &
189 !$omp& private(ijk,i,j,k) collapse (3)
190  do k = core_kstart,core_kend
191  do i = core_istart,core_iend
192  do j = core_jstart,core_jend
193  ijk = (j + c0 + i*c1 + k*c2)
194 
195  avar(ijk) = &
196  + a_m(ijk,-2) * var(ijk+increment_for_mp(3,class)) &
197  + a_m(ijk,-1) * var(ijk+increment_for_mp(1,class)) &
198  + a_m(ijk, 0) * var(ijk) &
199  + a_m(ijk, 1) * var(ijk+increment_for_mp(2,class)) &
200  + a_m(ijk, 2) * var(ijk+increment_for_mp(4,class))
201 
202  if (do_k) then
203  avar(ijk) = avar(ijk) + a_m(ijk,-3) * var(ijk+increment_for_mp(5,class))
204  avar(ijk) = avar(ijk) + a_m(ijk, 3) * var(ijk+increment_for_mp(6,class))
205  endif
206  enddo
207  enddo
208  enddo
209  endif
210 
211  j_start(1) = jstart
212  j_end(1) = jend
213  j_start(2) = 0 ! no iterations
214  j_end(2) = -1 ! no iterations
215 
216 !$omp parallel do default(none) shared(c0,c1,c2,avar,a_m,var,do_k,increment_for_mp,istart,jstart,kstart,iend,jend,kend,cell_class,core_istart,core_jstart,core_kstart,core_iend,core_jend,core_kend,use_corecell_loop) &
217 !$omp& private(ijk,i,j,k,class,interval) firstprivate(j_start,j_end) collapse (2)
218  do k = kstart,kend
219  do i = istart,iend
220 
221  if (use_corecell_loop) then
222  if (core_istart<= i .and. i <= core_iend .and. core_kstart <= k .and. k<=core_kend) then
223  j_start(1) = jstart
224  j_end(1) = core_jstart-1
225  j_start(2) = core_jend+1
226  j_end(2) = jend
227  else
228  j_start(1) = jstart
229  j_end(1) = jend
230  j_start(2) = 0 ! no iterations
231  j_end(2) = -1 ! no iterations
232  endif
233  endif
234 
235  do interval=1,2
236  do j = j_start(interval),j_end(interval)
237  ijk = (j + c0 + i*c1 + k*c2)
238  class = cell_class(ijk)
239 
240  avar(ijk) = &
241  + a_m(ijk,-2) * var(ijk+increment_for_mp(3,class)) &
242  + a_m(ijk,-1) * var(ijk+increment_for_mp(1,class)) &
243  + a_m(ijk, 0) * var(ijk) &
244  + a_m(ijk, 1) * var(ijk+increment_for_mp(2,class)) &
245  + a_m(ijk, 2) * var(ijk+increment_for_mp(4,class))
246 
247  if (do_k) then
248  avar(ijk) = avar(ijk) + a_m(ijk,-3) * var(ijk+increment_for_mp(5,class))
249  avar(ijk) = avar(ijk) + a_m(ijk, 3) * var(ijk+increment_for_mp(6,class))
250  endif
251  enddo
252  enddo
253  enddo
254  enddo
255 
256  ENDIF ! RE_INDEXING
257 
258  call send_recv(avar,nlayers_bicgs)
259  RETURN
260 
261  CONTAINS
262 
263  include 'functions.inc'
264 
265  END SUBROUTINE leq_matvec
266 
267 
268 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
269 ! C
270 ! Subroutine: LEQ_MSOLVE C
271 ! Purpose: C
272 ! Notes: if leq_method is biggs or cg then this subroutine is C
273 ! invoked when leq_pc='line'. if leq_method is gmres then C
274 ! this subroutine is invoked (leq_pc setting does not matter) C
275 ! C
276 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
277 ! Reviewer: Date: C
278 ! C
279 ! Literature/Document References: C
280 ! Variables referenced: C
281 ! Variables modified: C
282 ! Local variables: C
283 ! C
284 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
285 
286  SUBROUTINE leq_msolve(VNAME, B_m, A_M, Var, CMETHOD)
287 
288 !-----------------------------------------------
289 ! Modules
290 !-----------------------------------------------
291  USE param
292  USE param1
293  USE geometry
294  USE compar
295  USE indices
296  USE sendrecv
297  USE functions
298  IMPLICIT NONE
299 !-----------------------------------------------
300 ! Dummy arguments
301 !-----------------------------------------------
302 ! Variable name
303  CHARACTER(LEN=*), INTENT(IN) :: Vname
304 ! Vector b_m
305 ! DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
306  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3)
307 ! Septadiagonal matrix A_m
308 ! DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
309  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3)
310 ! Variable
311 ! DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
312  DOUBLE PRECISION, INTENT(INOUT) :: Var(dimension_3)
313 ! Sweep direction of leq solver (leq_sweep)
314 ! e.g., options = 'isis', 'rsrs' (default), 'asas'
315  CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
316 !-----------------------------------------------
317 ! Local parameters
318 !-----------------------------------------------
319  LOGICAL, PARAMETER :: USE_IKLOOP = .false.
320  LOGICAL, PARAMETER :: SETGUESS = .true.
321 !-----------------------------------------------
322 ! Local variables
323 !-----------------------------------------------
324 !
325  INTEGER :: ITER, NITER
326  INTEGER :: IJK, I , J, K
327  INTEGER :: I1, J1, K1, I2, J2, K2, IK, JK, IJ
328  INTEGER :: ISIZE, JSIZE, KSIZE
329  INTEGER :: ICASE
330 
331 ! CHARACTER(LEN=4), PARAMETER :: CMETHOD = 'II'
332  CHARACTER :: CH
333  LOGICAL :: DO_ISWEEP, DO_JSWEEP, DO_KSWEEP
334  LOGICAL :: DO_SENDRECV, DO_REDBLACK, DO_ALL
335 
336 !-----------------------------------------------
337 !!$ double precision omp_start, omp_end
338 !!$ double precision omp_get_wtime
339 ! by Tingwen
340 !!$ omp_start=omp_get_wtime()
341 
342  IF (setguess) THEN
343 
344 ! !!omp parallel do private(i,j,k,ijk) collapse(3)
345 ! do k = kstart3,kend3
346 ! do i = istart3,iend3
347 ! do j = jstart3,jend3
348 ! IJK = (J + C0 + I*C1 + K*C2)
349 ! VAR(IJK) = B_M(IJK)
350 ! enddo
351 ! enddo
352 ! enddo
353 
354  var(:) = b_m(:)
355 
356  call send_recv(var,nlayers_bicgs)
357  ENDIF
358 
359  niter = len( cmethod )
360 
361  DO iter=1,niter
362 
363 ! Perform sweeps
364  ch = cmethod( iter:iter )
365  do_isweep = (ch .EQ. 'I') .OR. (ch .EQ. 'i')
366  do_jsweep = (ch .EQ. 'J') .OR. (ch .EQ. 'j')
367  do_ksweep = (ch .EQ. 'K') .OR. (ch .EQ. 'k')
368  do_all = (ch .EQ. 'A') .OR. (ch .EQ. 'a')
369  do_redblack = (ch .EQ. 'R') .OR. (ch .EQ. 'r')
370  do_sendrecv = (ch .EQ. 'S') .OR. (ch .EQ. 's')
371 
372  IF (no_k) THEN ! two dimensional
373 ! 2D run no need to enable openmp parallel
374  IF ( do_isweep ) THEN
375 !!$omp parallel do private(I)
376  DO i=istart,iend,1
377  CALL leq_isweep( i, vname, var, a_m, b_m )
378  ENDDO
379  ENDIF
380 ! ----------------------------------------------------------------<<<
381 ! Handan Liu added 2D RSRS sweep and parallelized this loop on Jan 22 2013:
382  IF (do_redblack) THEN
383 !$omp parallel do private(I)
384  DO i=istart,iend,2
385  CALL leq_isweep( i, vname, var, a_m, b_m )
386  ENDDO
387 !$omp parallel do private(I)
388  DO i=istart+1,iend,2
389  CALL leq_isweep( i, vname, var, a_m, b_m )
390  ENDDO
391  ENDIF
392 ! ---------------------------------------------------------------->>>
393  ELSE ! three dimensional
394 
395 
396 ! do_all true only for leq_pc='asas'
397 ! ---------------------------------------------------------------->>>
398  IF(do_all) THEN ! redblack for all sweeps, not used by default
399 ! JK Loop
400 ! --------------------------------
401  j1 = jstart
402  k1 = kstart
403  j2 = jend
404  k2 = kend
405  jsize = j2-j1+1
406  ksize = k2-k1+1
407  DO icase = 1, 2
408 !!$omp parallel do private(K,J,JK)
409  DO jk=icase, ksize*jsize, 2
410  if (mod(jk,jsize).ne.0) then
411  k = int( jk/jsize ) + k1
412  else
413  k = int( jk/jsize ) + k1 -1
414  endif
415  j = (jk-1-(k-k1)*jsize) + j1 + mod(k,2)
416  if(j.gt.j2) j=j-j2 + j1 -1
417  CALL leq_jksweep(j, k, vname, var, a_m, b_m)
418  ENDDO
419 
420  ENDDO
421  call send_recv(var,nlayers_bicgs)
422 
423 ! IJ Loop
424 ! --------------------------------
425  i1 = istart
426  j1 = jstart
427  i2 = iend
428  j2 = jend
429  isize = i2-i1+1
430  jsize = j2-j1+1
431  DO icase = 1, 2
432 !!$omp parallel do private(J,I,IJ)
433  DO ij=icase, jsize*isize, 2
434  if (mod(ij,isize).ne.0) then
435  j = int( ij/isize ) + j1
436  else
437  j = int( ij/isize ) + j1 -1
438  endif
439  i = (ij-1-(j-j1)*isize) + i1 + mod(j,2)
440  if(i.gt.i2) i=i-i2 + i1 -1
441  CALL leq_ijsweep(i, j, vname, var, a_m, b_m)
442  ENDDO
443 
444  ENDDO
445  call send_recv(var,nlayers_bicgs)
446 
447 ! IK Loop
448 ! --------------------------------
449  i1 = istart
450  k1 = kstart
451  i2 = iend
452  k2 = kend
453  isize = i2-i1+1
454  ksize = k2-k1+1
455 
456  DO icase = 1, 2
457 !!$omp parallel do private(K,I,IK)
458  DO ik=icase, ksize*isize, 2
459  if (mod(ik,isize).ne.0) then
460  k = int( ik/isize ) + k1
461  else
462  k = int( ik/isize ) + k1 -1
463  endif
464  i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
465  if(i.gt.i2) i=i-i2 + i1 -1
466  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
467  ENDDO
468 
469  ENDDO
470  ENDIF ! end DO_ALL
471 ! ----------------------------------------------------------------<<<
472 
473 ! do_redblack only true leq_pc='rsrs'
474 ! ---------------------------------------------------------------->>>
475  IF(do_redblack) THEN
476 !i1 = istart
477 !k1 = kstart
478 !i2 = iend
479 !k2 = kend
480 !isize = i2-i1+1
481 !ksize = k2-k1+1
482 ! DO icase = 1, 2
483 !!$omp parallel do private(K,I,IK)
484 ! DO IK=icase, ksize*isize, 2
485 ! if (mod(ik,isize).ne.0) then
486 ! k = int( ik/isize ) + k1
487 ! else
488 ! k = int( ik/isize ) + k1 -1
489 ! endif
490 ! i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
491 ! if(i.gt.i2) i=i-i2 + i1 -1
492 ! CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
493 ! ENDDO
494 ! ENDDO
495 ! ELSE
496 ! Handan Liu split above loop for OpenMP at May 22 2013, modified at July 17
497 !$omp parallel do default(shared) private(I,K) schedule(auto)
498  DO k=kstart,kend
499  IF(mod(k,2).ne.0)THEN
500  DO i=istart+1,iend,2
501  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
502  ENDDO
503  ELSE
504  DO i=istart,iend,2
505  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
506  ENDDO
507  ENDIF
508  ENDDO
509 !$omp end parallel do
510 !$omp parallel do default(shared) private(I,K) schedule(auto)
511  DO k=kstart,kend
512  IF(mod(k,2).ne.0)THEN
513  DO i=istart,iend,2
514  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
515  ENDDO
516  ELSE
517  DO i=istart+1,iend,2
518  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
519  ENDDO
520  ENDIF
521  ENDDO
522 !$omp end parallel do
523 
524  ENDIF ! end if(do_redblack)
525 ! ----------------------------------------------------------------<<<
526 
527 ! Not sure the purpose of us_ikloop
528 ! The SMP directives below need review !Tingwen Jan 2012
529  IF(use_ikloop) THEN
530 ! use_ikloop is currently hard-wired to false (so goto else branch)
531 ! ---------------------------------------------------------------->>>
532  i1 = istart
533  k1 = kstart
534  i2 = iend
535  k2 = kend
536  isize = i2-i1+1
537  ksize = k2-k1+1
538  IF (do_isweep) THEN
539 !!$omp parallel do private(K,I,IK)
540  DO ik=1, ksize*isize
541  if (mod(ik,isize).ne.0) then
542  k = int( ik/isize ) + k1
543  else
544  k = int( ik/isize ) + k1 -1
545  endif
546  i = (ik-1-(k-k1)*isize) + i1
547  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
548  ENDDO
549  ENDIF
550  IF (do_ksweep) THEN
551 !!$omp parallel do private(K,I,IK)
552  DO ik=1, ksize*isize
553  if (mod(ik,ksize).ne.0) then
554  i = int( ik/ksize ) + i1
555  else
556  i = int( ik/ksize ) + i1 -1
557  endif
558  k = (ik-1-(i-i1)*ksize) + k1
559  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
560  ENDDO
561  ENDIF
562 ! ----------------------------------------------------------------<<<
563  ELSE ! else branch of if(use_ikloop)
564 ! Not sure the purpose of us_ikloop
565 ! The SMP directives below need review !Tingwen Jan 2012
566 ! ---------------------------------------------------------------->>>
567  IF (do_isweep) THEN
568 !!$omp parallel do private(K,I)
569  DO k=kstart,kend
570  DO i=istart,iend
571  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
572  ENDDO
573  ENDDO
574  ENDIF
575  IF (do_ksweep) THEN
576 !!$omp parallel do private(K,I)
577  DO i=istart,iend
578  DO k=kstart,kend
579  CALL leq_iksweep(i, k, vname, var, a_m, b_m)
580  ENDDO
581  ENDDO
582  ENDIF
583  ENDIF ! end if/else (use(ikloop)
584 ! ----------------------------------------------------------------<<<
585 
586  ENDIF ! end if/else (do_k)
587 
588 
589 ! this is called for all settings of leq_pc
590  IF (do_sendrecv) call send_recv(var,nlayers_bicgs)
591 
592 
593  ENDDO ! end do iter=1,niter
594 !!$ omp_end=omp_get_wtime()
595 !!$ write(*,*)'leq_msolve:',omp_end - omp_start
596 
597  RETURN
598  END SUBROUTINE leq_msolve
599 
600 
601 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
602 ! C
603 ! Subroutine: LEQ_MSOLVE0 C
604 ! Notes: do nothing or no preconditioning (leq_pc='none') C
605 ! C
606 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
607 ! Reviewer: Date: C
608 ! C
609 ! Literature/Document References: C
610 ! Variables referenced: C
611 ! Variables modified: C
612 ! Local variables: C
613 ! C
614 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
615 
616  SUBROUTINE leq_msolve0(VNAME, B_m, A_M, Var, CMETHOD)
617 
618 !-----------------------------------------------
619 ! Modules
620 !-----------------------------------------------
621  USE param
622  USE param1
623  USE geometry
624  USE compar
625  USE indices
626  USE sendrecv
627  use parallel
628  IMPLICIT NONE
629 !-----------------------------------------------
630 ! Dummy arguments
631 !-----------------------------------------------
632 ! Variable name
633  CHARACTER(LEN=*), INTENT(IN) :: Vname
634 ! Vector b_m
635 ! DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
636  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3)
637 ! Septadiagonal matrix A_m
638 ! DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
639  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3)
640 ! Variable
641 ! DOUBLE PRECISION, INTENT(OUT) :: Var(ijkstart3:ijkend3)
642  DOUBLE PRECISION, INTENT(OUT) :: Var(dimension_3)
643 ! sweep direction
644  CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
645 !-----------------------------------------------
646 ! Local variables
647 !-----------------------------------------------
648  integer :: ijk
649 !-----------------------------------------------
650 
651 ! do nothing or no preconditioning
652  if (use_doloop) then ! mfix.dat keyword default=false
653 !!$omp parallel do private(ijk)
654  do ijk=ijkstart3,ijkend3
655  var(ijk) = b_m(ijk)
656  enddo
657  else
658  var(:) = b_m(:)
659  endif
660  call send_recv(var,nlayers_bicgs)
661 
662  return
663  end subroutine leq_msolve0
664 
665 
666 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
667 ! C
668 ! Subroutine: LEQ_MSOLVE1 C
669 ! Notes: diagonal scaling (leq_pc='diag') C
670 ! C
671 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
672 ! Reviewer: Date: C
673 ! C
674 ! Literature/Document References: C
675 ! Variables referenced: C
676 ! Variables modified: C
677 ! Local variables: C
678 ! C
679 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
680 
681  SUBROUTINE leq_msolve1(VNAME, B_m, A_M, Var, CMETHOD)
682 
683 !-----------------------------------------------
684 ! Modules
685 !-----------------------------------------------
686  USE param
687  USE param1
688  USE geometry
689  USE compar
690  USE indices
691  USE sendrecv
692  use parallel
693 ! USE cutcell, only: RE_INDEXING,INTERIOR_CELL_AT
694  USE cutcell
695  USE functions
696  IMPLICIT NONE
697 !-----------------------------------------------
698 ! Dummy arguments
699 !-----------------------------------------------
700 ! Variable name
701  CHARACTER(LEN=*), INTENT(IN) :: Vname
702 ! Vector b_m
703 ! DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
704  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3)
705 ! Septadiagonal matrix A_m
706 ! DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
707  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3)
708 ! Variable
709 ! DOUBLE PRECISION, INTENT(OUT) :: Var(ijkstart3:ijkend3)
710  DOUBLE PRECISION, INTENT(OUT) :: Var(dimension_3)
711 ! sweep direction
712  CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
713 !-----------------------------------------------
714 ! Local variables
715 !-----------------------------------------------
716  integer :: i,j,k, ijk
717 !-----------------------------------------------
718 
719  if (use_doloop) then ! mfix.dat keyword default=false
720 !!$omp parallel do private(ijk)
721  do ijk=ijkstart3,ijkend3
722  var(ijk) = zero
723  enddo
724  else
725  var(:) = zero
726  endif
727 
728 ! diagonal scaling
729  IF(.NOT.re_indexing) THEN
730 !$omp parallel do private(i,j,k,ijk) collapse (3)
731  do k=kstart2,kend2
732  do i=istart2,iend2
733  do j=jstart2,jend2
734  ijk = funijk( i,j,k )
735  var(ijk) = b_m(ijk)/a_m(ijk,0)
736  enddo
737  enddo
738  enddo
739  ELSE
740 !$omp parallel do private(ijk) collapse (1)
741  DO ijk=ijkstart3,ijkend3
742  var(ijk) = b_m(ijk)/a_m(ijk,0)
743  ENDDO
744  ENDIF
745 
746  call send_recv(var,nlayers_bicgs)
747 
748  return
749  end subroutine leq_msolve1
750 
751 
752 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
753 ! C
754 ! Subroutine: LEQ_ISWEEP C
755 ! Purpose: Perform line sweep at coordinate I C
756 ! C
757 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
758 ! Reviewer: Date: C
759 ! C
760 ! Literature/Document References: C
761 ! Variables referenced: C
762 ! Variables modified: C
763 ! Local variables: C
764 ! C
765 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
766 
767  SUBROUTINE leq_isweep(I, Vname, VAR, A_M, B_M)
768 
769 !-----------------------------------------------
770 ! Modules
771 !-----------------------------------------------
772  USE param
773  USE param1
774  USE geometry
775  USE compar
776  USE indices
777  USE funits
778  USE sendrecv
779  USE mpi_utility
780  USE functions
781  IMPLICIT NONE
782 !-----------------------------------------------
783 ! Dummy arguments
784 !-----------------------------------------------
785 ! Line position
786  INTEGER, INTENT(IN) :: I
787 ! Variable name
788  CHARACTER(LEN=*), INTENT(IN) :: Vname
789 ! Variable
790  DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
791 ! Septadiagonal matrix A_m
792  DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
793 ! Vector b_m
794  DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
795 !-----------------------------------------------
796 ! Local variables
797 !-----------------------------------------------
798  DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC, DD, EE, BB
799  INTEGER :: NSTART, NEND, INFO
800  INTEGER :: IJK, J, K, IM1JK, IP1JK
801 !-----------------------------------------------
802 
803  nend = jend
804  nstart = jstart
805  k = 1
806 
807  DO j=nstart, nend
808  ijk = funijk(i,j,k)
809  im1jk = im_of(ijk)
810  ip1jk = ip_of(ijk)
811  dd(j) = a_m(ijk, 0)
812  cc(j) = a_m(ijk, -2)
813  ee(j) = a_m(ijk, 2)
814  bb(j) = b_m(ijk) - a_m(ijk,-1) * var( im1jk ) &
815  - a_m(ijk, 1) * var( ip1jk )
816  ENDDO
817 
818  cc(nstart) = zero
819  ee(nend) = zero
820  info = 0
821 ! CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
822  CALL dgtsv( jend-jstart+1, 1, cc(jstart+1), dd, ee, bb, jend-jstart+1, info )
823 
824  IF (info.NE.0) THEN
825  RETURN
826  ENDIF
827 
828  DO j=nstart, nend
829  ijk = funijk(i,j,k)
830  var(ijk) = bb(j)
831  ENDDO
832 
833  RETURN
834  END SUBROUTINE leq_isweep
835 
836 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
837 ! C
838 ! Subroutine: LEQ_IKSWEEP C
839 ! Purpose: Perform line sweep at coordinate I, K C
840 ! C
841 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
842 ! Reviewer: Date: C
843 ! C
844 ! Literature/Document References: C
845 ! Variables referenced: C
846 ! Variables modified: C
847 ! Local variables: C
848 ! C
849 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
850 
851  SUBROUTINE leq_iksweep(I, K, Vname, VAR, A_M, B_M)
852 
853 !-----------------------------------------------
854 ! Modules
855 !-----------------------------------------------
856  USE param
857  USE param1
858  USE geometry
859  USE compar
860  USE funits
861  USE indices
862  USE sendrecv
863  USE mpi_utility
864  USE functions
865  IMPLICIT NONE
866 !-----------------------------------------------
867 ! Dummy arguments
868 !-----------------------------------------------
869 ! Line position
870  INTEGER, INTENT(IN) :: I, K
871 ! Variable name
872  CHARACTER(LEN=*), INTENT(IN) :: Vname
873 ! Variable
874  DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
875 ! Septadiagonal matrix A_m
876  DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
877 ! Vector b_m
878  DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
879 !-----------------------------------------------
880 ! Local variables
881 !-----------------------------------------------
882  DOUBLE PRECISION, DIMENSION(JSTART:JEND) :: CC, DD, EE, BB
883  INTEGER :: NSTART, NEND, INFO
884  INTEGER :: IJK, J, CLASS
885 !-----------------------------------------------
886 
887  nend = jend
888  nstart = jstart
889 
890 !!$omp parallel do private(j,ijk,im1jk,ip1jk,ijkm1,ijkp1)
891  DO j=nstart, nend
892 ! IJK = FUNIJK(IMAP_C(I),JMAP_C(J),KMAP_C(K))
893  ijk = (j + c0 + i*c1 + k*c2)
894  CLASS = cell_class(ijk)
895  dd(j) = a_m(ijk, 0)
896  cc(j) = a_m(ijk, -2)
897  ee(j) = a_m(ijk, 2)
898  bb(j) = b_m(ijk) - a_m(ijk,-1) * var( ijk+increment_for_mp(1,class) ) &
899  - a_m(ijk, 1) * var( ijk+increment_for_mp(2,class) ) &
900  - a_m(ijk,-3) * var( ijk+increment_for_mp(5,class) ) &
901  - a_m(ijk, 3) * var( ijk+increment_for_mp(6,class) )
902  ENDDO
903 
904  cc(nstart) = zero
905  ee(nend) = zero
906  info = 0
907 ! CALL DGTSL( JEND-JSTART+1, CC, DD, EE, BB, INFO )
908  CALL dgtsv(nend-nstart+1, 1, cc(nstart+1), dd, ee, bb, nend-nstart+1, info)
909 
910  IF (info.NE.0) THEN
911  write(*,*) 'leq_iksweep',info, mype
912  IF(dmp_log)WRITE (unit_log,*) 'ROUTINE = ', ' IKSWEEP'
913  RETURN
914  ENDIF
915 
916  DO j=nstart, nend
917  var(j + c0 + i*c1 + k*c2) = bb(j)
918  ENDDO
919 
920  RETURN
921  END SUBROUTINE leq_iksweep
922 
923 
924 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
925 ! C
926 ! Subroutine: LEQ_JKSWEEP C
927 ! Purpose: Perform line sweep at coordinate I, K C
928 ! C
929 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
930 ! Reviewer: Date: C
931 ! C
932 ! Literature/Document References: C
933 ! Variables referenced: C
934 ! Variables modified: C
935 ! Local variables: C
936 ! C
937 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
938 
939  SUBROUTINE leq_jksweep(J, K, Vname, VAR, A_M, B_M)
940 
941 !-----------------------------------------------
942 ! Modules
943 !-----------------------------------------------
944  USE param
945  USE param1
946  USE geometry
947  USE funits
948  USE compar
949  USE indices
950  USE functions
951  IMPLICIT NONE
952 !-----------------------------------------------
953 ! Dummy arguments
954 !-----------------------------------------------
955 ! Line position
956  INTEGER, INTENT(IN) :: J, K
957 ! Variable name
958  CHARACTER(LEN=*), INTENT(IN) :: Vname
959 ! Variable
960  DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
961 ! Septadiagonal matrix A_m
962  DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
963 ! Vector b_m
964  DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
965 !-----------------------------------------------
966 ! Local variables
967 !-----------------------------------------------
968  DOUBLE PRECISION, DIMENSION (ISTART:IEND) :: CC, DD, EE, BB
969  INTEGER :: NSTART, NEND, INFO, IJK, I
970 !-----------------------------------------------
971 
972  nend = iend
973  nstart = istart
974 
975  DO i=nstart,nend
976  ijk = funijk(i,j,k)
977  dd(i) = a_m(ijk, 0)
978  cc(i) = a_m(ijk, -1)
979  ee(i) = a_m(ijk, 1)
980  bb(i) = b_m(ijk) - a_m(ijk,-2) * var( jm_of(ijk) ) &
981  - a_m(ijk, 2) * var( jp_of(ijk) ) &
982  - a_m(ijk,-3) * var( km_of(ijk) ) &
983  - a_m(ijk, 3) * var( kp_of(ijk) )
984  ENDDO
985 
986  cc(nstart) = zero
987  ee(nend) = zero
988  info = 0
989  CALL dgtsv(nend-nstart+1, 1, cc(nstart+1), dd, ee, bb, nend-nstart+1, info)
990 
991  IF (info.NE.0) THEN
992  IF(dmp_log)WRITE (unit_log,*) 'VNAME = ', vname
993  IF(dmp_log)WRITE (unit_log,*) 'ROUTINE = ', ' JKSWEEP'
994  IF(dmp_log)WRITE (unit_log,*) 'DGTSV RETURNS INFO = ', info
995  call mfix_exit(mype)
996  ENDIF
997 
998  DO i=nstart, nend
999  ijk = funijk(i,j,k)
1000  var(ijk) = bb(i)
1001  ENDDO
1002 
1003  RETURN
1004  END SUBROUTINE leq_jksweep
1005 
1006 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1007 ! C
1008 ! Subroutine: LEQ_IJSWEEP C
1009 ! Purpose: Perform line sweep at coordinate I, K C
1010 ! C
1011 ! Author: Ed D'Azevedo Date: 21-JAN-99 C
1012 ! Reviewer: Date: C
1013 ! C
1014 ! Literature/Document References: C
1015 ! Variables referenced: C
1016 ! Variables modified: C
1017 ! Local variables: C
1018 ! C
1019 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1020 
1021  SUBROUTINE leq_ijsweep(I, J, Vname, VAR, A_M, B_M)
1022 
1023 !-----------------------------------------------
1024 ! Modules
1025 !-----------------------------------------------
1026  USE param
1027  USE param1
1028  USE geometry
1029  USE funits
1030  USE compar
1031  USE indices
1032  USE functions
1033  IMPLICIT NONE
1034 !-----------------------------------------------
1035 ! Dummy arguments
1036 !-----------------------------------------------
1037 ! Line position
1038  INTEGER, INTENT(IN) :: I, J
1039 ! Variable name
1040  CHARACTER(LEN=*), INTENT(IN) :: Vname
1041 ! Variable
1042  DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
1043 ! Septadiagonal matrix A_m
1044  DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
1045 ! Vector b_m
1046  DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
1047 !-----------------------------------------------
1048 ! Local variables
1049 !-----------------------------------------------
1050  DOUBLE PRECISION, DIMENSION (KSTART:KEND) :: CC, DD, EE, BB
1051  INTEGER :: NEND, NSTART, INFO, IJK, K
1052 !-----------------------------------------------
1053 
1054  nend = kend
1055  nstart = kstart
1056 
1057  DO k=nstart, nend
1058  ijk = funijk(i,j,k)
1059  dd(k) = a_m(ijk, 0)
1060  cc(k) = a_m(ijk, -3)
1061  ee(k) = a_m(ijk, 3)
1062  bb(k) = b_m(ijk) - a_m(ijk,-2) * var( jm_of(ijk) ) &
1063  - a_m(ijk, 2) * var( jp_of(ijk) ) &
1064  - a_m(ijk,-1) * var( im_of(ijk) ) &
1065  - a_m(ijk, 1) * var( ip_of(ijk) )
1066  ENDDO
1067 
1068  cc(nstart) = zero
1069  ee(nend) = zero
1070  info = 0
1071  CALL dgtsv(nend-nstart+1, 1, cc(nstart+1), dd, ee, bb, nend-nstart+1, info)
1072 
1073  IF (info.NE.0) THEN
1074  IF(dmp_log)WRITE (unit_log,*) 'VNAME = ', vname
1075  IF(dmp_log)WRITE (unit_log,*) 'ROUTINE = ', ' IJSWEEP'
1076  IF(dmp_log)WRITE (unit_log,*) 'DGTSV RETURNS INFO = ', info
1077  call mfix_exit(mype)
1078  ENDIF
1079 
1080  DO k=nstart, nend
1081  ijk = funijk(i,j,k)
1082  var(ijk) = bb(k)
1083  ENDDO
1084 
1085  RETURN
1086  END SUBROUTINE leq_ijsweep
1087 
1088 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1089 ! C
1090 ! C
1091 ! C
1092 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1093 
1094  double precision function dot_product_par(r1,r2)
1095 
1096 !-----------------------------------------------
1097 ! Modules
1098 !-----------------------------------------------
1099  use compar
1100  use cutcell
1101  use functions
1102  use geometry
1103  use indices
1104  use parallel, only: is_serial
1105  use mpi_utility
1106  use param, only: dimension_3
1107  implicit none
1108 !-----------------------------------------------
1109 ! Dummy arguments
1110 !-----------------------------------------------
1111 ! double precision, intent(in), dimension(ijkstart3:ijkend3) :: r1,r2
1112  double precision, intent(in), dimension(DIMENSION_3) :: r1,r2
1113 !-----------------------------------------------
1114 ! Local parameters
1115 !-----------------------------------------------
1116  logical, parameter :: do_global_sum = .true.
1117 !-----------------------------------------------
1118 ! Local variables
1119 !-----------------------------------------------
1120  DOUBLE PRECISION, allocatable, Dimension(:) :: r1_g, r2_g
1121  double precision :: prod
1122  integer :: i, j, k, ijk
1123 !-----------------------------------------------
1124 
1125  if (numpes.eq.1.and.is_serial) then
1126 !!$ if (.false.) then
1127  dot_product_par = dot_product(r1,r2)
1128 !!$ endif
1129 
1130 !!$ dot_product_par = 0
1131 !!$omp parallel do private(ijk) reduction(+:dot_product_par)
1132 !!$ do ijk = 1, DIMENSION_3
1133 !!$ dot_product_par = dot_product_par + r1(ijk)*r2(ijk)
1134 !!$ enddo
1135  return
1136  endif
1137 
1138  if(do_global_sum) then
1139  prod = 0.0d0
1140 
1141  IF(re_indexing) THEN
1142 ! IF(.FALSE.) THEN
1143 ! Somehow, looping in this order leads to smaller time step than k,i,j nested loop below ....
1144  DO ijk = ijkstart3,ijkend3
1145  IF(interior_cell_at(ijk)) prod = prod + r1(ijk)*r2(ijk)
1146  ENDDO
1147 
1148  call global_all_sum(prod, dot_product_par)
1149 
1150  ELSE
1151 
1152 !$omp parallel do private(i,j,k,ijk) reduction(+:prod) collapse (3)
1153  do k = kstart1, kend1
1154  do i = istart1, iend1
1155  do j = jstart1, jend1
1156  ijk = funijk_map_c(i,j,k)
1157  prod = prod + r1(ijk)*r2(ijk)
1158  enddo
1159  enddo
1160  enddo
1161 
1162  call global_all_sum(prod, dot_product_par)
1163 
1164  ENDIF
1165 
1166  else
1167  if(mype.eq.root) then
1168  allocate (r1_g(1:ijkmax3))
1169  allocate (r2_g(1:ijkmax3))
1170  else
1171  allocate (r1_g(10))
1172  allocate (r2_g(10))
1173  endif
1174  call gather(r1,r1_g)
1175  call gather(r2,r2_g)
1176 
1177  if(mype.eq.root) then
1178  prod = 0.0d0
1179 
1180 !$omp parallel do private(i,j,k,ijk) reduction(+:prod) collapse (3)
1181  do k = kmin1, kmax1
1182  do i = imin1, imax1
1183  do j = jmin1, jmax1
1184  ijk = funijk_gl(imap_c(i),jmap_c(j),kmap_c(k))
1185 ! ijk = funijk_gl (i,j,k)
1186  prod = prod + r1_g(ijk)*r2_g(ijk)
1187  enddo
1188  enddo
1189  enddo
1190 
1191  endif
1192  call bcast( prod)
1193 
1194  dot_product_par = prod
1195 
1196  deallocate (r1_g)
1197  deallocate (r2_g)
1198 
1199  endif
1200 
1201  end function dot_product_par
1202 
1203 
1204 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1205 ! C
1206 ! C
1207 ! C
1208 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1209 
1210  function dot_product_par2(r1,r2,r3,r4)
1211 
1212 !-----------------------------------------------
1213 ! Modules
1214 !-----------------------------------------------
1215  use mpi_utility
1216  use geometry
1217  use compar
1218  use indices
1219  use functions
1220  implicit none
1221 !-----------------------------------------------
1222 ! Dummy arguments
1223 !-----------------------------------------------
1224  double precision, intent(in), dimension(ijkstart3:ijkend3) :: r1,r2,r3,r4
1225 !-----------------------------------------------
1226 ! Local parameters
1227 !-----------------------------------------------
1228  logical, parameter :: do_global_sum = .true.
1229 !-----------------------------------------------
1230 ! Local variables
1231 !-----------------------------------------------
1232  DOUBLE PRECISION, allocatable, Dimension(:,:) :: r_temp, rg_temp
1233  double precision, Dimension(2) :: prod, dot_product_par2
1234  integer :: i, j, k, ijk
1235 !-----------------------------------------------
1236 
1237  if(do_global_sum) then
1238 
1239  prod(:) = 0.0d0
1240 
1241 !$omp parallel do private(i,j,k,ijk) reduction(+:prod) collapse (3)
1242  do k = kstart1, kend1
1243  do i = istart1, iend1
1244  do j = jstart1, jend1
1245 
1246  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
1247 
1248  ijk = funijk_map_c(i,j,k)
1249  prod(1) = prod(1) + r1(ijk)*r2(ijk)
1250  prod(2) = prod(2) + r3(ijk)*r4(ijk)
1251  enddo
1252  enddo
1253  enddo
1254 
1255  call global_all_sum(prod, dot_product_par2)
1256 
1257  else
1258  allocate (r_temp(ijkstart3:ijkend3,4))
1259  r_temp(:,1) = r1
1260  r_temp(:,2) = r2
1261  r_temp(:,3) = r3
1262  r_temp(:,4) = r4
1263 
1264  if(mype.eq.root) then
1265  allocate (rg_temp(1:ijkmax3,4))
1266  else
1267  allocate (rg_temp(10,4))
1268  endif
1269  call gather(r_temp,rg_temp)
1270 
1271  if(mype.eq.root) then
1272  prod = 0.0d0
1273 !$omp parallel do private(i,j,k,ijk) reduction(+:prod) collapse (3)
1274  do k = kmin1, kmax1
1275  do i = imin1, imax1
1276  do j = jmin1, jmax1
1277  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
1278  ijk = funijk_gl(imap_c(i),jmap_c(j),kmap_c(k))
1279 ! ijk = funijk_gl (i,j,k)
1280  prod(1) = prod(1) + rg_temp(ijk,1)*rg_temp(ijk,2)
1281  prod(2) = prod(2) + rg_temp(ijk,3)*rg_temp(ijk,4)
1282  enddo
1283  enddo
1284  enddo
1285  endif
1286  call bcast( prod)
1287 
1288  dot_product_par2 = prod
1289 
1290  deallocate (r_temp)
1291  deallocate (rg_temp)
1292 
1293  endif
1294 
1295  end function dot_product_par2
1296 
1297 END MODULE leqsol
compar::jend2
integer jend2
Definition: compar_mod.f:80
mpi_utility
Definition: mpi_utility_mod.f:6
compar::jmap_c
integer, dimension(:), allocatable jmap_c
Definition: compar_mod.f:78
cutcell::re_indexing
logical re_indexing
Definition: cutcell_mod.f:16
funits::dmp_log
logical dmp_log
Definition: funits_mod.f:6
compar::kmap_c
integer, dimension(:), allocatable kmap_c
Definition: compar_mod.f:78
compar::nlayers_bicgs
integer nlayers_bicgs
Definition: compar_mod.f:45
compar::c0
integer c0
Definition: compar_mod.f:104
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
compar::kend1
integer kend1
Definition: compar_mod.f:80
compar::istart1
integer istart1
Definition: compar_mod.f:80
leqsol::leq_adjust
logical leq_adjust
Definition: leqsol_mod.f:8
leqsol::leq_matvec
subroutine leq_matvec(VNAME, VAR, A_M, Avar)
Definition: leqsol_mod.f:104
leqsol::leq_iksweep
subroutine leq_iksweep(I, K, Vname, VAR, A_M, B_M)
Definition: leqsol_mod.f:852
compar::iend1
integer iend1
Definition: compar_mod.f:80
leqsol::icheck_bicgs
integer icheck_bicgs
Definition: leqsol_mod.f:35
compar::iend
integer iend
Definition: compar_mod.f:87
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
geometry::core_jend
integer core_jend
Definition: geometry_mod.f:249
compar::istart2
integer istart2
Definition: compar_mod.f:80
geometry::core_jstart
integer core_jstart
Definition: geometry_mod.f:249
leqsol::iter_tot
integer, dimension(dim_eqs) iter_tot
Definition: leqsol_mod.f:17
compar::iend2
integer iend2
Definition: compar_mod.f:80
error_manager::ival
Definition: error_manager_mod.f:17
param::dim_eqs
integer, parameter dim_eqs
Definition: param_mod.f:97
leqsol::leq_isweep
subroutine leq_isweep(I, Vname, VAR, A_M, B_M)
Definition: leqsol_mod.f:768
leqsol::leq_msolve
subroutine leq_msolve(VNAME, B_m, A_M, Var, CMETHOD)
Definition: leqsol_mod.f:287
geometry::core_kend
integer core_kend
Definition: geometry_mod.f:250
leqsol::leq_sweep
character(len=4), dimension(dim_eqs) leq_sweep
Definition: leqsol_mod.f:20
compar::kend2
integer kend2
Definition: compar_mod.f:80
leqsol::do_transpose
logical do_transpose
Definition: leqsol_mod.f:32
compar::kstart2
integer kstart2
Definition: compar_mod.f:80
compar::kstart
integer kstart
Definition: compar_mod.f:87
compar::kstart1
integer kstart1
Definition: compar_mod.f:80
mpi_utility::gather
Definition: mpi_utility_mod.f:46
compar::numpes
integer numpes
Definition: compar_mod.f:24
leqsol::leq_msolve0
subroutine leq_msolve0(VNAME, B_m, A_M, Var, CMETHOD)
Definition: leqsol_mod.f:617
geometry::ijkmax3
integer ijkmax3
Definition: geometry_mod.f:82
leqsol::leq_jksweep
subroutine leq_jksweep(J, K, Vname, VAR, A_M, B_M)
Definition: leqsol_mod.f:940
geometry::core_kstart
integer core_kstart
Definition: geometry_mod.f:250
sendrecv::send_recv
Definition: sendrecv_mod.f:75
compar::c2
integer c2
Definition: compar_mod.f:104
geometry::kmax1
integer kmax1
Definition: geometry_mod.f:58
leqsol::leq_ijsweep
subroutine leq_ijsweep(I, J, Vname, VAR, A_M, B_M)
Definition: leqsol_mod.f:1022
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
leqsol::leq_it
integer, dimension(dim_eqs) leq_it
Definition: leqsol_mod.f:11
compar::dead_cell_at
logical, dimension(:,:,:), allocatable dead_cell_at
Definition: compar_mod.f:127
geometry::imax1
integer imax1
Definition: geometry_mod.f:54
leqsol::dot_product_par
double precision function dot_product_par(r1, r2)
Definition: leqsol_mod.f:1095
compar::root
integer root
Definition: compar_mod.f:41
geometry::core_iend
integer core_iend
Definition: geometry_mod.f:248
exit
Definition: exit.f:2
compar::jstart2
integer jstart2
Definition: compar_mod.f:80
indices::increment_for_mp
integer, dimension(6, max_class) increment_for_mp
Definition: indices_mod.f:27
leqsol::dot_product_par2
double precision function, dimension(2) dot_product_par2(r1, r2, r3, r4)
Definition: leqsol_mod.f:1211
dgtsv
subroutine dgtsv(N, NRHS, DL, D, DU, B, LDB, INFO)
Definition: DGTSV.f:3
parallel::is_serial
logical is_serial
Definition: parallel_mod.f:11
mpi_utility::global_all_sum
Definition: mpi_utility_mod.f:66
leqsol::leq_tol
double precision, dimension(dim_eqs) leq_tol
Definition: leqsol_mod.f:23
geometry::jmax1
integer jmax1
Definition: geometry_mod.f:56
funits::unit_log
integer, parameter unit_log
Definition: funits_mod.f:21
compar::imap_c
integer, dimension(:), allocatable imap_c
Definition: compar_mod.f:78
leqsol::report_solver_stats
subroutine report_solver_stats(TNIT, STEPS)
Definition: leqsol_mod.f:51
param
Definition: param_mod.f:2
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
compar::kend
integer kend
Definition: compar_mod.f:87
geometry::no_k
logical no_k
Definition: geometry_mod.f:28
leqsol::leq_msolve1
subroutine leq_msolve1(VNAME, B_m, A_M, Var, CMETHOD)
Definition: leqsol_mod.f:682
leqsol::minimize_dotproducts
logical minimize_dotproducts
Definition: leqsol_mod.f:29
geometry::use_corecell_loop
logical use_corecell_loop
Definition: geometry_mod.f:246
cutcell::interior_cell_at
logical, dimension(:), allocatable interior_cell_at
Definition: cutcell_mod.f:40
geometry::jmin1
integer jmin1
Definition: geometry_mod.f:42
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
compar::jstart
integer jstart
Definition: compar_mod.f:87
compar::mype
integer mype
Definition: compar_mod.f:24
parallel::use_doloop
logical use_doloop
Definition: parallel_mod.f:10
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
error_manager::err_msg
character(len=line_length), dimension(line_count) err_msg
Definition: error_manager_mod.f:29
compar::istart
integer istart
Definition: compar_mod.f:87
compar::jend
integer jend
Definition: compar_mod.f:87
error_manager
Definition: error_manager_mod.f:8
leqsol::leq_method
integer, dimension(dim_eqs) leq_method
Definition: leqsol_mod.f:14
mpi_utility::bcast
Definition: mpi_utility_mod.f:53
geometry::imin1
integer imin1
Definition: geometry_mod.f:40
geometry::kmin1
integer kmin1
Definition: geometry_mod.f:44
leqsol::solver_statistics
logical solver_statistics
Definition: leqsol_mod.f:41
leqsol::opt_parallel
logical opt_parallel
Definition: leqsol_mod.f:38
geometry::core_istart
integer core_istart
Definition: geometry_mod.f:248
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
compar::c1
integer c1
Definition: compar_mod.f:104
compar::jend1
integer jend1
Definition: compar_mod.f:80
error_manager::flush_err_msg
subroutine flush_err_msg(DEBUG, HEADER, FOOTER, ABORT, LOG, CALL_TREE)
Definition: error_manager_mod.f:305
compar::jstart1
integer jstart1
Definition: compar_mod.f:80
indices::cell_class
integer, dimension(:), allocatable cell_class
Definition: indices_mod.f:42
leqsol::leq_pc
character(len=4), dimension(dim_eqs) leq_pc
Definition: leqsol_mod.f:26