MFIX  2016-1
leq_cg.f
Go to the documentation of this file.
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
compar::jend2
integer jend2
Definition: compar_mod.f:80
mpi_utility
Definition: mpi_utility_mod.f:6
funits::dmp_log
logical dmp_log
Definition: funits_mod.f:6
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
leqsol::leq_matvec
subroutine leq_matvec(VNAME, VAR, A_M, Avar)
Definition: leqsol_mod.f:104
param1::one
double precision, parameter one
Definition: param1_mod.f:29
compar::istart2
integer istart2
Definition: compar_mod.f:80
leqsol::iter_tot
integer, dimension(dim_eqs) iter_tot
Definition: leqsol_mod.f:17
compar::iend2
integer iend2
Definition: compar_mod.f:80
leqsol::leq_msolve
subroutine leq_msolve(VNAME, B_m, A_M, Var, CMETHOD)
Definition: leqsol_mod.f:287
compar::kend2
integer kend2
Definition: compar_mod.f:80
compar::kstart2
integer kstart2
Definition: compar_mod.f:80
compar::numpes
integer numpes
Definition: compar_mod.f:24
leqsol::leq_msolve0
subroutine leq_msolve0(VNAME, B_m, A_M, Var, CMETHOD)
Definition: leqsol_mod.f:617
compar::pe_io
integer pe_io
Definition: compar_mod.f:30
sendrecv::send_recv
Definition: sendrecv_mod.f:75
leqsol::dot_product_par
double precision function dot_product_par(r1, r2)
Definition: leqsol_mod.f:1095
compar::jstart2
integer jstart2
Definition: compar_mod.f:80
parallel::is_serial
logical is_serial
Definition: parallel_mod.f:11
funits::unit_log
integer, parameter unit_log
Definition: funits_mod.f:21
param
Definition: param_mod.f:2
leqsol::leq_msolve1
subroutine leq_msolve1(VNAME, B_m, A_M, Var, CMETHOD)
Definition: leqsol_mod.f:682
compar::mype
integer mype
Definition: compar_mod.f:24
parallel::use_doloop
logical use_doloop
Definition: parallel_mod.f:10
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
leq_cg0
subroutine leq_cg0(VNAME, VNO, VAR, A_M, B_m, cmethod, TOL, ITMAX, MATVEC, MSOLVE, IER)
Definition: leq_cg.f:98
leq_cg
subroutine leq_cg(VNAME, VNO, VAR, A_M, B_m, cmethod, TOL, PC, ITMAX, IER)
Definition: leq_cg.f:19
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27