File: /nfs/home/0/users/jenkins/mfix.git/model/leq_bicgs.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24 SUBROUTINE LEQ_BICGS(VNAME, VNO, VAR, A_M, B_m, cmethod, &
25 TOL, PC, ITMAX, IER)
26
27
28
29
30 USE param
31 USE param1
32 USE matrix
33 USE geometry
34 USE compar
35 USE indices
36 USE leqsol
37 USE funits
38 IMPLICIT NONE
39
40
41
42
43 CHARACTER(LEN=*), INTENT(IN) :: Vname
44
45 INTEGER, INTENT(IN) :: VNO
46
47
48
49
50
51 DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
52
53
54 DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
55
56
57
58 DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
59
60
61
62 CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
63
64 DOUBLE PRECISION, INTENT(IN) :: TOL
65
66
67 CHARACTER(LEN=4), INTENT(IN) :: PC
68
69 INTEGER, INTENT(IN) :: ITMAX
70
71 INTEGER, INTENT(INOUT) :: IER
72
73
74
75
76
77
78
79
80 EXTERNAL LEQ_MATVEC, LEQ_MSOLVE, LEQ_MSOLVE0, LEQ_MSOLVE1
81
82
83
84 if(PC.eq.'LINE') then
85 call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m, &
86 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE, IER )
87 elseif(PC.eq.'DIAG') then
88 call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m, &
89 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE1, IER )
90 elseif(PC.eq.'NONE') then
91 call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m, &
92 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE0, IER )
93 else
94 IF(DMP_LOG)WRITE (UNIT_LOG,*) &
95 'preconditioner option not found - check mfix.dat and readme'
96 call mfix_exit(myPE)
97 endif
98
99 return
100 END SUBROUTINE LEQ_BICGS
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118 SUBROUTINE LEQ_BICGS0(VNAME, VNO, VAR, A_M, B_m, cmethod, &
119 TOL, ITMAX, MATVEC, MSOLVE, IER )
120
121
122
123
124 USE param
125 USE param1
126 USE parallel
127 USE matrix
128 USE geometry
129 USE compar
130 USE mpi_utility
131 USE sendrecv
132 USE indices
133 USE leqsol
134 USE cutcell
135 USE functions
136
137 IMPLICIT NONE
138
139
140
141
142 CHARACTER(LEN=*), INTENT(IN) :: Vname
143
144 INTEGER, INTENT(IN) :: VNO
145
146
147
148
149
150 DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
151
152
153 DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
154
155
156 DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
157
158
159 CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
160
161 DOUBLE PRECISION, INTENT(IN) :: TOL
162
163 INTEGER, INTENT(IN) :: ITMAX
164
165 INTEGER, INTENT(INOUT) :: IER
166
167
168
169
170
171
172 EXTERNAL MATVEC, MSOLVE
173
174
175
176 INTEGER, PARAMETER :: idebugl = 0
177 DOUBLE PRECISION, PARAMETER :: ratiotol = 0.2
178 LOGICAL, PARAMETER :: do_unit_scaling = .true.
179
180
181
182
183 DOUBLE PRECISION, DIMENSION(:), allocatable :: R,Rtilde, P,Phat, Svec, Shat, Tvec,V
184
185 DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: &
186 alpha, beta, omega, rho
187 DOUBLE PRECISION :: TxS, TxT, RtildexV, RtildexR, &
188 aijmax, oam
189 DOUBLE PRECISION :: Rnorm, Rnorm0, Snorm, TOLMIN, pnorm
190 LOGICAL :: isconverged
191 INTEGER :: i, j, k, ijk
192 INTEGER :: iter
193 DOUBLE PRECISION, DIMENSION(2) :: TxS_TxT
194
195
196
197
198 INTERFACE
199 DOUBLE PRECISION FUNCTION DOT_PRODUCT_PAR( R1, R2 )
200 use compar
201 USE param
202
203 DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_3) :: R1,R2
204 END FUNCTION DOT_PRODUCT_PAR
205 END INTERFACE
206
207 INTERFACE
208 FUNCTION DOT_PRODUCT_PAR2( R1, R2, R3, R4 )
209 use compar
210 USE param
211 USE functions
212 DOUBLE PRECISION, INTENT(IN), DIMENSION(ijkstart3:ijkend3) :: &
213 R1, R2, R3, R4
214 DOUBLE PRECISION, DIMENSION(2) :: DOT_PRODUCT_PAR2
215 END FUNCTION DOT_PRODUCT_PAR2
216 END INTERFACE
217
218
219
220 allocate(R(DIMENSION_3))
221 allocate(Rtilde(DIMENSION_3))
222 allocate(P(DIMENSION_3))
223 allocate(Phat(DIMENSION_3))
224 allocate(Svec(DIMENSION_3))
225 allocate(Shat(DIMENSION_3))
226 allocate(Tvec(DIMENSION_3))
227 allocate(V(DIMENSION_3))
228
229 is_serial = numPEs.eq.1.and.is_serial
230
231
232 = ZERO
233 rnorm0 = ZERO
234 snorm = ZERO
235 pnorm = ZERO
236
237
238 (:) = zero
239 beta(:) = zero
240 omega(:) = zero
241 rho(:) = zero
242
243 use_doloop = .TRUE.
244
245
246
247
248 if (use_doloop) then
249
250 do ijk=ijkstart3,ijkend3
251 R(ijk) = zero
252 Rtilde(ijk) = zero
253 P(ijk) = zero
254 Phat(ijk) = zero
255 Svec(ijk) = zero
256 Shat(ijk) = zero
257 Tvec(ijk) = zero
258 V(ijk) = zero
259 enddo
260
261 else
262 R(:) = zero
263 Rtilde(:) = zero
264 P(:) = zero
265 Phat(:) = zero
266 Svec(:) = zero
267 Shat(:) = zero
268 Tvec(:) = zero
269 V(:) = zero
270 endif
271
272 TOLMIN = EPSILON( one )
273
274
275
276 if (do_unit_scaling) then
277
278
279 IF(RE_INDEXING) THEN
280
281 DO IJK = IJKSTART3,IJKEND3
282 aijmax = maxval(abs(A_M(ijk,:)) )
283 OAM = one/aijmax
284 A_M(IJK,:) = A_M(IJK,:)*OAM
285 B_M(IJK) = B_M(IJK)*OAM
286 ENDDO
287
288 ELSE
289
290
291 do k = kstart2,kend2
292 do i = istart2,iend2
293 do j = jstart2,jend2
294 IJK = funijk(i,j,k)
295 aijmax = maxval(abs(A_M(ijk,:)) )
296 OAM = one/aijmax
297 A_M(IJK,:) = A_M(IJK,:)*OAM
298 B_M(IJK) = B_M(IJK)*OAM
299 enddo
300 enddo
301 enddo
302
303 ENDIF
304
305 endif
306
307
308
309
310
311
312
313 call MATVEC(Vname, Var, A_M, R)
314
315 if (use_doloop) then
316
317 do ijk=ijkstart3,ijkend3
318 R(ijk) = B_m(ijk) - R(ijk)
319 enddo
320 else
321 R(:) = B_m(:) - R(:)
322 endif
323
324 call send_recv(R,nlayers_bicgs)
325
326 if(is_serial) then
327 Rnorm0 = zero
328 if (use_doloop) then
329
330 do ijk=ijkstart3,ijkend3
331 Rnorm0 = Rnorm0 + R(ijk)*R(ijk)
332 enddo
333 else
334 Rnorm0 = dot_product(R,R)
335 endif
336 Rnorm0 = sqrt( Rnorm0 )
337 else
338 Rnorm0 = sqrt( dot_product_par( R, R ) )
339 endif
340
341
342
343
344
345
346 call random_number(Rtilde(:))
347
348
349 IF(RE_INDEXING) CALL SHIFT_DP_ARRAY(Rtilde)
350
351
352 if (use_doloop) then
353
354 do ijk=ijkstart3,ijkend3
355 Rtilde(ijk) = R(ijk) + (2.0d0*Rtilde(ijk)-1.0d0)*1.0d-6*Rnorm0
356 enddo
357 else
358 Rtilde(:) = R(:) + (2.0d0*Rtilde(:)-1.0d0)*1.0d-6*Rnorm0
359 endif
360
361 if (idebugl >= 1) then
362 if(myPE.eq.0) print*,'leq_bicgs, initial: ', Vname,' resid ', Rnorm0
363 endif
364
365
366
367
368
369 = 1
370 do i=1,itmax
371
372 if(is_serial) then
373 if (use_doloop) then
374 = zero
375
376 do ijk=ijkstart3,ijkend3
377 RtildexR = RtildexR + Rtilde(ijk) * R(ijk)
378 enddo
379 rho(i-1) = RtildexR
380 else
381 rho(i-1) = dot_product( Rtilde, R )
382 endif
383 else
384 rho(i-1) = dot_product_par( Rtilde, R )
385 endif
386
387
388
389 if (rho(i-1) .eq. zero) then
390 if(i /= 1)then
391
392
393
394 = -2
395 else
396
397
398 = 0
399 endif
400 call send_recv(var,2)
401 return
402 endif
403
404 if (i .eq. 1) then
405 if (use_doloop) then
406
407
408 do ijk=ijkstart3,ijkend3
409 P(ijk) = R(ijk)
410 enddo
411 else
412 P(:) = R(:)
413 endif
414 else
415 beta(i-1) = ( rho(i-1)/rho(i-2) )*( alpha(i-1) / omega(i-1) )
416 if (use_doloop) then
417
418 do ijk=ijkstart3,ijkend3
419 P(ijk) = R(ijk) + beta(i-1)*( P(ijk) - omega(i-1)*V(ijk) )
420 enddo
421 else
422 P(:) = R(:) + beta(i-1)*( P(:) - omega(i-1)*V(:) )
423 endif
424 endif
425
426
427
428
429
430 call MSOLVE(Vname, P, A_m, Phat, CMETHOD)
431 call MATVEC(Vname, Phat, A_m, V)
432
433 if(is_serial) then
434 if (use_doloop) then
435 RtildexV = zero
436
437 do ijk=ijkstart3,ijkend3
438 RtildexV = RtildexV + Rtilde(ijk) * V(ijk)
439 enddo
440 else
441 RtildexV = dot_product( Rtilde, V )
442 endif
443 else
444 RtildexV = dot_product_par( Rtilde, V )
445 endif
446
447
448
449
450
451 (i) = rho(i-1) / RtildexV
452
453
454
455 if (use_doloop) then
456
457
458 do ijk=ijkstart3,ijkend3
459 Svec(ijk) = R(ijk) - alpha(i) * V(ijk)
460 enddo
461
462 else
463 Svec(:) = R(:) - alpha(i) * V(:)
464 endif
465
466
467
468
469
470 if(.not.minimize_dotproducts) then
471 if(is_serial) then
472 if (use_doloop) then
473 Snorm = zero
474
475 do ijk=ijkstart3,ijkend3
476 Snorm = Snorm + Svec(ijk) * Svec(ijk)
477 enddo
478 else
479 Snorm = dot_product( Svec, Svec )
480 endif
481 Snorm = sqrt( Snorm )
482 else
483 Snorm = sqrt( dot_product_par( Svec, Svec ) )
484 endif
485
486
487
488 if (Snorm <= TOLMIN) then
489 if (use_doloop) then
490
491 do ijk=ijkstart3,ijkend3
492 Var(ijk) = Var(ijk) + alpha(i)*Phat(ijk)
493 enddo
494 else
495 Var(:) = Var(:) + alpha(i)*Phat(:)
496 endif
497
498
499
500 if (idebugl >= 1) then
501 call MATVEC(Vname, Var, A_m, R)
502
503
504
505 if (use_doloop) then
506
507 do ijk=ijkstart3,ijkend3
508 R(ijk) = B_m(ijk) - R(ijk)
509 enddo
510 else
511 R(:) = B_m(:) - R(:)
512 endif
513
514 if(is_serial) then
515 if (use_doloop) then
516 Rnorm = zero
517
518 do ijk=ijkstart3,ijkend3
519 Rnorm = Rnorm + R(ijk)*R(ijk)
520 enddo
521 else
522 Rnorm = dot_product( R, R )
523 endif
524 Rnorm = sqrt( Rnorm )
525 else
526 Rnorm = sqrt( dot_product_par( R, R ) )
527 endif
528
529 endif
530 = .TRUE.
531 EXIT
532 endif
533 endif
534
535
536
537
538
539
540 call MSOLVE( Vname, Svec, A_m, Shat, CMETHOD)
541 call MATVEC( Vname, Shat, A_m, Tvec )
542
543 if(is_serial) then
544 if (use_doloop) then
545 TxS = zero
546 TxT = zero
547
548 do ijk=ijkstart3,ijkend3
549 TxS = TxS + Tvec(ijk) * Svec(ijk)
550 TxT = TxT + Tvec(ijk) * Tvec(ijk)
551 enddo
552 else
553 TxS = dot_product( Tvec, Svec )
554 TxT = dot_product( Tvec, Tvec )
555 endif
556 else
557 if(.not.minimize_dotproducts) then
558 TxS = dot_product_par( Tvec, Svec )
559 TxT = dot_product_par( Tvec, Tvec )
560 else
561 TxS_TxT = dot_product_par2(Tvec, Svec, Tvec, Tvec )
562 TxS = TxS_TxT(1)
563 TxT = TxS_TxT(2)
564 endif
565 endif
566
567 IF(TxT.eq.Zero) TxT = SMALL_NUMBER
568
569
570
571 (i) = TxS / TxT
572
573
574
575 if (use_doloop) then
576
577 do ijk=ijkstart3,ijkend3
578 Var(ijk) = Var(ijk) + &
579 alpha(i)*Phat(ijk) + omega(i)*Shat(ijk)
580 R(ijk) = Svec(ijk) - omega(i)*Tvec(ijk)
581 enddo
582 else
583 Var(:) = Var(:) + &
584 alpha(i)*Phat(:) + omega(i)*Shat(:)
585 R(:) = Svec(:) - omega(i)*Tvec(:)
586 endif
587
588
589
590 if(.not.minimize_dotproducts.or.(mod(iter,icheck_bicgs).eq.0)) then
591 if(is_serial) then
592 if (use_doloop) then
593 Rnorm = zero
594
595 do ijk=ijkstart3,ijkend3
596 Rnorm = Rnorm + R(ijk) * R(ijk)
597 enddo
598 else
599 Rnorm = dot_product(R, R )
600 endif
601 Rnorm = sqrt( Rnorm )
602 else
603 Rnorm = sqrt( dot_product_par(R, R) )
604 endif
605
606 if (idebugl.ge.1) then
607 if (myPE.eq.PE_IO) then
608 print*,'iter, Rnorm ', iter, Rnorm, Snorm
609 print*,'alpha(i), omega(i) ', alpha(i), omega(i)
610 print*,'TxS, TxT ', TxS, TxT
611 print*,'RtildexV, rho(i-1) ', RtildexV, rho(i-1)
612 endif
613 endif
614
615
616
617
618
619 = (Rnorm <= TOL*Rnorm0)
620
621 if (isconverged) then
622 iter_tot(vno) = iter_tot(vno) + iter + 1
623 EXIT
624 endif
625 endif
626
627
628 = iter + 1
629
630 enddo
631
632
633
634
635 if (idebugl >= 1) then
636 call MATVEC(Vname, Var, A_m, R)
637 if (use_doloop) then
638
639 do ijk=ijkstart3,ijkend3
640 R(ijk) = R(ijk) - B_m(ijk)
641 enddo
642 else
643 R(:) = R(:) - B_m(:)
644 endif
645
646 if(is_serial) then
647 if (use_doloop) then
648 Rnorm = zero
649
650 do ijk=ijkstart3,ijkend3
651 Rnorm = Rnorm + R(ijk) * R(ijk)
652 enddo
653 else
654 Rnorm = dot_product( R,R)
655 endif
656 Rnorm = sqrt( Rnorm )
657 else
658 Rnorm = sqrt( dot_product_par( R,R) )
659 endif
660
661 if(myPE.eq.0) print*,'leq_bicgs: final Rnorm ', Rnorm
662
663 if(myPE.eq.0) print*,'leq_bicgs ratio : ', Vname,' ',iter, &
664 ' L-2', Rnorm/Rnorm0
665 endif
666
667
668 if(.NOT.isConverged) isconverged = (real(Rnorm) <= TOL*Rnorm0);
669
670 = 0
671 if (.not.isconverged) then
672 IER = -1
673 iter_tot(vno) = iter_tot(vno) + iter
674 if (real(Rnorm) >= ratiotol*real(Rnorm0)) then
675 IER = -2
676 endif
677 endif
678
679 call send_recv(var,2)
680
681 deallocate(R)
682 deallocate(Rtilde)
683 deallocate(P)
684 deallocate(Phat)
685 deallocate(Svec)
686 deallocate(Shat)
687 deallocate(Tvec)
688 deallocate(V)
689
690 return
691 end subroutine LEQ_BICGS0
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708 SUBROUTINE LEQ_ISWEEP(I, Vname, VAR, A_M, B_M)
709
710
711
712
713 USE param
714 USE param1
715 USE matrix
716 USE geometry
717 USE compar
718 USE indices
719 USE funits
720 USE sendrecv
721 USE mpi_utility
722 USE functions
723 IMPLICIT NONE
724
725
726
727
728 INTEGER, INTENT(IN) :: I
729
730 CHARACTER(LEN=*), INTENT(IN) :: Vname
731
732 DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
733
734 DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
735
736 DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
737
738
739
740 DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC, DD, EE, BB
741 INTEGER :: NSTART, NEND, INFO
742 INTEGER :: IJK, J, K, IM1JK, IP1JK
743
744
745 = JEND
746 NSTART = JSTART
747 K = 1
748
749 DO J=NSTART, NEND
750 IJK = FUNIJK(I,J,K)
751 IM1JK = IM_OF(IJK)
752 IP1JK = IP_OF(IJK)
753 DD(J) = A_M(IJK, 0)
754 CC(J) = A_M(IJK, -2)
755 EE(J) = A_M(IJK, 2)
756 BB(J) = B_M(IJK) - A_M(IJK,-1) * Var( IM1JK ) &
757 - A_M(IJK, 1) * Var( IP1JK )
758 ENDDO
759
760 CC(NSTART) = ZERO
761 EE(NEND) = ZERO
762 INFO = 0
763
764 CALL DGTSV( JEND-JSTART+1, 1, CC(JSTART+1), DD, EE, BB, JEND-JSTART+1, INFO )
765
766 IF (INFO.NE.0) THEN
767 RETURN
768 ENDIF
769
770 DO J=NSTART, NEND
771 IJK = FUNIJK(I,J,K)
772 Var(IJK) = BB(J)
773 ENDDO
774
775 RETURN
776 END SUBROUTINE LEQ_ISWEEP
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793 SUBROUTINE LEQ_IKSWEEP(I, K, Vname, VAR, A_M, B_M)
794
795
796
797
798 USE param
799 USE param1
800 USE matrix
801 USE geometry
802 USE compar
803 USE funits
804 USE indices
805 USE sendrecv
806 USE mpi_utility
807 USE functions
808 IMPLICIT NONE
809
810
811
812
813 INTEGER, INTENT(IN) :: I, K
814
815 CHARACTER(LEN=*), INTENT(IN) :: Vname
816
817 DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
818
819 DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
820
821 DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
822
823
824
825 DOUBLE PRECISION, DIMENSION(JSTART:JEND) :: CC, DD, EE, BB
826 INTEGER :: NSTART, NEND, INFO
827 INTEGER :: IJK, J, IM1JK, IP1JK, IJKM1, IJKP1
828
829
830 = JEND
831 NSTART = JSTART
832
833
834 DO J=NSTART, NEND
835
836 = FUNIJK(I,J,K)
837 IM1JK = IM_OF(IJK)
838 IP1JK = IP_OF(IJK)
839 IJKM1 = KM_OF(IJK)
840 IJKP1 = KP_OF(IJK)
841 DD(J) = A_M(IJK, 0)
842 CC(J) = A_M(IJK, -2)
843 EE(J) = A_M(IJK, 2)
844 BB(J) = B_M(IJK) - A_M(IJK,-1) * Var( IM1JK ) &
845 - A_M(IJK, 1) * Var( IP1JK ) &
846 - A_M(IJK,-3) * Var( IJKM1 ) &
847 - A_M(IJK, 3) * Var( IJKP1 )
848 ENDDO
849
850 CC(NSTART) = ZERO
851 EE(NEND) = ZERO
852 INFO = 0
853
854 CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
855
856 IF (INFO.NE.0) THEN
857 write(*,*) 'leq_iksweep',INFO, myPE
858 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' IKSWEEP'
859 RETURN
860 ENDIF
861
862 DO J=NSTART, NEND
863 IJK = FUNIJK(I,J,K)
864 Var(IJK) = BB(J)
865 ENDDO
866
867 RETURN
868 END SUBROUTINE LEQ_IKSWEEP
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886 SUBROUTINE LEQ_JKSWEEP(J, K, Vname, VAR, A_M, B_M)
887
888
889
890
891 USE param
892 USE param1
893 USE matrix
894 USE geometry
895 USE funits
896 USE compar
897 USE indices
898 USE functions
899 IMPLICIT NONE
900
901
902
903
904 INTEGER, INTENT(IN) :: J, K
905
906 CHARACTER(LEN=*), INTENT(IN) :: Vname
907
908 DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
909
910 DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
911
912 DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
913
914
915
916 DOUBLE PRECISION, DIMENSION (ISTART:IEND) :: CC, DD, EE, BB
917 INTEGER :: NSTART, NEND, INFO, IJK, I
918
919
920 = IEND
921 NSTART = ISTART
922
923 DO I=NSTART,NEND
924 IJK = FUNIJK(I,J,K)
925 DD(I) = A_M(IJK, 0)
926 CC(I) = A_M(IJK, -1)
927 EE(I) = A_M(IJK, 1)
928 BB(I) = B_M(IJK) - A_M(IJK,-2) * Var( JM_OF(IJK) ) &
929 - A_M(IJK, 2) * Var( JP_OF(IJK) ) &
930 - A_M(IJK,-3) * Var( KM_OF(IJK) ) &
931 - A_M(IJK, 3) * Var( KP_OF(IJK) )
932 ENDDO
933
934 CC(NSTART) = ZERO
935 EE(NEND) = ZERO
936 INFO = 0
937 CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
938
939 IF (INFO.NE.0) THEN
940 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
941 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' JKSWEEP'
942 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
943 call mfix_exit(myPE)
944 ENDIF
945
946 DO I=NSTART, NEND
947 IJK = FUNIJK(I,J,K)
948 Var(IJK) = BB(I)
949 ENDDO
950
951 RETURN
952 END SUBROUTINE LEQ_JKSWEEP
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970 SUBROUTINE LEQ_IJSWEEP(I, J, Vname, VAR, A_M, B_M)
971
972
973
974
975 USE param
976 USE param1
977 USE matrix
978 USE geometry
979 USE funits
980 USE compar
981 USE indices
982 USE functions
983 IMPLICIT NONE
984
985
986
987
988 INTEGER, INTENT(IN) :: I, J
989
990 CHARACTER(LEN=*), INTENT(IN) :: Vname
991
992 DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
993
994 DOUBLE PRECISION, INTENT(IN) :: A_m(ijkstart3:ijkend3, -3:3)
995
996 DOUBLE PRECISION, INTENT(IN) :: B_m(ijkstart3:ijkend3)
997
998
999
1000 DOUBLE PRECISION, DIMENSION (KSTART:KEND) :: CC, DD, EE, BB
1001 INTEGER :: NEND, NSTART, INFO, IJK, K
1002
1003
1004 = KEND
1005 NSTART = KSTART
1006
1007 DO K=NSTART, NEND
1008 IJK = FUNIJK(I,J,K)
1009 DD(K) = A_M(IJK, 0)
1010 CC(K) = A_M(IJK, -3)
1011 EE(K) = A_M(IJK, 3)
1012 BB(K) = B_M(IJK) - A_M(IJK,-2) * Var( JM_OF(IJK) ) &
1013 - A_M(IJK, 2) * Var( JP_OF(IJK) ) &
1014 - A_M(IJK,-1) * Var( IM_OF(IJK) ) &
1015 - A_M(IJK, 1) * Var( IP_OF(IJK) )
1016 ENDDO
1017
1018 CC(NSTART) = ZERO
1019 EE(NEND) = ZERO
1020 INFO = 0
1021 CALL DGTSV(NEND-NSTART+1, 1, CC(NSTART+1), DD, EE, BB, NEND-NSTART+1, INFO)
1022
1023 IF (INFO.NE.0) THEN
1024 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
1025 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'ROUTINE = ', ' IJSWEEP'
1026 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1027 call mfix_exit(myPE)
1028 ENDIF
1029
1030 DO K=NSTART, NEND
1031 IJK = FUNIJK(I,J,K)
1032 Var(IJK) = BB(K)
1033 ENDDO
1034
1035 RETURN
1036 END SUBROUTINE LEQ_IJSWEEP
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055 SUBROUTINE LEQ_MATVEC(VNAME, VAR, A_M, Avar)
1056
1057
1058
1059
1060 USE param
1061 USE param1
1062 USE matrix
1063 USE geometry
1064 USE compar
1065 USE indices
1066 USE sendrecv
1067 USE mpi_utility
1068 USE cutcell
1069 USE functions
1070 IMPLICIT NONE
1071
1072
1073
1074
1075 CHARACTER(LEN=*), INTENT(IN) :: Vname
1076
1077
1078 DOUBLE PRECISION, INTENT(IN) :: Var(DIMENSION_3)
1079
1080
1081 DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1082
1083
1084 DOUBLE PRECISION, INTENT(OUT) :: AVar(DIMENSION_3)
1085
1086
1087
1088
1089 INTEGER :: I, J, K, IJK
1090 integer :: im1jk, ip1jk, ijm1k, ijp1k, ijkm1, ijkp1
1091
1092
1093 IF(RE_INDEXING) THEN
1094
1095 DO IJK = IJKSTART3, IJKEND3
1096
1097 = im_of(ijk)
1098 ip1jk = ip_of(ijk)
1099 ijm1k = jm_of(ijk)
1100 ijp1k = jp_of(ijk)
1101
1102
1103 AVar(ijk) = A_m(ijk,-2) * Var(ijm1k) &
1104 + A_m(ijk,-1) * Var(im1jk) &
1105 + A_m(ijk, 0) * Var(ijk) &
1106 + A_m(ijk, 1) * Var(ip1jk) &
1107 + A_m(ijk, 2) * Var(ijp1k)
1108
1109 if (do_k) then
1110 ijkm1 = km_of(ijk)
1111 ijkp1 = kp_of(ijk)
1112
1113
1114 AVar(ijk) = AVar(ijk) + A_m(ijk,-3) * Var(ijkm1) &
1115 + A_m(ijk, 3) * Var(ijkp1)
1116
1117 endif
1118
1119 enddo
1120
1121
1122 ELSE
1123
1124
1125
1126 if (do_k) then
1127
1128
1129
1130
1131 do k = kstart,kend
1132 do i = istart,iend
1133 do j = jstart,jend
1134 IJK = funijk(i,j,k)
1135 im1jk = im_of(ijk)
1136 ip1jk = ip_of(ijk)
1137 ijm1k = jm_of(ijk)
1138 ijp1k = jp_of(ijk)
1139 ijkm1 = km_of(ijk)
1140 ijkp1 = kp_of(ijk)
1141 AVar(ijk) = A_m(ijk,-3) * Var(ijkm1) &
1142 + A_m(ijk,-2) * Var(ijm1k) &
1143 + A_m(ijk,-1) * Var(im1jk) &
1144 + A_m(ijk, 0) * Var(ijk) &
1145 + A_m(ijk, 1) * Var(ip1jk) &
1146 + A_m(ijk, 2) * Var(ijp1k) &
1147 + A_m(ijk, 3) * Var(ijkp1)
1148 enddo
1149 enddo
1150 enddo
1151
1152 else
1153 k = 1
1154
1155 do i = istart,iend
1156 do j = jstart,jend
1157 IJK = funijk(i,j,k)
1158 im1jk = im_of(ijk)
1159 ip1jk = ip_of(ijk)
1160 ijm1k = jm_of(ijk)
1161 ijp1k = jp_of(ijk)
1162 AVar(ijk) = A_m(ijk,-2) * Var(ijm1k) &
1163 + A_m(ijk,-1) * Var(im1jk) &
1164 + A_m(ijk, 0) * Var(ijk) &
1165 + A_m(ijk, 1) * Var(ip1jk) &
1166 + A_m(ijk, 2) * Var(ijp1k)
1167 enddo
1168 enddo
1169
1170 endif
1171
1172 ENDIF
1173
1174 call send_recv(Avar,nlayers_bicgs)
1175 RETURN
1176 END SUBROUTINE LEQ_MATVEC
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197 SUBROUTINE LEQ_MSOLVE(VNAME, B_m, A_M, Var, CMETHOD)
1198
1199
1200
1201
1202 USE param
1203 USE param1
1204 USE matrix
1205 USE geometry
1206 USE compar
1207 USE indices
1208 USE sendrecv
1209 USE functions
1210 IMPLICIT NONE
1211
1212
1213
1214
1215 CHARACTER(LEN=*), INTENT(IN) :: Vname
1216
1217
1218 DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
1219
1220
1221 DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1222
1223
1224 DOUBLE PRECISION, INTENT(INOUT) :: Var(DIMENSION_3)
1225
1226
1227 CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
1228
1229
1230
1231 LOGICAL, PARAMETER :: USE_IKLOOP = .FALSE.
1232 LOGICAL, PARAMETER :: SETGUESS = .TRUE.
1233
1234
1235
1236
1237 INTEGER :: ITER, NITER
1238 INTEGER :: IJK, I , J, K
1239 INTEGER :: I1, J1, K1, I2, J2, K2, IK, JK, IJ
1240 INTEGER :: ISIZE, JSIZE, KSIZE
1241 INTEGER :: ICASE
1242
1243
1244 CHARACTER :: CH
1245 LOGICAL :: DO_ISWEEP, DO_JSWEEP, DO_KSWEEP
1246 LOGICAL :: DO_SENDRECV, DO_REDBLACK, DO_ALL
1247
1248
1249
1250
1251
1252
1253
1254 IF (SETGUESS) THEN
1255
1256 do k = kstart3,kend3
1257 do i = istart3,iend3
1258 do j = jstart3,jend3
1259 IJK = funijk(i,j,k)
1260 VAR(IJK) = B_M(IJK)
1261 enddo
1262 enddo
1263 enddo
1264
1265 call send_recv(var,nlayers_bicgs)
1266 ENDIF
1267
1268 NITER = LEN( CMETHOD )
1269
1270 DO ITER=1,NITER
1271
1272
1273 = CMETHOD( ITER:ITER )
1274 DO_ISWEEP = (CH .EQ. 'I') .OR. (CH .EQ. 'i')
1275 DO_JSWEEP = (CH .EQ. 'J') .OR. (CH .EQ. 'j')
1276 DO_KSWEEP = (CH .EQ. 'K') .OR. (CH .EQ. 'k')
1277 DO_ALL = (CH .EQ. 'A') .OR. (CH .EQ. 'a')
1278 DO_REDBLACK = (CH .EQ. 'R') .OR. (CH .EQ. 'r')
1279 DO_SENDRECV = (CH .EQ. 'S') .OR. (CH .EQ. 's')
1280
1281 IF (NO_K) THEN
1282
1283 IF ( DO_ISWEEP ) THEN
1284
1285 DO I=istart,iend,1
1286 CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
1287 ENDDO
1288 ENDIF
1289
1290
1291 IF (DO_REDBLACK) THEN
1292
1293 DO I=istart,iend,2
1294 CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
1295 ENDDO
1296
1297 DO I=istart+1,iend,2
1298 CALL LEQ_ISWEEP( I, Vname, Var, A_m, B_m )
1299 ENDDO
1300 ENDIF
1301
1302 ELSE
1303
1304
1305
1306
1307 IF(DO_ALL) THEN
1308
1309
1310 = jstart
1311 k1 = kstart
1312 j2 = jend
1313 k2 = kend
1314 jsize = j2-j1+1
1315 ksize = k2-k1+1
1316 DO icase = 1, 2
1317
1318 DO JK=icase, ksize*jsize, 2
1319 if (mod(jk,jsize).ne.0) then
1320 k = int( jk/jsize ) + k1
1321 else
1322 k = int( jk/jsize ) + k1 -1
1323 endif
1324 j = (jk-1-(k-k1)*jsize) + j1 + mod(k,2)
1325 if(j.gt.j2) j=j-j2 + j1 -1
1326 CALL LEQ_JKSWEEP(J, K, Vname, Var, A_m, B_m)
1327 ENDDO
1328
1329 ENDDO
1330 call send_recv(var,nlayers_bicgs)
1331
1332
1333
1334 = istart
1335 j1 = jstart
1336 i2 = iend
1337 j2 = jend
1338 isize = i2-i1+1
1339 jsize = j2-j1+1
1340 DO icase = 1, 2
1341
1342 DO IJ=icase, jsize*isize, 2
1343 if (mod(ij,isize).ne.0) then
1344 j = int( ij/isize ) + j1
1345 else
1346 j = int( ij/isize ) + j1 -1
1347 endif
1348 i = (ij-1-(j-j1)*isize) + i1 + mod(j,2)
1349 if(i.gt.i2) i=i-i2 + i1 -1
1350 CALL LEQ_IJSWEEP(I, J, Vname, Var, A_m, B_m)
1351 ENDDO
1352
1353 ENDDO
1354 call send_recv(var,nlayers_bicgs)
1355
1356
1357
1358 = istart
1359 k1 = kstart
1360 i2 = iend
1361 k2 = kend
1362 isize = i2-i1+1
1363 ksize = k2-k1+1
1364
1365 DO icase = 1, 2
1366
1367 DO IK=icase, ksize*isize, 2
1368 if (mod(ik,isize).ne.0) then
1369 k = int( ik/isize ) + k1
1370 else
1371 k = int( ik/isize ) + k1 -1
1372 endif
1373 i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
1374 if(i.gt.i2) i=i-i2 + i1 -1
1375 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1376 ENDDO
1377
1378 ENDDO
1379 ENDIF
1380
1381
1382
1383
1384 IF(DO_REDBLACK) THEN
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407 DO k=kstart,kend
1408 IF(mod(k,2).ne.0)THEN
1409 DO I=istart+1,iend,2
1410 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1411 ENDDO
1412 ELSE
1413 DO I=istart,iend,2
1414 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1415 ENDDO
1416 ENDIF
1417 ENDDO
1418
1419
1420 DO k=kstart,kend
1421 IF(mod(k,2).ne.0)THEN
1422 DO I=istart,iend,2
1423 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1424 ENDDO
1425 ELSE
1426 DO I=istart+1,iend,2
1427 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1428 ENDDO
1429 ENDIF
1430 ENDDO
1431
1432
1433 ENDIF
1434
1435
1436
1437
1438 IF(USE_IKLOOP) THEN
1439
1440
1441 = istart
1442 k1 = kstart
1443 i2 = iend
1444 k2 = kend
1445 isize = i2-i1+1
1446 ksize = k2-k1+1
1447 IF (DO_ISWEEP) THEN
1448
1449 DO IK=1, ksize*isize
1450 if (mod(ik,isize).ne.0) then
1451 k = int( ik/isize ) + k1
1452 else
1453 k = int( ik/isize ) + k1 -1
1454 endif
1455 i = (ik-1-(k-k1)*isize) + i1
1456 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1457 ENDDO
1458 ENDIF
1459 IF (DO_KSWEEP) THEN
1460
1461 DO IK=1, ksize*isize
1462 if (mod(ik,ksize).ne.0) then
1463 i = int( ik/ksize ) + i1
1464 else
1465 i = int( ik/ksize ) + i1 -1
1466 endif
1467 k = (ik-1-(i-i1)*ksize) + k1
1468 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1469 ENDDO
1470 ENDIF
1471
1472 ELSE
1473
1474
1475
1476 IF (DO_ISWEEP) THEN
1477
1478 DO K=kstart,kend
1479 DO I=istart,iend
1480 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1481 ENDDO
1482 ENDDO
1483 ENDIF
1484 IF (DO_KSWEEP) THEN
1485
1486 DO I=istart,iend
1487 DO K=kstart,kend
1488 CALL LEQ_IKSWEEP(I, K, Vname, Var, A_m, B_m)
1489 ENDDO
1490 ENDDO
1491 ENDIF
1492 ENDIF
1493
1494
1495 ENDIF
1496
1497
1498
1499 IF (DO_SENDRECV) call send_recv(var,nlayers_bicgs)
1500
1501
1502 ENDDO
1503
1504
1505
1506 RETURN
1507 END SUBROUTINE LEQ_MSOLVE
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525 SUBROUTINE LEQ_MSOLVE0(VNAME, B_m, A_M, Var, CMETHOD)
1526
1527
1528
1529
1530 USE param
1531 USE param1
1532 USE matrix
1533 USE geometry
1534 USE compar
1535 USE indices
1536 USE sendrecv
1537 use parallel
1538 IMPLICIT NONE
1539
1540
1541
1542
1543 CHARACTER(LEN=*), INTENT(IN) :: Vname
1544
1545
1546 DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
1547
1548
1549 DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1550
1551
1552 DOUBLE PRECISION, INTENT(OUT) :: Var(DIMENSION_3)
1553
1554 CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
1555
1556
1557
1558 integer :: ijk
1559
1560
1561
1562 if (use_doloop) then
1563
1564 do ijk=ijkstart3,ijkend3
1565 var(ijk) = b_m(ijk)
1566 enddo
1567 else
1568 var(:) = b_m(:)
1569 endif
1570 call send_recv(var,nlayers_bicgs)
1571
1572 return
1573 end subroutine leq_msolve0
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591 SUBROUTINE LEQ_MSOLVE1(VNAME, B_m, A_M, Var, CMETHOD)
1592
1593
1594
1595
1596 USE param
1597 USE param1
1598 USE matrix
1599 USE geometry
1600 USE compar
1601 USE indices
1602 USE sendrecv
1603 use parallel
1604
1605 USE cutcell
1606 USE functions
1607 IMPLICIT NONE
1608
1609
1610
1611
1612 CHARACTER(LEN=*), INTENT(IN) :: Vname
1613
1614
1615 DOUBLE PRECISION, INTENT(IN) :: B_m(DIMENSION_3)
1616
1617
1618 DOUBLE PRECISION, INTENT(IN) :: A_m(DIMENSION_3, -3:3)
1619
1620
1621 DOUBLE PRECISION, INTENT(OUT) :: Var(DIMENSION_3)
1622
1623 CHARACTER(LEN=4), INTENT(IN) :: CMETHOD
1624
1625
1626
1627 integer :: i,j,k, ijk
1628
1629
1630 if (use_doloop) then
1631
1632 do ijk=ijkstart3,ijkend3
1633 var(ijk) = zero
1634 enddo
1635 else
1636 var(:) = ZERO
1637 endif
1638
1639
1640 IF(.NOT.RE_INDEXING) THEN
1641
1642 do k=kstart2,kend2
1643 do i=istart2,iend2
1644 do j=jstart2,jend2
1645 ijk = funijk( i,j,k )
1646 var(ijk) = b_m(ijk)/A_m(ijk,0)
1647 enddo
1648 enddo
1649 enddo
1650 ELSE
1651
1652 DO IJK=IJKSTART3,IJKEND3
1653 var(ijk) = b_m(ijk)/A_m(ijk,0)
1654 ENDDO
1655 ENDIF
1656
1657 call send_recv(var,nlayers_bicgs)
1658
1659 return
1660 end subroutine leq_msolve1
1661
1662
1663
1664
1665
1666
1667
1668
1669 double precision function dot_product_par(r1,r2)
1670
1671
1672
1673
1674 use mpi_utility
1675 use geometry
1676 use compar
1677 use indices
1678 use cutcell
1679 use functions
1680 implicit none
1681
1682
1683
1684
1685 double precision, intent(in), dimension(DIMENSION_3) :: r1,r2
1686
1687
1688
1689 logical, parameter :: do_global_sum = .true.
1690
1691
1692
1693 DOUBLE PRECISION, allocatable, Dimension(:) :: r1_g, r2_g
1694 double precision :: prod
1695 integer :: i, j, k, ijk
1696
1697
1698 if(do_global_sum) then
1699 prod = 0.0d0
1700
1701 IF(RE_INDEXING) THEN
1702
1703
1704 DO IJK = IJKSTART3,IJKEND3
1705 IF(INTERIOR_CELL_AT(IJK)) prod = prod + r1(ijk)*r2(ijk)
1706 ENDDO
1707
1708 call global_all_sum(prod, dot_product_par)
1709
1710
1711 ELSE
1712
1713
1714
1715
1716 do k = kstart1, kend1
1717 do i = istart1, iend1
1718 do j = jstart1, jend1
1719 ijk = funijk_map_c (i,j,k)
1720 prod = prod + r1(ijk)*r2(ijk)
1721 enddo
1722 enddo
1723 enddo
1724
1725 call global_all_sum(prod, dot_product_par)
1726
1727 ENDIF
1728
1729 else
1730 if(myPE.eq.root) then
1731 allocate (r1_g(1:ijkmax3))
1732 allocate (r2_g(1:ijkmax3))
1733 else
1734 allocate (r1_g(10))
1735 allocate (r2_g(10))
1736 endif
1737 call gather(r1,r1_g)
1738 call gather(r2,r2_g)
1739
1740 if(myPE.eq.root) then
1741 prod = 0.0d0
1742
1743
1744 do k = kmin1, kmax1
1745 do i = imin1, imax1
1746 do j = jmin1, jmax1
1747 ijk = funijk_gl (imap_c(i),jmap_c(j),kmap_c(k))
1748
1749 = prod + r1_g(ijk)*r2_g(ijk)
1750 enddo
1751 enddo
1752 enddo
1753
1754 endif
1755 call bcast( prod)
1756
1757 dot_product_par = prod
1758
1759 deallocate (r1_g)
1760 deallocate (r2_g)
1761
1762 endif
1763
1764 end function dot_product_par
1765
1766
1767
1768
1769
1770
1771
1772
1773 function dot_product_par2(r1,r2,r3,r4)
1774
1775
1776
1777
1778 use mpi_utility
1779 use geometry
1780 use compar
1781 use indices
1782 use functions
1783 implicit none
1784
1785
1786
1787 double precision, intent(in), dimension(ijkstart3:ijkend3) :: r1,r2,r3,r4
1788
1789
1790
1791 logical, parameter :: do_global_sum = .true.
1792
1793
1794
1795 DOUBLE PRECISION, allocatable, Dimension(:,:) :: r_temp, rg_temp
1796 double precision, Dimension(2) :: prod, dot_product_par2
1797 integer :: i, j, k, ijk
1798
1799
1800 if(do_global_sum) then
1801
1802 prod(:) = 0.0d0
1803
1804
1805 do k = kstart1, kend1
1806 do i = istart1, iend1
1807 do j = jstart1, jend1
1808
1809 IF (DEAD_CELL_AT(I,J,K)) CYCLE
1810
1811 = funijk_map_c (i,j,k)
1812 prod(1) = prod(1) + r1(ijk)*r2(ijk)
1813 prod(2) = prod(2) + r3(ijk)*r4(ijk)
1814 enddo
1815 enddo
1816 enddo
1817
1818 call global_all_sum(prod, dot_product_par2)
1819
1820 else
1821 allocate (r_temp(ijkstart3:ijkend3,4))
1822 r_temp(:,1) = r1
1823 r_temp(:,2) = r2
1824 r_temp(:,3) = r3
1825 r_temp(:,4) = r4
1826
1827 if(myPE.eq.root) then
1828 allocate (rg_temp(1:ijkmax3,4))
1829 else
1830 allocate (rg_temp(10,4))
1831 endif
1832 call gather(r_temp,rg_temp)
1833
1834 if(myPE.eq.root) then
1835 prod = 0.0d0
1836
1837 do k = kmin1, kmax1
1838 do i = imin1, imax1
1839 do j = jmin1, jmax1
1840 IF (DEAD_CELL_AT(I,J,K)) CYCLE
1841 = funijk_gl (imap_c(i),jmap_c(j),kmap_c(k))
1842
1843 (1) = prod(1) + rg_temp(ijk,1)*rg_temp(ijk,2)
1844 prod(2) = prod(2) + rg_temp(ijk,3)*rg_temp(ijk,4)
1845 enddo
1846 enddo
1847 enddo
1848 endif
1849 call bcast( prod)
1850
1851 dot_product_par2 = prod
1852
1853 deallocate (r_temp)
1854 deallocate (rg_temp)
1855
1856 endif
1857
1858 end function dot_product_par2
1859
1860
1861
1862