File: N:\mfix\model\leq_cg.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: LEQ_CG                                                  C
4     !  Purpose: Solve system of linear system using CG method              C
5     !           conjugate gradients                                        C
6     !                                                                      C
7     !  Author: S. Pannala                                 Date: 18-JUL-07  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !  Variables referenced:                                               C
12     !  Variables modified:                                                 C
13     !  Local variables:                                                    C
14     !                                                                      C
15     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16     !
17           SUBROUTINE LEQ_CG(VNAME, VNO, VAR, A_M, B_m, cmethod, &
18                             TOL, PC, ITMAX, IER)
19     
20     !-----------------------------------------------
21     ! Modules
22     !-----------------------------------------------
23           USE param
24           USE param1
25           USE geometry
26           USE compar
27           USE indices
28           USE leqsol
29           USE funits
30           IMPLICIT NONE
31     !-----------------------------------------------
32     ! Dummy arguments
33     !-----------------------------------------------
34     ! variable name
35           CHARACTER(LEN=*), INTENT(IN) :: Vname
36     ! variable number (not really used here; see calling subroutine)
37           INTEGER, INTENT(IN) :: VNO
38     ! variable
39     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
40     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
41     !     e_Turb_G
42           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: Var
43     ! Septadiagonal matrix A_m
44           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3,-3:3), INTENT(INOUT) :: A_m
45     ! Vector b_m
46           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: B_m
47     ! Sweep direction of leq solver (leq_sweep)
48     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
49     ! Note: this setting only seems to matter when leq_pc='line'
50           CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
51     ! convergence tolerance (generally leq_tol)
52           DOUBLE PRECISION, INTENT(IN) :: TOL
53     ! preconditioner (leq_pc)
54     !     options = 'line' (default), 'diag', 'none'
55           CHARACTER(LEN=4), INTENT(IN) ::  PC
56     ! maximum number of iterations (generally leq_it)
57           INTEGER, INTENT(IN) :: ITMAX
58     ! error indicator
59           INTEGER, INTENT(INOUT) :: IER
60     !-------------------------------------------------
61     
62           if(PC.eq.'LINE') then   ! default
63              call LEQ_CG0( Vname, Vno, Var, A_m, B_m,  &
64                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE, IER )
65           elseif(PC.eq.'DIAG') then
66              call LEQ_CG0( Vname, Vno, Var, A_m, B_m,  &
67                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE1, IER )
68           elseif(PC.eq.'NONE') then
69              call LEQ_CG0( Vname, Vno, Var, A_m, B_m,  &
70                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE0, IER )
71           else
72              IF(DMP_LOG)WRITE (UNIT_LOG,*) &
73                'preconditioner option not found - check mfix.dat and readme'
74              call mfix_exit(myPE)
75           endif
76     
77           return
78           END SUBROUTINE LEQ_CG
79     
80     
81     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
82     !                                                                      C
83     !  Subroutine: LEQ_CG0                                                 C
84     !  Purpose: Compute residual of linear system                          C
85     !                                                                      C
86     !  Author: S. Pannala                                 Date: 18-JUL-07  C
87     !  Reviewer:                                          Date:            C
88     !                                                                      C
89     !  Literature/Document References:                                     C
90     !  Variables referenced:                                               C
91     !  Variables modified:                                                 C
92     !  Local variables:                                                    C
93     !                                                                      C
94     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
95     
96           SUBROUTINE LEQ_CG0(VNAME, VNO, VAR, A_M, B_m, cmethod, &
97                              TOL, ITMAX, MATVEC, MSOLVE, IER )
98     
99     !-----------------------------------------------
100     ! Modules
101     !-----------------------------------------------
102           USE param
103           USE param1
104           USE parallel
105           USE geometry
106           USE compar
107           USE mpi_utility
108           USE sendrecv
109           USE indices
110           USE leqsol
111           USE functions
112           IMPLICIT NONE
113     !-----------------------------------------------
114     ! Dummy arguments/procedure
115     !-----------------------------------------------
116     ! variable name
117           CHARACTER(LEN=*), INTENT(IN) :: Vname
118     ! variable number (not really used here-see calling subroutine)
119           INTEGER, INTENT(IN) :: VNO
120     ! variable
121     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
122     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
123     !     e_Turb_G
124           DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
125     ! Septadiagonal matrix A_m
126           DOUBLE PRECISION, INTENT(INOUT) :: A_m(ijkstart3:ijkend3,-3:3)
127     ! Vector b_m
128           DOUBLE PRECISION, INTENT(INOUT) :: B_m(ijkstart3:ijkend3)
129     ! Sweep direction of leq solver (leq_sweep)
130     !     options = 'isis', 'rsrs' (default), 'asas'
131           CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
132     ! convergence tolerance (generally leq_tol)
133           DOUBLE PRECISION, INTENT(IN) :: TOL
134     ! maximum number of iterations (generally leq_it)
135           INTEGER, INTENT(IN) :: ITMAX
136     ! error indicator
137           INTEGER, INTENT(INOUT) :: IER
138     ! dummy arguments/procedures set as indicated
139     !     matvec->leq_matvec
140     ! for preconditioner (leq_pc)
141     !    'line' msolve->leq_msolve  (default)
142     !    'diag' msolve->leq_msolve1
143     !    'none' msolve->leq_msolve0
144     !-----------------------------------------------
145     ! Local parameters
146     !-----------------------------------------------
147           INTEGER, PARAMETER :: idebugl = 0
148           DOUBLE PRECISION, PARAMETER :: ratiotol = 0.2
149           logical, parameter :: do_unit_scaling = .true.
150     !-----------------------------------------------
151     ! Local variables
152     !-----------------------------------------------
153           DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3) :: &
154                                       R, P, Zvec, Q
155           DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: &
156                                       alpha, beta, rho
157           DOUBLE PRECISION :: RxZ, PxQ, oam
158           DOUBLE PRECISION :: Rnorm, Rnorm0, TOLMIN
159           LOGICAL :: isconverged
160           INTEGER :: i, j, k, ijk
161           INTEGER :: iter
162     !-----------------------------------------------
163     
164           is_serial = numPEs.eq.1.and.is_serial
165     
166     ! these scalars should not be necessary to initialize but done as failsafe
167           rnorm = ZERO
168           rnorm0 = ZERO
169     
170     ! initializing
171           alpha(:)  = zero
172           beta(:)   = zero
173           rho(:)    = zero
174     
175     
176     ! zero out R,Zvec, P and Q
177     ! --------------------------------
178           if (use_doloop) then
179     !!$omp  parallel do private(ijk)
180              do ijk=ijkstart3,ijkend3
181                 R(ijk) = zero
182                 Zvec(ijk) = zero
183                 P(ijk) = zero
184                 Q(ijk) = zero
185              enddo
186           else
187              R(:) = zero
188              Zvec(:) = zero
189              P(:) = zero
190              Q(:) = zero
191           endif
192     
193           TOLMIN = EPSILON( one )
194     
195     ! Scale matrix to have unit diagonal
196     ! ---------------------------------------------------------------->>>
197           if (do_unit_scaling) then
198     !!$omp parallel do private(ijk,i,j,k,oam,aijmax)
199              do k = kstart2,kend2
200                 do i = istart2,iend2
201                    do j = jstart2,jend2
202                       IJK = funijk(i,j,k)
203     !                  aijmax = maxval(abs(A_M(ijk,:)) )
204     !                  if(aijmax.ne.abs(A_M(ijk,0))) &
205     !                  write(*,*) 'Not positive definite', k,i,j,(A_M(ijk,:))
206                       OAM = one/A_M(ijk,0)
207                       A_M(IJK,:) = A_M(IJK,:)*OAM
208                       B_M(IJK) = B_M(IJK)*OAM
209                    enddo
210                 enddo
211              enddo
212           endif
213     ! ----------------------------------------------------------------<<<
214     
215     
216     ! assume initial guess in Var + some small random number
217     ! r = b - A*x : Line 1
218     !     call random_number(Xinit(:))
219     !     if (use_doloop) then
220     !!$omp   parallel do private(ijk)
221     !        do ijk=ijkstart3,ijkend3
222     !           Xinit(ijk) = Var(ijk)*(ONE + (2.0d0*Xinit(ijk)-1.0d0)*1.0d-6)
223     !        enddo
224     !     else
225     !        Xinit(:) = Var(:)* (ONE + (2.0d0*Xinit(:)-1.0d0)*1.0d-6)
226     !     endif
227     !     Xinit(:) = Zero
228     
229           if (idebugl >= 1) then
230              if(myPE.eq.0) print*,'leq_cg, initial: ', Vname,' resid ', Rnorm0
231           endif
232     
233     
234     ! Compute initial residual, R = b - A*x
235     ! ---------------------------------------------------------------->>>
236           call MATVEC(Vname, Var, A_M, R)   ! returns R=A_M*Var
237     
238           if (use_doloop) then
239     !!$omp   parallel do private(ijk)
240              do ijk=ijkstart3,ijkend3
241                 R(ijk) = B_m(ijk) - R(ijk)
242              enddo
243           else
244              R(:) = B_m(:) - R(:)
245           endif
246     
247           if(is_serial) then
248              Rnorm0 = zero
249              if (use_doloop) then
250     !!$omp          parallel do private(ijk) reduction(+:Rnorm0)
251                 do ijk=ijkstart3,ijkend3
252                    Rnorm0 = Rnorm0 + R(ijk)*R(ijk)
253                 enddo
254              else
255                 Rnorm0 = dot_product(R,R)
256              endif
257              Rnorm0 = sqrt( Rnorm0 )
258           else
259              Rnorm0 = sqrt( dot_product_par( R, R ) )
260           endif
261     
262           if (idebugl >= 1) then
263              if(myPE.eq.0) print*,'leq_cg, initial: ', Vname,' resid ', Rnorm0
264           endif
265     ! ----------------------------------------------------------------<<<
266     
267     
268     ! Main loop : Line 2
269     ! ---------------------------------------------------------------->>>
270           iter = 1
271           do i=1,itmax
272     ! Solve M Zvec(:) = R(:) : Line 3
273     ! --------------------------------
274              call MSOLVE( Vname, R, A_m, Zvec, CMETHOD)   ! returns Zvec
275     
276     ! Solve Rho = RxZ : Line 4
277     ! --------------------------------
278              if(is_serial) then
279                 if (use_doloop) then
280                    RxZ = zero
281     !!$omp        parallel do private(ijk) reduction(+:RxZ)
282                    do ijk=ijkstart3,ijkend3
283                       RxZ = RxZ + R(ijk) * Zvec(ijk)
284                    enddo
285                    rho(i-1) = RxZ
286                 else
287                    rho(i-1) = dot_product( R, Zvec )
288                 endif
289              else
290                 rho(i-1) = dot_product_par( R, Zvec )
291              endif                  ! is_serial
292     
293     
294              if (rho(i-1) .eq. zero) then
295                 if(i /= 1)then
296     ! Method fails
297     ! --------------------------------
298                    ier = -2
299                 else
300     ! converged.  residual is already zero
301     ! --------------------------------
302                    ier = 0
303                 endif
304                 call send_recv(var,2)
305                 return
306              endif                  ! rho(i-1).eq.0
307     
308              if (i .eq. 1) then
309     ! P_1 = Z_0 : Line 6
310     ! --------------------------------
311                 if (use_doloop) then
312     !!$omp        parallel do private(ijk)
313                    do ijk=ijkstart3,ijkend3
314                       P(ijk) = Zvec(ijk)
315                    enddo
316                 else
317                    P(:) = Zvec(:)
318                 endif
319              else
320     ! beta = rho(i-1)/rho(i-2) : Line 8
321     ! P = Z + beta*P : Line 9
322     ! --------------------------------
323                 beta(i-1) = ( rho(i-1)/rho(i-2) )
324                 if (use_doloop) then
325     !!!$omp        parallel do private(ijk)
326                    do ijk=ijkstart3,ijkend3
327                       P(ijk) = Zvec(ijk) + beta(i-1)* P(ijk)
328                    enddo
329                 else
330                    P(:) = Zvec(:) + beta(i-1)*P(:)
331                 endif
332              endif                  ! i.eq.1
333     
334     ! Q(:) = A*P(:) : Line 10
335     ! --------------------------------
336              call MATVEC(Vname, P, A_m, Q)   ! Returns Q = A_m*P
337     
338              if(is_serial) then
339                 if (use_doloop) then
340                    PxQ = zero
341     !!$omp         parallel do private(ijk) reduction(+:PxQ)
342                    do ijk=ijkstart3,ijkend3
343                       PxQ = PxQ + P(ijk) * Q(ijk)
344                    enddo
345                 else
346                    PxQ = dot_product( P, Q )
347                 endif
348              else
349                 PxQ = dot_product_par( P, Q )
350              endif                  ! is_serial
351     
352     !  alpha = rho/PxQ : Line 11
353     ! --------------------------------
354              alpha(i) = rho(i-1)/PxQ
355     
356     ! x = x + alpha*p : Line 12
357     ! r = r - alpha*q : Line 13
358     ! --------------------------------
359              if (use_doloop) then
360     !!$omp     parallel do private(ijk)
361                 do ijk=ijkstart3,ijkend3
362                    R(ijk) = R(ijk) - alpha(i) * Q(ijk)
363                    Var(ijk) = Var(ijk) + alpha(i) * P(ijk)
364                 enddo
365              else
366                 R(:) = R(:) - alpha(i) * Q(:)
367                 Var(:) = Var(:) + alpha(i) * P(:)
368              endif                  ! use_doloop
369     
370     ! Check norm of R(:); if small enough, Exit
371     ! --------------------------------
372              if(is_serial) then
373                 if (use_doloop) then
374                    Rnorm = zero
375     !!$omp       parallel do private(ijk) reduction(+:Rnorm)
376                    do ijk=ijkstart3,ijkend3
377                       Rnorm = Rnorm + R(ijk) * R(ijk)
378                    enddo
379                 else
380                    Rnorm = dot_product( R, R )
381                 endif
382                 Rnorm = sqrt( Rnorm )
383              else
384                 Rnorm = sqrt( dot_product_par( R, R ) )
385              endif                  ! is_serial
386     
387              if (idebugl.ge.1) then
388                 print*,'leq_cs, initial: ', Vname,' Vnorm ', Rnorm
389                 if (myPE.eq.PE_IO) then
390                    print*,'iter, Rnorm ', iter, Rnorm
391                    print*,'PxQ, rho(i-1) ', PxQ, rho(i-1)
392                    print*,'alpha(i), beta(i-1) ', alpha(i), beta(i-1)
393                 endif
394              endif
395     
396              isconverged = (Rnorm <= TOL*Rnorm0)
397     
398              if (isconverged) then
399                 iter_tot(vno) = iter_tot(vno) + iter + 1
400                 EXIT
401              endif
402     
403     ! Advance the iteration count
404              iter = iter + 1
405     
406           enddo
407     ! end of linear solver loop
408     ! ----------------------------------------------------------------<<<
409     
410     
411           if (idebugl >= 1) then
412              call MATVEC( Vname, Var, A_m, R )
413              if (use_doloop) then
414     !!$omp  parallel do private(ijk)
415                 do ijk=ijkstart3,ijkend3
416                    R(ijk) = R(ijk) - B_m(ijk)
417                 enddo
418              else
419                 R(:) = R(:) - B_m(:)
420              endif
421     
422              if(is_serial) then
423                 if (use_doloop) then
424                    Rnorm = zero
425     !!$omp         parallel do private(ijk) reduction(+:Rnorm)
426                    do ijk=ijkstart3,ijkend3
427                       Rnorm = Rnorm + R(ijk) * R(ijk)
428                    enddo
429                 else
430                    Rnorm = dot_product( R,R)
431                 endif
432                 Rnorm = sqrt( Rnorm )
433              else
434                 Rnorm = sqrt( dot_product_par( R,R) )
435              endif
436     
437              if(myPE.eq.0) print*,'leq_cg: final Rnorm ', Rnorm
438     
439              if(myPE.eq.0)  print*,'leq_cg ratio : ', Vname,' ',iter,     &
440              ' L-2', Rnorm/Rnorm0
441           endif
442     
443           IER = 0
444           if (.not.isconverged) then
445              IER = -1
446              iter_tot(vno) = iter_tot(vno) + iter
447              if (real(Rnorm) >= ratiotol*real(Rnorm0)) then
448                 IER = -2
449              endif
450           endif
451     
452           call send_recv(var,2)
453     
454           return
455           end subroutine LEQ_CG0
456