File: /nfs/home/0/users/jenkins/mfix.git/model/leq_gmres.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 SUBROUTINE LEQ_GMRES(VNAME, VNO, VAR, A_M, B_M, &
18 cmethod, TOL, ITMAX, MAX_IT, IER)
19
20
21
22
23 USE PARAM
24 USE PARAM1
25 USE MATRIX
26 USE GEOMETRY
27 USE INDICES
28 USE debug
29 USE compar
30 USE mpi_utility
31 IMPLICIT NONE
32
33
34
35
36 CHARACTER(LEN=*), INTENT(IN) :: Vname
37
38 INTEGER, INTENT(IN) :: VNO
39
40
41
42
43 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: Var
44
45 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3,-3:3), INTENT(INOUT) :: A_m
46
47 DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: B_m
48
49
50 CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
51
52 DOUBLE PRECISION, INTENT(IN) :: TOL
53
54 INTEGER, INTENT(IN) :: ITMAX
55
56
57 INTEGER, INTENT(IN) :: MAX_IT
58
59 INTEGER, INTENT(INOUT) :: IER
60
61
62
63 LOGICAL :: IS_BM_ZERO, ALL_IS_BM_ZERO
64
65
66
67
68
69 EXTERNAL LEQ_MATVEC, LEQ_MSOLVE
70
71
72 = (MAXVAL( ABS(B_M(:)) ) .EQ. ZERO)
73 CALL GLOBAL_ALL_AND( IS_BM_ZERO, ALL_IS_BM_ZERO )
74
75 IF (ALL_IS_BM_ZERO) THEN
76 VAR(:) = ZERO
77 RETURN
78 ENDIF
79
80 CALL LEQ_GMRES0( VNAME, VNO, VAR, A_M, B_M, &
81 CMETHOD, TOL, ITMAX, MAX_IT, LEQ_MATVEC, LEQ_MSOLVE, IER )
82
83 RETURN
84 END SUBROUTINE LEQ_GMRES
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102 SUBROUTINE LEQ_GMRES0(VNAME, VNO, VAR, A_M, B_M, &
103 CMETHOD, TOL, ITMAX, MAX_IT, &
104 MATVEC, MSOLVE, IER )
105
106
107
108
109 USE PARAM
110 USE PARAM1
111 USE PARALLEL
112 USE MATRIX
113 USE GEOMETRY
114 USE INDICES
115 USE debug
116 USE funits
117 USE gridmap
118 USE mpi_utility
119 USE functions
120 IMPLICIT NONE
121
122
123
124 INTERFACE
125 DOUBLE PRECISION FUNCTION DOT_PRODUCT_PAR( R1, R2 )
126 use compar
127 use param
128 DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_3) :: R1, R2
129 END FUNCTION DOT_PRODUCT_PAR
130 END INTERFACE
131
132
133
134
135 CHARACTER(LEN=*), INTENT(IN) :: Vname
136
137 INTEGER, INTENT(IN) :: VNO
138
139
140
141
142 DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
143
144 DOUBLE PRECISION, INTENT(INOUT) :: A_m(ijkstart3:ijkend3,-3:3)
145
146 DOUBLE PRECISION, INTENT(INOUT) :: B_m(ijkstart3:ijkend3)
147
148
149 CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
150
151 DOUBLE PRECISION, INTENT(IN) :: TOL
152
153 INTEGER, INTENT(IN) :: ITMAX
154
155 INTEGER, INTENT(IN) :: MAX_IT
156
157 INTEGER, INTENT(INOUT) :: IER
158
159
160
161
162 EXTERNAL MATVEC, MSOLVE
163
164
165
166 INTEGER JDEBUG
167 PARAMETER(JDEBUG=0)
168 DOUBLE PRECISION TOLTRIG
169 PARAMETER(TOLTRIG=1.0D-4)
170
171
172
173
174 DOUBLE PRECISION V(IJKSTART3:IJKEND3,ITMAX+1)
175 DOUBLE PRECISION H(ITMAX+1,ITMAX)
176 DOUBLE PRECISION CS(ITMAX)
177 DOUBLE PRECISION SN(ITMAX)
178 DOUBLE PRECISION Y(ITMAX)
179
180 DOUBLE PRECISION R(IJKSTART3:IJKEND3)
181 DOUBLE PRECISION TEMP(IJKSTART3:IJKEND3)
182 DOUBLE PRECISION WW(IJKSTART3:IJKEND3)
183
184 DOUBLE PRECISION E1(ITMAX+2)
185 DOUBLE PRECISION SS(ITMAX+2)
186
187 DOUBLE PRECISION BNRM2, ERROR, ERROR0
188 DOUBLE PRECISION NORM_R0, NORM_R, NORM_W
189 DOUBLE PRECISION INV_NORM_R, NORM_S, NORM_Y
190 DOUBLE PRECISION YII
191
192 INTEGER IJK, II, JJ, KK, I, K
193 INTEGER RESTRT, M, ITER, MDIM
194 DOUBLE PRECISION DTEMP
195 DOUBLE PRECISION MINH, MAXH, CONDH, ALL_MINH, ALL_MAXH
196 DOUBLE PRECISION INV_H_IP1_I
197
198 LOGICAL IS_BM_ZERO, IS_ERROR, IS_CONVERGED
199 LOGICAL ALL_IS_BM_ZERO, ALL_IS_ERROR, ALL_IS_CONVERGED
200
201 CHARACTER(LEN=40) :: NAME
202
203
204
205
206 = 'LEQ_GMRES0 ' // TRIM(VNAME)
207 RESTRT = ITMAX
208 M = RESTRT
209 ITER = 0
210 IER = 0
211
212
213
214
215 (:) = ZERO
216 TEMP(:) = ZERO
217
218 H(:,:) = ZERO
219 CS(:) = ZERO
220 E1(:) = ZERO
221 E1(1) = ONE
222 SS(:) = ZERO
223 SN(:) = ZERO
224 WW(:) = ZERO
225 V(:,:) = ZERO
226
227
228 BNRM2 = DOT_PRODUCT_PAR( B_M, B_M )
229 BNRM2 = sqrt( BNRM2 )
230
231 if (idebug.ge.1) then
232 call write_debug(name, 'bnrm2 = ', bnrm2 )
233 endif
234
235 IS_BM_ZERO = BNRM2 .EQ. ZERO
236 CALL GLOBAL_ALL_AND( IS_BM_ZERO, ALL_IS_BM_ZERO )
237 IF (ALL_IS_BM_ZERO) THEN
238 BNRM2 = ONE
239 ENDIF
240
241
242
243
244
245 CALL MATVEC(VNAME, VAR, A_M, R)
246
247
248 DO KK=KSTART,KEND
249 DO JJ=JSTART,JEND
250 DO II=ISTART,IEND
251 IJK = FUNIJK( II,JJ,KK )
252 TEMP(IJK) = B_M(IJK) - R(IJK)
253 ENDDO
254 ENDDO
255 ENDDO
256
257
258 CALL MSOLVE(VNAME, TEMP, A_M, R, CMETHOD)
259
260 = DOT_PRODUCT_PAR( R, R )
261 NORM_R = SQRT( NORM_R )
262
263 ERROR = NORM_R/BNRM2
264 ERROR0 = ERROR
265 NORM_R0 = NORM_R
266
267 IF (JDEBUG.GE.1) THEN
268 CALL WRITE_DEBUG( NAME, ' INITIAL ERROR ', ERROR0 )
269 CALL WRITE_DEBUG( NAME, ' INITIAL RESIDUAL ', NORM_R0)
270 ENDIF
271
272
273
274 DO ITER=1,MAX_IT
275
276
277
278
279 CALL MATVEC( VNAME, VAR, A_M, R )
280
281
282
283 DO KK=KSTART,KEND
284 DO JJ=JSTART,JEND
285 DO II=ISTART,IEND
286 IJK = FUNIJK( II,JJ,KK )
287 TEMP(IJK) = B_M(IJK) - R(IJK)
288 ENDDO
289 ENDDO
290 ENDDO
291
292
293 CALL MSOLVE(VNAME, TEMP, A_M, R, CMETHOD)
294
295 = DOT_PRODUCT_PAR( R, R )
296 NORM_R = SQRT( NORM_R )
297 INV_NORM_R = ONE / NORM_R
298
299
300
301 DO KK=KSTART,KEND
302 DO JJ=JSTART,JEND
303 DO II=ISTART,IEND
304 IJK = FUNIJK( II,JJ,KK )
305 V(IJK,1) = R(IJK) * INV_NORM_R
306 ENDDO
307 ENDDO
308 ENDDO
309
310 SS(:) = NORM_R * E1(:)
311
312
313
314
315 DO I=1,M
316
317
318 CALL MATVEC(VNAME, V(:,I), A_M, TEMP)
319
320 CALL MSOLVE(VNAME, TEMP, A_M, WW, CMETHOD)
321
322 DO K=1,I
323 DTEMP = DOT_PRODUCT_PAR( WW(:), V(:,K) )
324 H(K,I) = DTEMP
325
326
327
328 DO KK=KSTART,KEND
329 DO JJ=JSTART,JEND
330 DO II=ISTART,IEND
331 IJK = FUNIJK( II,JJ,KK )
332 WW(IJK) = WW(IJK) - H(K,I)*V(IJK,K)
333 ENDDO
334 ENDDO
335 ENDDO
336 ENDDO
337
338 NORM_W = DOT_PRODUCT_PAR( WW(:), WW(:) )
339 NORM_W = SQRT( NORM_W )
340 H(I+1,I) = NORM_W
341
342 = ONE / H(I+1,I)
343
344
345 DO KK=KSTART,KEND
346 DO JJ=JSTART,JEND
347 DO II=ISTART,IEND
348 IJK = FUNIJK( II,JJ,KK )
349 V(IJK, I+1) = WW(IJK) * INV_H_IP1_I
350 ENDDO
351 ENDDO
352 ENDDO
353
354
355
356 DO K=1,I-1
357 DTEMP = CS(K)*H(K,I) + SN(K)*H(K+1,I)
358 H(K+1,I) = -SN(K)*H(K,I) + CS(K)*H(K+1,I)
359 H(K,I) = DTEMP
360 ENDDO
361
362
363
364 CALL ROTMAT( H(I,I), H(I+1,I), CS(I), SN(I) )
365
366 DTEMP = CS(I)*SS(I)
367 SS(I+1) = -SN(I)*SS(I)
368 SS(I) = DTEMP
369 H(I,I) = CS(I)*H(I,I) + SN(I)*H(I+1,I)
370 H(I+1,I) = ZERO
371 ERROR = ABS( SS(I+1) ) / BNRM2
372
373 IS_CONVERGED = (ERROR .LE. TOL*ERROR0)
374 CALL GLOBAL_ALL_AND( IS_CONVERGED, ALL_IS_CONVERGED )
375 IF (ALL_IS_CONVERGED) THEN
376
377
378
379
380
381 = I
382 DO II=1,MDIM
383 Y(II) = SS(II)
384 ENDDO
385
386 DO II=MDIM,1,-1
387 YII = Y(II)/H(II,II)
388 Y(II) = YII
389 DO JJ=1,II-1
390 Y(JJ) = Y(JJ) - H(JJ,II)*YII
391 ENDDO
392 ENDDO
393
394
395
396 IF (JDEBUG.GE.1) THEN
397 MAXH = ABS(H(1,1))
398 MINH = ABS(H(1,1))
399 DO II=1,MDIM
400 MAXH = MAX( MAXH, ABS(H(II,II)) )
401 MINH = MIN( MINH, ABS(H(II,II)) )
402 ENDDO
403
404 CALL GLOBAL_ALL_MAX( MAXH, ALL_MAXH )
405 CALL GLOBAL_ALL_MIN( MINH, ALL_MINH )
406 MAXH = ALL_MAXH
407 MINH = ALL_MINH
408 CONDH = MAXH/MINH
409
410 TEMP(1:MDIM) = SS(1:MDIM) - &
411 MATMUL( H(1:MDIM,1:MDIM ), Y(1:MDIM) )
412 DTEMP = DOT_PRODUCT( TEMP(1:MDIM), TEMP(1:MDIM) )
413 DTEMP = SQRT( DTEMP )
414
415 NORM_S = DOT_PRODUCT( SS(1:MDIM),SS(1:MDIM) )
416 NORM_S = SQRT( NORM_S )
417
418 NORM_Y = DOT_PRODUCT( Y(1:MDIM),Y(1:MDIM) )
419 NORM_Y = SQRT( NORM_Y )
420
421 IS_ERROR = (DTEMP .GT. CONDH*NORM_S)
422 CALL GLOBAL_ALL_OR( IS_ERROR, ALL_IS_ERROR )
423 IF (ALL_IS_ERROR) THEN
424 CALL WRITE_DEBUG(NAME, &
425 'DTEMP, NORM_S ', DTEMP, NORM_S )
426 CALL WRITE_DEBUG(NAME, &
427 'CONDH, NORM_Y ', CONDH, NORM_Y )
428 CALL WRITE_DEBUG(NAME, &
429 '** STOP IN LEQ_GMRES ** ')
430 IER = 999
431 RETURN
432 ENDIF
433 ENDIF
434
435
436
437
438 DO KK=KSTART,KEND
439 DO JJ=JSTART,JEND
440 DO II=ISTART,IEND
441 IJK = FUNIJK( II,JJ,KK )
442 VAR(IJK) = VAR(IJK) + &
443 DOT_PRODUCT( V(IJK,1:I), Y(1:I) )
444 ENDDO
445 ENDDO
446 ENDDO
447
448 EXIT
449 ENDIF
450 ENDDO
451
452
453 = ( ERROR .LE. TOL*ERROR0)
454 CALL GLOBAL_ALL_AND( IS_CONVERGED, ALL_IS_CONVERGED )
455 IF ( ALL_IS_CONVERGED ) THEN
456 EXIT
457 ENDIF
458
459
460
461
462
463
464
465
466 = M
467 DO II=1,MDIM
468 Y(II) = SS(II)
469 ENDDO
470
471 DO II=MDIM,1,-1
472 YII = Y(II)/H(II,II)
473 Y(II) = YII
474 DO JJ=1,II-1
475 Y(JJ) = Y(JJ) - H(JJ,II)*YII
476 ENDDO
477 ENDDO
478
479
480
481 IF (JDEBUG.GE.1) THEN
482 MAXH = ABS(H(1,1))
483 MINH = ABS(H(1,1))
484 DO II=1,MDIM
485 MAXH = MAX( MAXH, ABS(H(II,II)) )
486 MINH = MIN( MINH, ABS(H(II,II)) )
487 ENDDO
488 CONDH = MAXH/MINH
489
490 TEMP(1:MDIM) = SS(1:MDIM) - &
491 MATMUL( H(1:MDIM,1:MDIM ), Y(1:MDIM) )
492
493 DTEMP = DOT_PRODUCT( TEMP(1:MDIM), TEMP(1:MDIM) )
494 DTEMP = SQRT( DTEMP )
495
496 NORM_S = DOT_PRODUCT( SS(1:MDIM),SS(1:MDIM) )
497 NORM_S = SQRT( NORM_S )
498
499 NORM_Y = DOT_PRODUCT( Y(1:MDIM),Y(1:MDIM) )
500 NORM_Y = SQRT( NORM_Y )
501
502 IS_ERROR = (DTEMP .GT. CONDH*NORM_S)
503 CALL GLOBAL_ALL_OR( IS_ERROR, ALL_IS_ERROR )
504 IF (ALL_IS_ERROR) THEN
505 CALL WRITE_DEBUG(NAME, &
506 'DTEMP, NORM_S ', DTEMP, NORM_S )
507 CALL WRITE_DEBUG(NAME, &
508 'CONDH, NORM_Y ', CONDH, NORM_Y )
509 CALL WRITE_DEBUG(NAME, &
510 '** STOP IN LEQ_GMRES ** ')
511 IER = 999
512 RETURN
513 ENDIF
514 ENDIF
515
516
517
518 DO KK=KSTART,KEND
519 DO JJ=JSTART,JEND
520 DO II=ISTART,IEND
521 IJK = FUNIJK( II,JJ,KK )
522 VAR(IJK) = VAR(IJK) + &
523 DOT_PRODUCT( V(IJK,1:M), Y(1:M) )
524 ENDDO
525 ENDDO
526 ENDDO
527
528 CALL MATVEC(VNAME, VAR, A_M, R)
529
530
531
532 DO KK=KSTART,KEND
533 DO JJ=JSTART,JEND
534 DO II=ISTART,IEND
535 IJK = FUNIJK( II,JJ,KK )
536 TEMP(IJK) = B_M(IJK) - R(IJK)
537 ENDDO
538 ENDDO
539 ENDDO
540
541
542 CALL MSOLVE( VNAME, TEMP, A_M, R, CMETHOD)
543 = DOT_PRODUCT_PAR( R, R )
544 NORM_R = SQRT( NORM_R )
545
546 SS(I+1) = NORM_R
547 ERROR = SS(I+1) / BNRM2
548
549 IF (JDEBUG.GE.1) THEN
550 CALL WRITE_DEBUG(NAME, 'LEQ_GMRES: I, ITER ', I, ITER)
551 CALL WRITE_DEBUG(NAME, 'LEQ_GMRES: ERROR, NORM_R ', &
552 ERROR, NORM_R )
553 ENDIF
554
555 IS_CONVERGED = (ERROR .LE. TOL*ERROR0)
556 CALL GLOBAL_ALL_AND( IS_CONVERGED, ALL_IS_CONVERGED )
557 IF (ALL_IS_CONVERGED) THEN
558 EXIT
559 ENDIF
560
561 ENDDO
562
563
564 = (ERROR .GT. TOL*ERROR0)
565 CALL GLOBAL_ALL_OR( IS_ERROR, ALL_IS_ERROR )
566 IF (ALL_IS_ERROR) THEN
567 IER = 1
568 ENDIF
569
570 IF (JDEBUG.GE.1) THEN
571 CALL MATVEC(VNAME, VAR, A_M, R)
572
573
574
575 DO KK=KSTART,KEND
576 DO JJ=JSTART,JEND
577 DO II=ISTART,IEND
578 IJK = FUNIJK( II,JJ,KK )
579 R(IJK) = R(IJK) - B_M(IJK)
580 ENDDO
581 ENDDO
582 ENDDO
583
584 NORM_R=DOT_PRODUCT_PAR(R,R)
585 NORM_R = SQRT( NORM_R )
586 CALL WRITE_DEBUG(NAME, 'ITER, I ', ITER, I )
587 CALL WRITE_DEBUG(NAME, 'VNAME = ' // VNAME )
588 CALL WRITE_DEBUG(NAME, 'IER = ', IER )
589 CALL WRITE_DEBUG(NAME, 'ERROR0 = ', ERROR0 )
590 CALL WRITE_DEBUG(NAME, 'NORM_R0 = ', NORM_R0 )
591 CALL WRITE_DEBUG(NAME, 'FINAL ERROR = ', ERROR )
592 CALL WRITE_DEBUG(NAME, 'FINAL RESIDAUL ', NORM_R )
593 CALL WRITE_DEBUG(NAME, 'ERR RATIO ', ERROR/ERROR0 )
594 CALL WRITE_DEBUG(NAME, 'RESID RATIO ', NORM_R/NORM_R0 )
595 ENDIF
596
597 RETURN
598 END SUBROUTINE LEQ_GMRES0
599
600
601
602
603
604
605
606
607 SUBROUTINE ROTMAT( A, B, C, S )
608
609
610
611
612 IMPLICIT NONE
613
614
615
616 DOUBLE PRECISION A,B,C,S
617
618
619
620 DOUBLE PRECISION ONE,ZERO
621 PARAMETER(ONE=1.0D0,ZERO=0.0D0)
622
623
624
625 DOUBLE PRECISION TEMP
626
627
628 IF (B.EQ.ZERO) THEN
629 C = ONE
630 S = ZERO
631 ELSEIF (ABS(B) .GT. ABS(A)) THEN
632 TEMP = A / B
633 S = ONE / SQRT( ONE + TEMP*TEMP )
634 C = TEMP * S
635 ELSE
636 TEMP = B / A
637 C = ONE / SQRT( ONE + TEMP*TEMP )
638 S = TEMP * C
639 ENDIF
640
641 RETURN
642 END
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659 subroutine leq_check( vname, A_m )
660
661
662
663
664 USE PARAM
665 USE PARAM1
666 USE MATRIX
667 USE GEOMETRY
668 USE INDICES
669 USE debug
670 USE compar
671 USE mpi_utility
672 USE functions
673 IMPLICIT NONE
674
675
676
677
678 DOUBLE PRECISION, INTENT(INOUT) :: A_M(IJKSTART3:IJKEND3, -3:3)
679
680 CHARACTER(LEN=*), INTENT(IN) :: VNAME
681
682
683
684 logical, parameter :: do_reset = .true.
685
686
687
688 integer, dimension(-3:3) :: ijktable
689 integer :: istartl, iendl, jstartl, jendl, kstartl, kendl, elstart, elend
690 integer :: i,j,k,el, ijk,ijk2, ii,jj,kk
691 integer :: nerror, all_nerror
692 logical :: is_in_k, is_in_j, is_in_i, is_in
693 logical :: is_bc_k, is_bc_j, is_bc_i, is_bc, is_ok
694
695
696 = kstart2
697 kendl = kend2
698 jstartl = jstart2
699 jendl = jend2
700 istartl = istart2
701 iendl = iend2
702
703 if (no_k) then
704 kstartl = 1
705 kendl = 1
706 endif
707
708 nerror = 0
709
710 do k=kstartl,kendl
711 do j=jstartl,jendl
712 do i=istartl,iendl
713 is_in_k = (kstart1 <= k) .and. (k <= kend1)
714 is_in_j = (jstart1 <= j) .and. (j <= jend1)
715 is_in_i = (istart1 <= i) .and. (i <= iend1)
716
717 is_in = is_in_k .and. is_in_j .and. is_in_i
718 if (is_in) cycle
719
720 ijk = funijk(i,j,k)
721 ijktable( -2 ) = jm_of(ijk)
722 ijktable( -1 ) = im_of(ijk)
723 ijktable( 0 ) = ijk
724 ijktable( 1 ) = ip_of(ijk)
725 ijktable( 2 ) = jp_of(ijk)
726
727 if (.not. no_k) then
728 ijktable( -3 ) = km_of(ijk)
729 ijktable( 3 ) = kp_of(ijk)
730 endif
731
732 elstart = -3
733 elend = 3
734 if (no_k) then
735 elstart = -2
736 elend = 2
737 endif
738
739 do el=elstart,elend
740 ijk2 = ijktable( el )
741 kk = k_of(ijk2)
742 jj = j_of(ijk2)
743 ii = i_of(ijk2)
744 is_bc_k = (kk < kstart2) .or. (kk > kend2)
745 is_bc_j = (jj < jstart2) .or. (jj > jend2)
746 is_bc_i = (ii < istart2) .or. (ii > iend2)
747 is_bc = is_bc_k .or. is_bc_j .or. is_bc_i
748 if (is_bc) then
749 is_ok = (A_m(ijk,el).eq.zero)
750 if (.not.is_ok) then
751 nerror = nerror + 1
752 if (do_reset) A_m(ijk,el) = zero
753 endif
754 endif
755 enddo
756 enddo
757 enddo
758 enddo
759
760 call global_sum( nerror, all_nerror )
761 nerror = all_nerror
762
763 if ((nerror >= 1) .and. (myPE.eq.PE_IO)) then
764 if (do_reset) then
765 call write_debug( 'leq_check: ' // trim( vname ), &
766 'A_m is reset. nerror = ', nerror )
767 else
768 call write_debug( 'leq_check: ' // trim( vname ), &
769 'A_m is not reset. nerror = ', nerror )
770 endif
771 endif
772
773 return
774 end subroutine leq_check
775