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