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