File: /nfs/home/0/users/jenkins/mfix.git/model/leq_bicgst.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 SUBROUTINE LEQ_BICGSt(VNAME, VNO, VAR, A_M, B_m, cmethod, TOL, PC, ITMAX,IER)
21
22
23
24
25 USE param
26 USE param1
27 USE matrix
28 USE geometry
29 USE compar
30 USE indices
31 USE leqsol
32 USE funits
33 IMPLICIT NONE
34
35
36
37
38
39 INTEGER :: IER
40
41 INTEGER :: ITMAX
42
43 INTEGER :: VNO
44
45 DOUBLE PRECISION :: TOL
46
47 CHARACTER(LEN=4) :: PC
48
49 DOUBLE PRECISION, DIMENSION(-3:3,ijkstart3:ijkend3) :: A_m
50
51 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: B_m
52
53 CHARACTER(LEN=*) :: Vname
54
55 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: Var
56
57 CHARACTER(LEN=*) :: CMETHOD
58
59
60 EXTERNAL LEQ_MATVECt, LEQ_MSOLVEt, LEQ_MSOLVE0t, LEQ_MSOLVE1t
61
62
63 if(PC.eq.'LINE') then
64 call LEQ_BICGS0t( Vname, Vno, Var, A_m, B_m, &
65 cmethod, TOL, ITMAX, LEQ_MATVECt, LEQ_MSOLVEt, IER )
66 elseif(PC.eq.'DIAG') then
67 call LEQ_BICGS0t( Vname, Vno, Var, A_m, B_m, &
68 cmethod, TOL, ITMAX, LEQ_MATVECt, LEQ_MSOLVE1t, IER )
69 elseif(PC.eq.'NONE') then
70 call LEQ_BICGS0t( Vname, Vno, Var, A_m, B_m, &
71 cmethod, TOL, ITMAX, LEQ_MATVECt, LEQ_MSOLVE0t, IER )
72 else
73 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'preconditioner option not found - check mfix.dat and readme'
74 call mfix_exit(myPE)
75 endif
76
77 return
78 END SUBROUTINE LEQ_BICGSt
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102 SUBROUTINE LEQ_BICGS0t(VNAME, VNO, VAR, A_M, B_m, cmethod, TOL, ITMAX, &
103 MATVECt, MSOLVEt, IER )
104
105
106
107 USE param
108 USE param1
109 USE parallel
110 USE matrix
111 USE geometry
112 USE compar
113 USE mpi_utility
114 USE sendrecv
115 USE indices
116 USE leqsol
117 USE functions
118 IMPLICIT NONE
119
120
121
122
123
124
125 INTEGER :: IER
126
127 INTEGER :: ITMAX
128
129 INTEGER :: VNO
130
131 DOUBLE PRECISION :: TOL
132
133 DOUBLE PRECISION, DIMENSION(-3:3,ijkstart3:ijkend3) :: A_m
134
135 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: B_m
136
137 CHARACTER(LEN=*) :: Vname
138
139 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: Var
140
141 CHARACTER(LEN=*) :: CMETHOD
142
143
144
145
146 INTEGER, PARAMETER :: idebugl = 1
147 DOUBLE PRECISION :: ratiotol = 0.2
148
149 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: &
150 R,Rtilde, P,Phat, Svec, Shat, Tvec,V
151 DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: alpha,beta,omega,rho
152 DOUBLE PRECISION :: TxS, TxT, oam,RtildexV, &
153 RtildexR, aijmax, Rnorm=0, Rnorm0, Snorm, TOLMIN, pnorm
154 LOGICAL :: isconverged
155 INTEGER :: i, j, k, iter
156
157 INTEGER :: ijk2
158 DOUBLE PRECISION, DIMENSION(2) :: TxS_TxT
159
160
161
162
163
164 EXTERNAL MATVECt, MSOLVEt
165
166 INTERFACE
167 DOUBLE PRECISION FUNCTION DOT_PRODUCT_PAR( R1, R2 )
168 use compar
169 use param
170 DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_3) :: R1,R2
171 END FUNCTION DOT_PRODUCT_PAR
172 END INTERFACE
173
174 INTERFACE
175 FUNCTION DOT_PRODUCT_PAR2( R1, R2, R3, R4 )
176 use compar
177 DOUBLE PRECISION, INTENT(IN), DIMENSION(ijkstart3:ijkend3) :: &
178 R1,R2, R3, R4
179 DOUBLE PRECISION, DIMENSION(2) :: DOT_PRODUCT_PAR2
180 END FUNCTION DOT_PRODUCT_PAR2
181 END INTERFACE
182
183 logical, parameter :: do_unit_scaling = .true.
184
185 is_serial = numPEs.eq.1.and.is_serial
186
187 alpha(:) = zero
188 beta(:) = zero
189 omega(:) = zero
190 rho(:) = zero
191
192
193
194
195
196 if (use_doloop) then
197
198
199
200
201 do ijk2=ijkstart3,ijkend3
202 R(ijk2) = zero
203 Rtilde(ijk2) = zero
204 P(ijk2) = zero
205 Phat(ijk2) = zero
206 Svec(ijk2) = zero
207 Shat(ijk2) = zero
208 Tvec(ijk2) = zero
209 V(ijk2) = zero
210 enddo
211
212 else
213
214 R(:) = zero
215 Rtilde(:) = zero
216 P(:) = zero
217 Phat(:) = zero
218 Svec(:) = zero
219 Shat(:) = zero
220 Tvec(:) = zero
221 V(:) = zero
222
223 endif
224
225
226 TOLMIN = EPSILON( one )
227
228 if (do_unit_scaling) then
229
230
231
232
233
234
235 do k = kstart2,kend2
236 do i = istart2,iend2
237 do j = jstart2,jend2
238 IJK2 = funijk(i,j,k)
239
240 aijmax = maxval(abs(A_M(:,ijk2)) )
241
242 OAM = one/aijmax
243
244 A_M(:,IJK2) = A_M(:,IJK2)*OAM
245
246 B_M(IJK2) = B_M(IJK2)*OAM
247
248 enddo
249 enddo
250 enddo
251 endif
252
253
254
255
256
257
258
259 call MATVECt( Vname, Var, A_M, R )
260
261
262 if (use_doloop) then
263
264
265
266 do ijk2=ijkstart3,ijkend3
267 R(ijk2) = B_m(ijk2) - R(ijk2)
268 enddo
269 else
270 R(:) = B_m(:) - R(:)
271 endif
272
273 if(is_serial) then
274 Rnorm0 = zero
275 if (use_doloop) then
276
277
278
279 do ijk2=ijkstart3,ijkend3
280 Rnorm0 = Rnorm0 + R(ijk2)*R(ijk2)
281 enddo
282 else
283 Rnorm0 = dot_product(R,R)
284 endif
285 Rnorm0 = sqrt( Rnorm0 )
286 else
287 Rnorm0 = sqrt( dot_product_par( R, R ) )
288 endif
289
290 call random_number(Rtilde(:))
291
292 if (use_doloop) then
293
294
295 do ijk2=ijkstart3,ijkend3
296 Rtilde(ijk2) = R(ijk2) + (2.0d0*Rtilde(ijk2)-1.0d0)*1.0d-6*Rnorm0
297 enddo
298 else
299 Rtilde(:) = R(:) + (2.0d0*Rtilde(:)-1.0d0)*1.0d-6*Rnorm0
300 endif
301
302 if (idebugl >= 1) then
303 if(myPE.eq.0) print*,'leq_bicgs, initial: ', Vname,' resid ', Rnorm0
304 endif
305
306
307
308 = 1
309 do i=1,itmax
310
311 if(is_serial) then
312 if (use_doloop) then
313 RtildexR = zero
314
315
316 do ijk2=ijkstart3,ijkend3
317 RtildexR = RtildexR + Rtilde(ijk2) * R(ijk2)
318 enddo
319 rho(i-1) = RtildexR
320 else
321 rho(i-1) = dot_product( Rtilde, R )
322 endif
323 else
324 rho(i-1) = dot_product_par( Rtilde, R )
325 endif
326
327
328
329 if (rho(i-1) .eq. zero) then
330 if(i /= 1)then
331
332
333
334
335 = -2
336 else
337
338
339
340 = 0
341 endif
342 call send_recv(var,2)
343 return
344 endif
345
346 if (i .eq. 1) then
347 if (use_doloop) then
348
349
350 do ijk2=ijkstart3,ijkend3
351 P(ijk2) = R(ijk2)
352 enddo
353 else
354 P(:) = R(:)
355 endif
356 else
357 beta(i-1) = ( rho(i-1)/rho(i-2) )*( alpha(i-1) / omega(i-1) )
358 if (use_doloop) then
359
360
361 do ijk2=ijkstart3,ijkend3
362 P(ijk2) = R(ijk2) + beta(i-1)*( P(ijk2) - omega(i-1)*V(ijk2) )
363 enddo
364 else
365 P(:) = R(:) + beta(i-1)*( P(:) - omega(i-1)*V(:) )
366 endif
367 endif
368
369
370
371
372
373
374 call MSOLVEt( Vname, P, A_m, Phat, CMETHOD)
375
376 call MATVECt( Vname, Phat, A_m, V )
377
378 if(is_serial) then
379 if (use_doloop) then
380 RtildexV = zero
381
382
383 do ijk2=ijkstart3,ijkend3
384 RtildexV = RtildexV + Rtilde(ijk2) * V(ijk2)
385 enddo
386 else
387 RtildexV = dot_product( Rtilde, V )
388 endif
389 else
390 RtildexV = dot_product_par( Rtilde, V )
391 endif
392
393
394
395 (i) = rho(i-1) / RtildexV
396
397 if (use_doloop) then
398
399
400 do ijk2=ijkstart3,ijkend3
401 Svec(ijk2) = R(ijk2) - alpha(i) * V(ijk2)
402 enddo
403 else
404 Svec(:) = R(:) - alpha(i) * V(:)
405 endif
406
407 if(.not.minimize_dotproducts) then
408
409
410
411
412 if(is_serial) then
413 if (use_doloop) then
414 Snorm = zero
415
416
417 do ijk2=ijkstart3,ijkend3
418 Snorm = Snorm + Svec(ijk2) * Svec(ijk2)
419 enddo
420 else
421 Snorm = dot_product( Svec, Svec )
422 endif
423 Snorm = sqrt( Snorm )
424 else
425 Snorm = sqrt( dot_product_par( Svec, Svec ) )
426 endif
427
428
429
430 if (Snorm <= TOLMIN) then
431 if (use_doloop) then
432
433
434 do ijk2=ijkstart3,ijkend3
435 Var(ijk2) = Var(ijk2) + alpha(i)*Phat(ijk2)
436 enddo
437 else
438 Var(:) = Var(:) + alpha(i)*Phat(:)
439 endif
440
441 if (idebugl >= 1) then
442
443
444
445 call MATVECt( Vname, Var, A_m, R )
446
447
448
449
450 if (use_doloop) then
451
452
453 do ijk2=ijkstart3,ijkend3
454 R(ijk2) = B_m(ijk2) - R(ijk2)
455 enddo
456 else
457 R(:) = B_m(:) - R(:)
458 endif
459
460 if(is_serial) then
461 if (use_doloop) then
462 Rnorm = zero
463
464
465 do ijk2=ijkstart3,ijkend3
466 Rnorm = Rnorm + R(ijk2)*R(ijk2)
467 enddo
468 else
469 Rnorm = dot_product( R, R )
470 endif
471 Rnorm = sqrt( Rnorm )
472 else
473 Rnorm = sqrt( dot_product_par( R, R ) )
474 endif
475
476 endif
477
478 EXIT
479 endif
480
481 endif
482
483
484
485
486
487 call MSOLVEt( Vname, Svec, A_m, Shat, CMETHOD)
488
489 call MATVECt( Vname, Shat, A_m, Tvec )
490
491 if(is_serial) then
492 if (use_doloop) then
493 TxS = zero
494 TxT = zero
495
496
497 do ijk2=ijkstart3,ijkend3
498 TxS = TxS + Tvec(ijk2) * Svec(ijk2)
499 TxT = TxT + Tvec(ijk2) * Tvec(ijk2)
500 enddo
501 else
502 TxS = dot_product( Tvec, Svec )
503 TxT = dot_product( Tvec, Tvec )
504 endif
505 else
506 if(.not.minimize_dotproducts) then
507 TxS = dot_product_par( Tvec, Svec )
508 TxT = dot_product_par( Tvec, Tvec )
509 else
510 TxS_TxT = dot_product_par2(Tvec, Svec, Tvec, Tvec )
511 TxS = TxS_TxT(1)
512 TxT = TxS_TxT(2)
513 endif
514 endif
515 IF(TxT.eq.Zero) TxT = SMALL_NUMBER
516 omega(i) = TxS / TxT
517
518
519 if (use_doloop) then
520
521
522 do ijk2=ijkstart3,ijkend3
523 Var(ijk2) = Var(ijk2) + &
524 alpha(i)*Phat(ijk2) + omega(i)*Shat(ijk2)
525 R(ijk2) = Svec(ijk2) - omega(i)*Tvec(ijk2)
526 enddo
527 else
528 Var(:) = Var(:) + &
529 alpha(i)*Phat(:) + omega(i)*Shat(:)
530 R(:) = Svec(:) - omega(i)*Tvec(:)
531 endif
532
533 if(.not.minimize_dotproducts.or.(mod(iter,5).eq.0)) then
534 if(is_serial) then
535 if (use_doloop) then
536 Rnorm = zero
537
538
539 do ijk2=ijkstart3,ijkend3
540 Rnorm = Rnorm + R(ijk2) * R(ijk2)
541 enddo
542 else
543 Rnorm = dot_product(R, R )
544 endif
545 Rnorm = sqrt( Rnorm )
546 else
547 Rnorm = sqrt( dot_product_par(R, R) )
548 endif
549
550 if (idebugl.ge.1) then
551 if (myPE.eq.PE_IO) then
552 print*,'iter, Rnorm ', iter, Rnorm, Snorm
553 print*,'alpha(i), omega(i) ', alpha(i), omega(i)
554 print*,'TxS, TxT ', TxS, TxT
555 print*,'RtildexV, rho(i-1) ', RtildexV, rho(i-1)
556 endif
557 endif
558
559
560
561
562
563
564 = (Rnorm <= TOL*Rnorm0)
565
566 if (isconverged) then
567 iter_tot(vno) = iter_tot(vno) + iter + 1
568 EXIT
569 endif
570
571 endif
572
573
574 = iter + 1
575
576 enddo
577
578 if (idebugl >= 1) then
579 call MATVECt( Vname, Var, A_m, R )
580 if (use_doloop) then
581
582
583 do ijk2=ijkstart3,ijkend3
584 R(ijk2) = R(ijk2) - B_m(ijk2)
585 enddo
586 else
587 R(:) = R(:) - B_m(:)
588 endif
589
590 if(is_serial) then
591 if (use_doloop) then
592 Rnorm = zero
593
594
595 do ijk2=ijkstart3,ijkend3
596 Rnorm = Rnorm + R(ijk2) * R(ijk2)
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(myPE.eq.0) print*,'leq_bicgs: final Rnorm ', Rnorm
607
608 if(myPE.eq.0) print*,'leq_bicgs ratio : ', Vname,' ',iter, &
609 ' L-2', Rnorm/Rnorm0
610 endif
611
612 isconverged = (real(Rnorm) <= TOL*Rnorm0);
613
614 = 0
615 if (.not.isconverged) then
616 iter_tot(vno) = iter_tot(vno) + iter
617 IER = -1
618 if (real(Rnorm) >= ratiotol*real(Rnorm0)) then
619 IER = -2
620 endif
621 endif
622
623 call send_recv(var,2)
624
625 return
626 end subroutine LEQ_BICGS0t
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647 SUBROUTINE LEQ_ISWEEPt(I,Vname, VAR, A_M, B_M)
648
649
650
651
652 USE param
653 USE param1
654 USE matrix
655 USE geometry
656 USE compar
657 USE indices
658 USE funits
659 USE sendrecv
660 USE mpi_utility
661 USE functions
662 IMPLICIT NONE
663
664
665
666
667
668
669
670 INTEGER I
671
672
673
674 DOUBLE PRECISION A_m(-3:3,ijkstart3:ijkend3)
675
676
677 DOUBLE PRECISION B_m(ijkstart3:ijkend3)
678
679
680 CHARACTER(LEN=*) Vname
681
682
683
684 DOUBLE PRECISION Var(ijkstart3:ijkend3)
685
686
687
688
689
690
691
692
693
694
695
696 INTEGER :: NSTART, NEND
697 DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC,DD,EE,BB
698 INTEGER :: INFO, IJK, J, K, IM1JK, IP1JK
699
700 NEND = JEND
701 NSTART = JSTART
702 K = 1
703
704 DO J=NSTART, NEND
705
706 = FUNIJK(I,J,K)
707 IM1JK = IM_OF(IJK)
708 IP1JK = IP_OF(IJK)
709
710 DD(J) = A_M(0, IJK)
711 CC(J) = A_M(-2, IJK)
712 EE(J) = A_M(2, IJK)
713 BB(J) = B_M(IJK) - A_M(-1, IJK) * Var( IM1JK ) &
714 - A_M(1, IJK) * Var( IP1JK )
715
716 END DO
717
718 CC(NSTART) = ZERO
719 EE(NEND) = ZERO
720 INFO = 0
721
722
723 CALL DGTSV( JEND-JSTART+1, 1, CC(JSTART+1), DD, EE, BB, JEND-JSTART+1, INFO )
724
725 IF (INFO.NE.0) THEN
726 RETURN
727 ENDIF
728
729 DO J=NSTART, NEND
730 IJK = FUNIJK(I,J,K)
731 Var(IJK) = BB(J)
732 ENDDO
733
734 RETURN
735 END SUBROUTINE LEQ_ISWEEPt
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754 SUBROUTINE LEQ_IKSWEEPt(I,K,Vname, VAR, A_M, B_M )
755
756
757
758
759 USE param
760 USE param1
761 USE matrix
762 USE geometry
763 USE compar
764 USE funits
765 USE indices
766 USE sendrecv
767 USE mpi_utility
768 USE functions
769 IMPLICIT NONE
770
771
772
773
774
775
776
777 INTEGER I,K
778
779
780 DOUBLE PRECISION A_m(-3:3,ijkstart3:ijkend3)
781
782
783 DOUBLE PRECISION B_m(ijkstart3:ijkend3)
784
785
786 CHARACTER(LEN=*) Vname
787
788
789
790 DOUBLE PRECISION Var(ijkstart3:ijkend3)
791
792
793
794
795
796
797
798
799
800
801
802
803 DOUBLE PRECISION, DIMENSION (JSTART:JEND) :: CC,DD,EE, BB
804 INTEGER :: NSTART, NEND, INFO, IJK, J, IM1JK, IP1JK, IJKM1, IJKP1
805
806 NEND = JEND
807 NSTART = JSTART
808
809
810 DO J=NSTART, NEND
811
812
813 = FUNIJK(I,J,K)
814 IM1JK = IM_OF(IJK)
815 IP1JK = IP_OF(IJK)
816 IJKM1 = KM_OF(IJK)
817 IJKP1 = KP_OF(IJK)
818
819 DD(J) = A_M(0, IJK)
820 CC(J) = A_M(-2, IJK)
821 EE(J) = A_M(2, IJK)
822 BB(J) = B_M(IJK) - A_M(-1, IJK) * Var( IM1JK ) &
823 - A_M(1, IJK) * Var( IP1JK ) &
824 - A_M(-3, IJK) * Var( IJKM1 ) &
825 - A_M(3, IJK) * Var( IJKP1 )
826
827 ENDDO
828
829 CC(NSTART) = ZERO
830 EE(NEND) = ZERO
831 INFO = 0
832
833 CALL DGTSV( JEND-JSTART+1, 1, CC(JSTART+1), DD, EE, BB, JEND-JSTART+1, INFO )
834
835 IF (INFO.NE.0) THEN
836 write(*,*) 'leq_iksweep',INFO, myPE
837 RETURN
838 ENDIF
839
840 DO J=NSTART, NEND
841
842 IJK = FUNIJK(I,J,K)
843 Var(IJK) = BB(J)
844
845 ENDDO
846
847 RETURN
848 END SUBROUTINE LEQ_IKSWEEPt
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869 SUBROUTINE LEQ_MATVECt(VNAME, VAR, A_M, Avar )
870
871
872
873
874 USE param
875 USE param1
876 USE matrix
877 USE geometry
878 USE compar
879 USE indices
880 USE sendrecv
881 USE mpi_utility
882 USE functions
883 IMPLICIT NONE
884
885
886
887
888
889
890
891
892
893 DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
894
895
896 DOUBLE PRECISION AVar(ijkstart3:ijkend3)
897
898
899 CHARACTER(LEN=*) Vname
900
901
902 DOUBLE PRECISION Var(ijkstart3:ijkend3)
903
904
905
906
907
908
909
910
911
912
913 INTEGER I, J, K, IJK2
914
915 integer :: im1jk,ip1jk, ijm1k,ijp1k, ijkm1,ijkp1
916 logical, parameter :: use_send_recv = .true.
917 logical, parameter :: need_distribute_Avar = .true.
918 logical, parameter :: use_funijk = .false.
919
920 if (do_k) then
921
922
923
924
925
926
927
928 do k = kstart,kend
929 do i = istart,iend
930 do j = jstart,jend
931
932 IJK2 = funijk(i,j,k)
933
934 im1jk = im_of(ijk2)
935 ip1jk = ip_of(ijk2)
936 ijm1k = jm_of(ijk2)
937 ijp1k = jp_of(ijk2)
938
939 = km_of(ijk2)
940 ijkp1 = kp_of(ijk2)
941
942
943 AVar(ijk2) = A_m(-3, ijk2) * Var(ijkm1) &
944 + A_m(-2, ijk2) * Var(ijm1k) &
945 + A_m(-1, ijk2) * Var(im1jk) &
946 + A_m(0, ijk2) * Var(ijk2) &
947 + A_m(1, ijk2) * Var(ip1jk) &
948 + A_m(2, ijk2) * Var(ijp1k) &
949 + A_m(3, ijk2) * Var(ijkp1)
950
951 enddo
952 enddo
953 enddo
954
955 else
956 k = 1
957
958
959
960 do i = istart,iend
961 do j = jstart,jend
962
963
964 IJK2 = funijk(i,j,k)
965
966
967 im1jk = im_of(ijk2)
968 ip1jk = ip_of(ijk2)
969 ijm1k = jm_of(ijk2)
970 ijp1k = jp_of(ijk2)
971 AVar(ijk2) = A_m(-2, ijk2) * Var(ijm1k) &
972 + A_m(-1, ijk2) * Var(im1jk) &
973 + A_m(0, ijk2) * Var(ijk2) &
974 + A_m(1, ijk2) * Var(ip1jk) &
975 + A_m(2, ijk2) * Var(ijp1k)
976
977 enddo
978 enddo
979
980 endif
981
982 call send_recv(Avar,nlayers_bicgs)
983 return
984 END SUBROUTINE LEQ_MATVECt
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003 SUBROUTINE LEQ_MSOLVEt(VNAME, B_m, A_M, Var, CMETHOD)
1004
1005
1006
1007
1008 USE param
1009 USE param1
1010 USE matrix
1011 USE geometry
1012 USE compar
1013 USE indices
1014 USE sendrecv
1015 USE functions
1016 IMPLICIT NONE
1017
1018
1019
1020
1021
1022
1023
1024
1025 DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1026
1027
1028 DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1029
1030
1031 CHARACTER(LEN=*) Vname
1032
1033
1034 DOUBLE PRECISION Var(ijkstart3:ijkend3)
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046 INTEGER :: I, J, K, ITER, NITER, IJK2
1047 INTEGER :: I1 , K1 , I2, K2, IK, ISIZE, KSIZE
1048 INTEGER :: ICASE
1049
1050
1051
1052 CHARACTER(LEN=4) :: CMETHOD
1053 CHARACTER :: CH
1054 LOGICAL :: DO_ISWEEP, DO_JSWEEP, DO_KSWEEP
1055 LOGICAL :: DO_SENDRECV, DO_REDBLACK
1056 LOGICAL, PARAMETER :: USE_IKLOOP = .FALSE.
1057
1058 LOGICAL, PARAMETER :: SETGUESS = .TRUE.
1059
1060
1061
1062
1063
1064
1065 IF (SETGUESS) THEN
1066
1067
1068
1069
1070 do k = kstart3,kend3
1071 do i = istart3,iend3
1072 do j = jstart3,jend3
1073
1074 IJK2 = funijk(i,j,k)
1075
1076 VAR(IJK2) = B_M(IJK2)
1077
1078 enddo
1079 enddo
1080 enddo
1081
1082 call send_recv(var,nlayers_bicgs)
1083
1084 ENDIF
1085
1086 NITER = LEN( CMETHOD )
1087
1088 DO ITER=1,NITER
1089
1090
1091
1092 = CMETHOD( ITER:ITER )
1093 DO_ISWEEP = (CH .EQ. 'I') .OR. (CH .EQ. 'i')
1094 DO_JSWEEP = (CH .EQ. 'J') .OR. (CH .EQ. 'j')
1095 DO_KSWEEP = (CH .EQ. 'K') .OR. (CH .EQ. 'k')
1096 DO_SENDRECV = (CH .EQ. 'S') .OR. (CH .EQ. 's')
1097 DO_REDBLACK = (CH .EQ. 'R') .OR. (CH .EQ. 'r')
1098
1099 IF (NO_K) THEN
1100
1101 IF ( DO_ISWEEP ) THEN
1102
1103 DO I=istart,iend,1
1104 CALL LEQ_ISWEEPt( I, Vname, Var, A_m, B_m )
1105 ENDDO
1106 ENDIF
1107
1108 ELSE
1109
1110 IF(DO_REDBLACK) THEN
1111
1112 i1 = istart
1113 k1 = kstart
1114 i2 = iend
1115 k2 = kend
1116 isize = i2-i1+1
1117 ksize = k2-k1+1
1118
1119 DO icase = 1, 2
1120
1121 DO IK=icase, ksize*isize, 2
1122 if (mod(ik,isize).ne.0) then
1123 k = int( ik/isize ) + k1
1124 else
1125 k = int( ik/isize ) + k1 -1
1126 endif
1127 i = (ik-1-(k-k1)*isize) + i1 + mod(k,2)
1128 if(i.gt.i2) i=i-i2 + i1 -1
1129 CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1130 ENDDO
1131 ENDDO
1132
1133 ENDIF
1134
1135
1136 IF(USE_IKLOOP) THEN
1137
1138 i1 = istart
1139 k1 = kstart
1140 i2 = iend
1141 k2 = kend
1142 isize = i2-i1+1
1143 ksize = k2-k1+1
1144
1145 IF (DO_ISWEEP) THEN
1146
1147 DO IK=1, ksize*isize
1148 if (mod(ik,isize).ne.0) then
1149 k = int( ik/isize ) + k1
1150 else
1151 k = int( ik/isize ) + k1 -1
1152 endif
1153 i = (ik-1-(k-k1)*isize) + i1
1154 CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1155 ENDDO
1156 ENDIF
1157
1158 IF (DO_KSWEEP) THEN
1159
1160 DO IK=1, ksize*isize
1161 if (mod(ik,ksize).ne.0) then
1162 i = int( ik/ksize ) + i1
1163 else
1164 i = int( ik/ksize ) + i1 -1
1165 endif
1166 k = (ik-1-(i-i1)*ksize) + k1
1167
1168 CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1169 ENDDO
1170 ENDIF
1171
1172 ELSE
1173
1174
1175 IF (DO_ISWEEP) THEN
1176
1177 DO K=kstart,kend
1178 DO I=istart,iend
1179 CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1180 ENDDO
1181 ENDDO
1182 ENDIF
1183
1184 IF (DO_KSWEEP) THEN
1185
1186 DO I=istart,iend
1187 DO K=kstart,kend
1188 CALL LEQ_IKSWEEPt( I,K, Vname, Var, A_m, B_m )
1189 ENDDO
1190 ENDDO
1191 ENDIF
1192
1193 ENDIF
1194
1195 IF (DO_SENDRECV) call send_recv(var,nlayers_bicgs)
1196
1197 ENDIF
1198
1199 ENDDO
1200
1201
1202 RETURN
1203 END SUBROUTINE LEQ_MSOLVEt
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224 SUBROUTINE LEQ_JKSWEEPt(J,K,Vname, VAR, A_M, B_M )
1225
1226
1227
1228
1229 USE param
1230 USE param1
1231 USE matrix
1232 USE geometry
1233 USE funits
1234 USE compar
1235 USE indices
1236 USE functions
1237 IMPLICIT NONE
1238
1239
1240
1241
1242
1243
1244
1245 INTEGER J,K
1246
1247
1248 DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1249
1250
1251 DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1252
1253
1254 CHARACTER(LEN=*) Vname
1255
1256
1257 DOUBLE PRECISION Var(ijkstart3:ijkend3)
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269 DOUBLE PRECISION, DIMENSION (IMAX2) :: CC,DD,EE,BB
1270 INTEGER :: NN, INFO, IJK, I
1271
1272 NN = IMAX2
1273
1274 DO I=1,NN
1275 IJK = FUNIJK(I,J,K)
1276
1277 DD(I) = A_M(0, IJK)
1278 CC(I) = A_M(-1, IJK)
1279 EE(I) = A_M(1, IJK)
1280 BB(I) = B_M(IJK) - A_M(-2,IJK) * Var( JM_OF(IJK) ) &
1281 - A_M(2, IJK) * Var( JP_OF(IJK) ) &
1282 - A_M(-3, IJK) * Var( KM_OF(IJK) ) &
1283 - A_M(3, IJK) * Var( KP_OF(IJK) )
1284
1285 ENDDO
1286
1287 CC(1) = ZERO
1288 EE(NN) = ZERO
1289
1290 = 0
1291 CALL DGTSL( NN, CC, DD, EE, BB, INFO )
1292
1293
1294 IF (INFO.NE.0) THEN
1295 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
1296 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1297 call mfix_exit(myPE)
1298 ENDIF
1299
1300 DO I=1,NN
1301 IJK = FUNIJK(I,J,K)
1302 Var(IJK) = BB(I)
1303 ENDDO
1304
1305 RETURN
1306 END SUBROUTINE LEQ_JKSWEEPt
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326 SUBROUTINE LEQ_IJSWEEPt(I,J,Vname, VAR, A_M, B_M )
1327
1328
1329
1330
1331 USE param
1332 USE param1
1333 USE matrix
1334 USE geometry
1335 USE funits
1336 USE compar
1337 USE indices
1338 USE functions
1339 IMPLICIT NONE
1340
1341
1342
1343
1344
1345
1346
1347 INTEGER I,J
1348
1349
1350 DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1351
1352
1353 DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1354
1355
1356 CHARACTER(LEN=*) Vname
1357
1358
1359
1360 DOUBLE PRECISION Var(ijkstart3:ijkend3)
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371 DOUBLE PRECISION, DIMENSION (KMAX2) :: CC,DD,EE,BB
1372 INTEGER :: NN, INFO, IJK, K
1373
1374 NN = KMAX2
1375
1376 DO K=1,NN
1377 IJK = FUNIJK(I,J,K)
1378
1379 DD(K) = A_M(0,IJK)
1380 CC(K) = A_M(-3,IJK)
1381 EE(K) = A_M(3,IJK)
1382 BB(K) = B_M(IJK) - A_M(-2,IJK) * Var( JM_OF(IJK) ) &
1383 - A_M(2,IJK) * Var( JP_OF(IJK) ) &
1384 - A_M(-1,IJK) * Var( IM_OF(IJK) ) &
1385 - A_M(1,IJK) * Var( IP_OF(IJK) )
1386
1387 ENDDO
1388
1389 CC(1) = ZERO
1390 EE(NN) = ZERO
1391
1392 = 0
1393 CALL DGTSL( NN, CC, DD, EE, BB, INFO )
1394
1395
1396 IF (INFO.NE.0) THEN
1397 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'VNAME = ', VNAME
1398 IF(DMP_LOG)WRITE (UNIT_LOG,*) 'DGTSV RETURNS INFO = ', INFO
1399 call mfix_exit(myPE)
1400 ENDIF
1401
1402 DO K=1,NN
1403 IJK = FUNIJK(I,J,K)
1404 Var(IJK) = BB(K)
1405 ENDDO
1406
1407 RETURN
1408 END SUBROUTINE LEQ_IJSWEEPt
1409
1410
1411 SUBROUTINE LEQ_MSOLVE0t(VNAME, B_m, A_M, Var, CMETHOD )
1412
1413
1414
1415
1416
1417 USE param
1418 USE param1
1419 USE matrix
1420 USE geometry
1421 USE compar
1422 USE indices
1423 USE sendrecv
1424
1425 use parallel
1426 IMPLICIT NONE
1427
1428
1429
1430
1431
1432
1433
1434
1435 DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1436
1437
1438 DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1439
1440
1441 CHARACTER(LEN=*) Vname
1442
1443
1444 DOUBLE PRECISION Var(ijkstart3:ijkend3)
1445 integer :: ijk
1446
1447 CHARACTER(LEN=4) :: CMETHOD
1448
1449
1450
1451 if (use_doloop) then
1452
1453
1454 do ijk=ijkstart3,ijkend3
1455 var(ijk) = b_m(ijk)
1456 enddo
1457 else
1458 var(:) = b_m(:)
1459 endif
1460 call send_recv(var,nlayers_bicgs)
1461
1462 return
1463 end subroutine leq_msolve0t
1464
1465
1466 SUBROUTINE LEQ_msolve1t(VNAME, B_m, A_M, Var, CMETHOD )
1467
1468
1469
1470
1471
1472 USE param
1473 USE param1
1474 USE matrix
1475 USE geometry
1476 USE compar
1477 USE indices
1478 USE sendrecv
1479 USE parallel
1480 USE functions
1481
1482 IMPLICIT NONE
1483
1484
1485
1486
1487
1488
1489
1490
1491 DOUBLE PRECISION A_m(-3:3, ijkstart3:ijkend3)
1492
1493
1494 DOUBLE PRECISION B_m(ijkstart3:ijkend3)
1495
1496
1497 CHARACTER(LEN=*) Vname
1498
1499
1500 DOUBLE PRECISION Var(ijkstart3:ijkend3)
1501
1502 CHARACTER(LEN=4) :: CMETHOD
1503
1504 integer :: i,j,k, ijk2
1505
1506 if (use_doloop) then
1507
1508
1509
1510 do ijk2=ijkstart3,ijkend3
1511 var(ijk2) = zero
1512 enddo
1513 else
1514 var(:) = ZERO
1515 endif
1516
1517
1518
1519
1520
1521
1522 do k=kstart2,kend2
1523 do i=istart2,iend2
1524 do j=jstart2,jend2
1525
1526 ijk2 = funijk( i,j,k )
1527 var(ijk2) = b_m(ijk2)/A_m(0,ijk2)
1528
1529 enddo
1530 enddo
1531 enddo
1532
1533 call send_recv(var,nlayers_bicgs)
1534
1535 return
1536 end subroutine leq_msolve1t
1537
1538