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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine LEQ_BICGS                                                C
4     !  Purpose: Solve system of linear system using BICGS method           C
5     !           Biconjugate gradients stabilized                           C
6     !                                                                      C
7     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  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     ! Handan Liu wrote below:               !Jan 22 2013
17     ! The modification is as below:
18     !       Adding a loop of 2D RSRS sweep and parallelizing for OpenMP.
19     !       Splitting the existing 3D RSRS loop into two loops for OpenMP
20     !               due to data dependency
21     !       Adding openmp directives in all loops in leq_bicgs0.
22     !
23           SUBROUTINE LEQ_BICGS(VNAME, VNO, VAR, A_M, B_m, cmethod, &
24                                TOL, PC, ITMAX, IER)
25     
26     !-----------------------------------------------
27     ! Modules
28     !-----------------------------------------------
29           USE param
30           USE param1
31           USE geometry
32           USE compar
33           USE indices
34           USE leqsol
35           USE funits
36           IMPLICIT NONE
37     !-----------------------------------------------
38     ! Dummy arguments
39     !-----------------------------------------------
40     ! variable name
41           CHARACTER(LEN=*), INTENT(IN) :: Vname
42     ! variable number (not really used here; see calling subroutine)
43           INTEGER, INTENT(IN) :: VNO
44     ! variable
45     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
46     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
47     !     e_Turb_G
48     !      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: Var
49           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
50     ! Septadiagonal matrix A_m
51     !      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3,-3:3), INTENT(INOUT) :: A_m
52          DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
53     
54     ! Vector b_m
55     !      DOUBLE PRECISION, DIMENSION(ijkstart3:ijkend3), INTENT(INOUT) :: B_m
56           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
57     ! Sweep direction of leq solver (leq_sweep)
58     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
59     ! Note: this setting only seems to matter when leq_pc='line'
60           CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
61     ! convergence tolerance (generally leq_tol)
62           DOUBLE PRECISION, INTENT(IN) :: TOL
63     ! preconditioner (leq_pc)
64     !     options = 'line' (default), 'diag', 'none'
65           CHARACTER(LEN=4), INTENT(IN) ::  PC
66     ! maximum number of iterations (generally leq_it)
67           INTEGER, INTENT(IN) :: ITMAX
68     ! error indicator
69           INTEGER, INTENT(INOUT) :: IER
70     !-------------------------------------------------
71     ! Local Variables
72     !-------------------------------------------------
73     
74           if(PC.eq.'LINE') then   ! default
75              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,  &
76                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE, .true., IER )
77           elseif(PC.eq.'DIAG') then
78              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,   &
79                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE1, .true., IER )
80           elseif(PC.eq.'NONE') then
81              call LEQ_BICGS0( Vname, Vno, Var, A_m, B_m,   &
82                 cmethod, TOL, ITMAX, LEQ_MATVEC, LEQ_MSOLVE0, .false., IER )
83           else
84              IF(DMP_LOG)WRITE (UNIT_LOG,*) &
85                'preconditioner option not found - check mfix.dat and readme'
86              call mfix_exit(myPE)
87           endif
88     
89           return
90           END SUBROUTINE LEQ_BICGS
91     
92     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
93     !                                                                      C
94     !  Subroutine: LEQ_BICGS0                                              C
95     !  Purpose: Compute residual of linear system                          C
96     !                                                                      C
97     !  Author: Ed D'Azevedo                               Date: 21-JAN-99  C
98     !  Reviewer:                                          Date:            C
99     !                                                                      C
100     !  Literature/Document References:                                     C
101     !  Variables referenced:                                               C
102     !  Variables modified:                                                 C
103     !  Local variables:                                                    C
104     !                                                                      C
105     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
106     
107           SUBROUTINE LEQ_BICGS0(VNAME, VNO, VAR, A_M, B_m, cmethod, &
108                                 TOL, ITMAX, MATVEC, MSOLVE, USE_PC, IER )
109     
110     !-----------------------------------------------
111     ! Modules
112     !-----------------------------------------------
113           USE param
114           USE param1
115           USE geometry
116           USE compar
117           USE mpi_utility
118           USE sendrecv
119           USE indices
120           USE leqsol
121           USE cutcell
122           USE functions
123     
124           IMPLICIT NONE
125     !-----------------------------------------------
126     ! Dummy arguments/procedure
127     !-----------------------------------------------
128     ! variable name
129           CHARACTER(LEN=*), INTENT(IN) :: Vname
130     ! variable number (not really used here-see calling subroutine)
131           INTEGER, INTENT(IN) :: VNO
132     ! variable
133     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
134     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
135     !     e_Turb_G
136     !      DOUBLE PRECISION, INTENT(INOUT) :: Var(ijkstart3:ijkend3)
137           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: Var
138     ! Septadiagonal matrix A_m
139     !      DOUBLE PRECISION, INTENT(INOUT) :: A_m(ijkstart3:ijkend3,-3:3)
140           DOUBLE PRECISION, DIMENSION(DIMENSION_3,-3:3), INTENT(INOUT) :: A_m
141     ! Vector b_m
142     !      DOUBLE PRECISION, INTENT(INOUT) :: B_m(ijkstart3:ijkend3)
143           DOUBLE PRECISION, DIMENSION(DIMENSION_3), INTENT(INOUT) :: B_m
144     ! Sweep direction of leq solver (leq_sweep)
145     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
146           CHARACTER(LEN=*), INTENT(IN) :: CMETHOD
147     ! convergence tolerance (generally leq_tol)
148           DOUBLE PRECISION, INTENT(IN) :: TOL
149     ! maximum number of iterations (generally leq_it)
150           INTEGER, INTENT(IN) :: ITMAX
151     ! indicate whether to use preconditioner
152           LOGICAL, INTENT(IN) :: USE_PC
153     ! error indicator
154           INTEGER, INTENT(INOUT) :: IER
155     ! dummy arguments/procedures set as indicated
156     !     matvec->leq_matvec
157     ! for preconditioner (leq_pc)
158     !    'line' msolve->leq_msolve  (default)
159     !    'diag' msolve->leq_msolve1
160     !    'none' msolve->leq_msolve0
161     !-----------------------------------------------
162     ! Local parameters
163     !-----------------------------------------------
164           INTEGER, PARAMETER :: idebugl = 0
165           DOUBLE PRECISION, PARAMETER :: ratiotol = 0.2
166           LOGICAL, PARAMETER :: do_unit_scaling = .true.
167     !-----------------------------------------------
168     ! Local variables
169     !-----------------------------------------------
170     
171           DOUBLE PRECISION, DIMENSION(:), allocatable :: R,Rtilde, Tvec,V
172           DOUBLE PRECISION, DIMENSION(:), allocatable, target :: P, P_preconditioned
173           DOUBLE PRECISION, DIMENSION(:), allocatable, target :: Svec, Svec_preconditioned
174     
175           ! Phat points to either preconditioned value of P, or P itself (to avoid copying for efficiency)
176           DOUBLE PRECISION, POINTER :: Phat(:), Shat(:)
177     
178           DOUBLE PRECISION, DIMENSION(0:ITMAX+1) :: &
179                             alpha, beta, omega, rho
180           DOUBLE PRECISION :: TxS, TxT, RtildexV, &
181                               aijmax, oam
182           DOUBLE PRECISION :: Rnorm, Rnorm0, Snorm, TOLMIN, pnorm
183           LOGICAL :: isconverged
184           INTEGER :: i, j, k, ijk
185           INTEGER :: iter
186           DOUBLE PRECISION, DIMENSION(2) :: TxS_TxT
187     !-----------------------------------------------
188     
189     ! Initialize the error flag.
190           IER=0
191     
192     ! Scale matrix to have unit diagonal
193     ! ---------------------------------------------------------------->>>
194           if (do_unit_scaling) then
195     
196              IF(RE_INDEXING) THEN  ! Loop only over active cells
197     !$omp parallel do default(shared) private(ijk,oam,aijmax)
198                 DO IJK = IJKSTART3,IJKEND3
199                    aijmax = maxval(abs(A_M(ijk,:)) )
200                    if (aijmax > 0.0)then
201                       OAM = one/aijmax
202                       A_M(IJK,:) = A_M(IJK,:)*OAM
203                       B_M(IJK) = B_M(IJK)*OAM
204                    else
205                       ier = -2
206                    endif
207                 ENDDO
208     
209              ELSE
210     
211     !$omp parallel do default(shared) private(ijk,i,j,k,oam,aijmax)
212                 do k = kstart2,kend2
213                    do i = istart2,iend2
214                       do j = jstart2,jend2
215                          IJK = funijk(i,j,k)
216                          aijmax = maxval(abs(A_M(ijk,:)) )
217                          if (aijmax > 0.0) then
218                             OAM = one/aijmax
219                             A_M(IJK,:) = A_M(IJK,:)*OAM
220                             B_M(IJK) = B_M(IJK)*OAM
221                          else
222                             ier = -2
223                          endif
224                       enddo
225                    enddo
226                 enddo
227     
228              ENDIF
229     
230           endif
231     
232     ! A singlular matrix was detected.
233           if(IER /= 0) RETURN
234     ! ----------------------------------------------------------------<<<
235     
236     
237     
238           allocate(R(DIMENSION_3))
239           allocate(Rtilde(DIMENSION_3))
240           allocate(P(DIMENSION_3))
241           allocate(P_preconditioned(DIMENSION_3))
242           allocate(Svec(DIMENSION_3))
243           allocate(Svec_preconditioned(DIMENSION_3))
244           allocate(Tvec(DIMENSION_3))
245           allocate(V(DIMENSION_3))
246     
247     ! these scalars should not be necessary to initialize but done as failsafe
248           rnorm = ZERO
249           rnorm0 = ZERO
250           snorm = ZERO
251           pnorm = ZERO
252     
253     ! initializing
254           alpha(:)  = zero
255           beta(:)   = zero
256           omega(:)  = zero
257           rho(:)    = zero
258     
259     !$omp parallel sections
260           R(:) = zero
261     !$omp section
262           Rtilde(:) = zero
263     !$omp section
264           P(:) = zero
265     !$omp section
266           P_preconditioned(:) = zero
267     !$omp section
268           Svec(:) = zero
269     !$omp section
270           Svec_preconditioned(:) = zero
271     !$omp section
272           Tvec(:) = zero
273     !$omp section
274           V(:) = zero
275     !$omp end parallel sections
276     
277           TOLMIN = EPSILON( one )
278     
279     ! Compute initial residual (R = b-A*x) for Ax=b
280     !    assume initial guess in Var
281     !    rtilde = r
282     ! ---------------------------------------------------------------->>>
283           call MATVEC(Vname, Var, A_M, R)   ! returns R=A*Var
284     
285     !$omp parallel workshare
286           R(:) = B_m(:) - R(:)
287     !$omp end parallel workshare
288     
289           call send_recv(R,nlayers_bicgs)
290     
291           Rnorm0 = sqrt( dot_product_par( R, R ) )
292     
293     ! determine an initial guess for the residual = residual + small random
294     ! number (so it could be set to anything). note that since random_number
295     ! is used to supply the guess, this line could potentially be the source
296     ! of small differences between runs.  the random number is shifted below
297     ! between -1 and 1 and then scaled by factor 1.0D-6*Rnorm0
298           call random_number(Rtilde(:))
299     
300     ! Shift random number array to be consistent with case when RE_INDEXING is .FALSE.
301            IF(RE_INDEXING) CALL SHIFT_DP_ARRAY(Rtilde)
302     
303     !$omp parallel workshare
304            Rtilde(:) = R(:) + (2.0d0*Rtilde(:)-1.0d0)*1.0d-6*Rnorm0
305     !$omp end parallel workshare
306     
307           if (idebugl >= 1) then
308              if(myPE.eq.0) print*,'leq_bicgs, initial: ', Vname,' resid ', Rnorm0
309           endif
310     ! ----------------------------------------------------------------<<<
311     
312     
313     ! Main loop
314     ! ---------------------------------------------------------------->>>
315           iter = 1
316           do i=1,itmax
317     
318              rho(i-1) = dot_product_par( Rtilde, R )
319     
320              if (rho(i-1) .eq. zero) then
321                 if(i /= 1)then
322     ! Method fails
323     ! --------------------------------
324                    ier = -2
325                 else
326     ! Method converged.  residual is already zero
327     ! --------------------------------
328                    ier = 0
329                 endif
330                 call send_recv(var,2)
331                 return
332              endif ! rho(i-1).eq.0
333     
334              if (i .eq. 1) then
335     !$omp parallel workshare
336                 P(:) = R(:)
337     !$omp end parallel workshare
338              else
339                 beta(i-1) = ( rho(i-1)/rho(i-2) )*( alpha(i-1) / omega(i-1) )
340     !$omp parallel workshare
341                 P(:) = R(:) + beta(i-1)*( P(:) - omega(i-1)*V(:) )
342     !$omp end parallel workshare
343              endif ! i.eq.1
344     
345     ! Solve A*Phat(:) = P(:)
346     ! V(:) = A*Phat(:)
347     ! --------------------------------
348              if (USE_PC) then
349                 call MSOLVE(Vname, P, A_m, P_preconditioned, CMETHOD) ! returns P_preconditioned
350                 Phat => P_preconditioned
351              else
352                 Phat => P
353              endif
354     
355              call MATVEC(Vname, Phat, A_m, V)   ! returns V=A*Phat
356     
357              RtildexV = dot_product_par( Rtilde, V )
358     
359     ! compute alpha
360     ! --------------------------------
361              alpha(i) = rho(i-1) / RtildexV
362     
363     ! compute Svec
364     ! --------------------------------
365     !$omp parallel workshare
366              Svec(:) = R(:) - alpha(i) * V(:)
367     !$omp end parallel workshare
368     
369     ! Check norm of Svec(:); if small enough:
370     ! set X(:) = X(:) + alpha(i)*Phat(:) and stop
371     ! --------------------------------
372              if(.not.minimize_dotproducts) then
373                 Snorm = sqrt( dot_product_par( Svec, Svec ) )
374     
375                 if (Snorm <= TOLMIN) then
376     !$omp parallel workshare
377                    Var(:) = Var(:) + alpha(i)*Phat(:)
378     !$omp end parallel workshare
379     
380     ! Recompute residual norm
381     ! --------------------------------
382                    if (idebugl >= 1) then
383                       call MATVEC(Vname, Var, A_m, R)   ! returns R=A*Var
384     !                  Rnorm = sqrt( dot_product_par( Var, Var ) )
385     !                  print*,'leq_bicgs, initial: ', Vname,' Vnorm ', Rnorm
386     
387     !$omp parallel workshare
388                       R(:) = B_m(:) - R(:)
389     !$omp end parallel workshare
390     
391                       Rnorm = sqrt( dot_product_par( R, R ) )
392                    endif            ! idebugl >= 1
393                    isConverged = .TRUE.
394                    EXIT
395                 endif               ! end if (Snorm <= TOLMIN)
396              endif                  ! end if (.not.minimize_dotproducts)
397     
398     ! Solve A*Shat(:) = Svec(:)
399     ! Tvec(:) = A*Shat(:)
400     ! --------------------------------
401     
402              if (USE_PC) then
403                 call MSOLVE(Vname, Svec, A_m, Svec_preconditioned, CMETHOD) ! returns S_preconditioned
404                 Shat => Svec_preconditioned
405              else
406                 Shat => Svec
407              endif
408     
409              call MATVEC( Vname, Shat, A_m, Tvec )   ! returns Tvec=A*Shat
410     
411              if(.not.minimize_dotproducts) then
412     !!     $omp parallel sections
413                 TxS = dot_product_par( Tvec, Svec )
414     !!     $omp section
415                 TxT = dot_product_par( Tvec, Tvec )
416     !!     $omp end parallel sections
417              else
418                 TxS_TxT = dot_product_par2(Tvec, Svec, Tvec, Tvec )
419                 TxS = TxS_TxT(1)
420                 TxT = TxS_TxT(2)
421              endif
422     
423              IF(TxT.eq.Zero) TxT = SMALL_NUMBER
424     
425     ! compute omega
426     ! --------------------------------
427              omega(i) = TxS / TxT
428     
429     ! compute new guess for Var
430     ! --------------------------------
431     
432     !$omp parallel sections
433                 Var(:) = Var(:) + alpha(i)*Phat(:) + omega(i)*Shat(:)
434     !$omp section
435                 R(:) = Svec(:) - omega(i)*Tvec(:)
436     !$omp end parallel sections
437     
438     ! --------------------------------
439              if(.not.minimize_dotproducts.or.(mod(iter,icheck_bicgs).eq.0)) then
440                 Rnorm = sqrt( dot_product_par(R, R) )
441     
442                 if (idebugl.ge.1) then
443                    if (myPE.eq.PE_IO) then
444                       print*,'iter, Rnorm ', iter, Rnorm, Snorm
445                       print*,'alpha(i), omega(i) ', alpha(i), omega(i)
446                       print*,'TxS, TxT ', TxS, TxT
447                       print*,'RtildexV, rho(i-1) ', RtildexV, rho(i-1)
448                    endif
449                 endif
450     
451     !           call mfix_exit(myPE)
452     
453     ! Check convergence; continue if necessary
454     ! for continuation, it is necessary that omega(i) .ne. 0
455                 isconverged = (Rnorm <= TOL*Rnorm0)
456     
457                 if (isconverged) then
458                    iter_tot(vno) = iter_tot(vno) + iter + 1
459                    EXIT
460                 endif
461              endif                  ! end if(.not.minimize_dotproducts)
462     
463     ! Advance the iteration count
464              iter = iter + 1
465     
466           enddo   ! end do i=1,itmax
467     ! end of linear solver loop
468     ! ----------------------------------------------------------------<<<
469     
470           if (idebugl >= 1) then
471              call MATVEC(Vname, Var, A_m, R)   ! returns R=A*Var
472     
473     !$omp parallel workshare
474              R(:) = R(:) - B_m(:)
475     !$omp end parallel workshare
476     
477              Rnorm = sqrt( dot_product_par( R,R) )
478     
479              if(myPE.eq.0) print*,'leq_bicgs: final Rnorm ', Rnorm
480     
481              if(myPE.eq.0)  print*,'leq_bicgs ratio : ', Vname,' ',iter,  &
482              ' L-2', Rnorm/Rnorm0
483           endif   ! end if(idebugl >=1)
484     
485     !      isconverged = (real(Rnorm) <= TOL*Rnorm0);
486           if(.NOT.isConverged) isconverged = (real(Rnorm) <= TOL*Rnorm0);
487     !     write(*,*) '***',iter, isconverged, Rnorm, TOL, Rnorm0, myPE
488           IER = 0
489           if (.not.isconverged) then
490              IER = -1
491              iter_tot(vno) = iter_tot(vno) + iter
492              if (real(Rnorm) >= ratiotol*real(Rnorm0)) then
493                 IER = -2
494              endif
495           endif
496     
497           call send_recv(var,2)
498     
499           deallocate(R)
500           deallocate(Rtilde)
501           deallocate(P)
502           deallocate(P_preconditioned)
503           deallocate(Svec)
504           deallocate(Svec_preconditioned)
505           deallocate(Tvec)
506           deallocate(V)
507     
508           return
509           end subroutine LEQ_BICGS0
510