MFIX  2016-1
conv_dif_phi.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: CONV_DIF_PHI C
4 ! Purpose: Determine convection diffusion terms for a scalar eqn. C
5 ! The off-diagonal coefficients calculated here must be positive. C
6 ! The center coefficient and the source vector are negative; C
7 ! See source_phi C
8 ! C
9 ! Diffusion at the flow boundaries is prevented by setting the C
10 ! diffusion coefficients at boundary cells to zero and then using a C
11 ! harmonic average to calculate the boundary diffusivity. The value C
12 ! diffusivities at the boundaries are checked in check_data_30. C
13 ! Ensure that harmonic avergaing is used in this routine. C
14 ! See source_phi C
15 ! C
16 ! Author: M. Syamlal Date: 21-APR-97 C
17 ! C
18 ! C
19 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20 
21  SUBROUTINE conv_dif_phi(PHI, DIF, DISC, UF, VF, WF, &
22  flux_e, flux_n, flux_t, m, a_m, b_m)
23 
24 
25 ! Modules
26 !---------------------------------------------------------------------//
27  USE param, only: dimension_3, dimension_m
28  USE run, only: def_cor
29  IMPLICIT NONE
30 
31 ! Dummy arguments
32 !---------------------------------------------------------------------//
33 ! Scalar
34  DOUBLE PRECISION Phi(dimension_3)
35 ! Gamma -- diffusion coefficient
36  DOUBLE PRECISION Dif(dimension_3)
37 ! Discretization index
38  INTEGER, INTENT(IN) :: Disc
39 ! Velocity components
40  DOUBLE PRECISION, INTENT(IN) :: Uf(dimension_3)
41  DOUBLE PRECISION, INTENT(IN) :: Vf(dimension_3)
42  DOUBLE PRECISION, INTENT(IN) :: Wf(dimension_3)
43 ! Mass flux components
44  DOUBLE PRECISION, INTENT(IN) :: Flux_E(dimension_3)
45  DOUBLE PRECISION, INTENT(IN) :: Flux_N(dimension_3)
46  DOUBLE PRECISION, INTENT(IN) :: Flux_T(dimension_3)
47 ! Phase index
48  INTEGER, INTENT(IN) :: M
49 ! Septadiagonal matrix A_m
50  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
51 ! Vector b_m
52  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
53 !---------------------------------------------------------------------//
54 
55 ! DEFERRED CORRECTION IS USED WITH THE SCALAR TRANSPORT EQN.
56  IF(def_cor)THEN
57  CALL conv_dif_phi0(phi, dif, uf, vf, wf, &
58  flux_e, flux_n, flux_t, m, a_m)
59  IF (disc > 1) CALL conv_dif_phi_dc(phi, disc, uf, vf, wf,&
60  flux_e, flux_n, flux_t, m, b_m)
61  ELSE
62 
63 ! DO NOT USE DEFERRED CORRECTION TO SOLVE THE SCALAR TRANSPORT EQN.
64  IF (disc == 0) THEN
65  CALL conv_dif_phi0(phi, dif, uf, vf, wf, &
66  flux_e, flux_n, flux_t, m, a_m)
67  ELSE
68  CALL conv_dif_phi1(phi, dif, disc, uf, vf, wf, &
69  flux_e, flux_n, flux_t, m, a_m)
70  ENDIF
71  ENDIF
72 
73  CALL dif_phi_is (dif, a_m, m)
74 
75  RETURN
76  END SUBROUTINE conv_dif_phi
77 
78 
79 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
80 ! C
81 ! Purpose: Calculate the components of diffusive flux through the C
82 ! faces of a scalar cell. Note the fluxes are calculated at C
83 ! all faces regardless of fuid_at condition of the west, south C
84 ! or bottom cell. C
85 ! C
86 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
87  SUBROUTINE get_phicell_diff_terms(Dif, D_FE, D_FW, D_FN, D_FS, &
88  d_ft, d_fb, ijk)
89 
90 ! Modules
91 !---------------------------------------------------------------------//
92  USE cutcell, only: cut_treatment_at, cut_cell_at
93 
94  USE functions, only: fluid_at
95  USE functions, only: east_of, north_of, top_of
96  USE functions, only: west_of, south_of, bottom_of
97  USE functions, only: im_of, jm_of, km_of
98  USE functions, only: ip_of, jp_of, kp_of
99 
100  USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
101 
102  USE geometry, only: odx_e, ody_n, odz_t
103  USE geometry, only: do_k
104  USE geometry, only: ox
105  USE geometry, only: dx, dy, dz
106  USE geometry, only: ayz, axz, axy
107 
108  USE indices, only: i_of, j_of, k_of
109  USE indices, only: im1, jm1, km1
110 
111  USE param, only: dimension_3
112  IMPLICIT NONE
113 
114 ! Dummy arguments
115 !---------------------------------------------------------------------//
116 ! Gamma -- diffusion coefficient
117  DOUBLE PRECISION, INTENT(IN) :: Dif(dimension_3)
118 
119 ! fluxes through faces of given ijk u-momentum cell
120  DOUBLE PRECISION, INTENT(OUT) :: d_fe, d_fw
121  DOUBLE PRECISION, INTENT(OUT) :: d_fn, d_fs
122  DOUBLE PRECISION, INTENT(OUT) :: d_ft, d_fb
123 ! ijk index
124  INTEGER, INTENT(IN) :: ijk
125 
126 ! Local variables
127 !---------------------------------------------------------------------//
128 ! indices
129  INTEGER :: imjk, ijmk, ijkm
130  INTEGER :: ipjk, ijpk, ijkp
131  INTEGER :: i, j, k, im, jm, km
132  INTEGER :: ijke, ijkw, ijkn, ijks, ijkt, ijkb
133 
134 ! area terms
135  DOUBLE PRECISION :: C_AE, C_AW, C_AN, C_AS, C_AT, C_AB
136 !---------------------------------------------------------------------//
137 
138  imjk = im_of(ijk)
139  ijmk = jm_of(ijk)
140  ijkm = km_of(ijk)
141 
142 
143  i = i_of(ijk)
144  j = j_of(ijk)
145  k = k_of(ijk)
146  im = im1(i)
147  jm = jm1(j)
148  km = km1(k)
149 
150  ijke = east_of(ijk)
151  ijkn = north_of(ijk)
152  ijks = south_of(ijk)
153  ijkw = west_of(ijk)
154 
155  c_ae = odx_e(i)*ayz(ijk)
156  c_aw = odx_e(im)*ayz(imjk)
157  c_an = ody_n(j)*axz(ijk)
158  c_as = ody_n(jm)*axz(ijmk)
159  c_at = ox(i)*odz_t(k)*axy(ijk)
160  c_ab = ox(i)*odz_t(km)*axy(ijkm)
161 
162  IF(cut_treatment_at(ijk).AND.cut_cell_at(ijk)) THEN
163  ipjk = ip_of(ijk)
164  ijpk = jp_of(ijk)
165  ijkp = kp_of(ijk)
166 
167  IF (.NOT.fluid_at(ipjk)) c_ae = odx_e(i)*dy(j)*dz(k)
168  IF (.NOT.fluid_at(imjk)) c_aw = odx_e(im)*dy(j)*dz(k)
169  IF (.NOT.fluid_at(ijpk)) c_an = ody_n(j)*dx(i)*dz(k)
170  IF (.NOT.fluid_at(ijmk)) c_as = ody_n(jm)*dx(i)*dz(k)
171  IF (.NOT.fluid_at(ijkp)) c_at = ox(i)*odz_t(k)*dx(i)*dy(j)
172  IF (.NOT.fluid_at(ijkm)) c_ab = ox(i)*odz_t(km)*dx(i)*dy(j)
173  ENDIF
174 
175 ! East face (i+1/2, j, k)
176  d_fe = avg_x_h(dif(ijk),dif(ijke),i)*c_ae
177 
178 ! West face (i-1/1, j, k)
179  d_fw = avg_x_h(dif(ijkw),dif(ijk),im)*c_aw
180 
181 ! North face (i, j+1/2, k)
182  d_fn = avg_y_h(dif(ijk),dif(ijkn),j)*c_an
183 
184 ! South face (i, j-1/2, k)
185  d_fs = avg_y_h(dif(ijks),dif(ijk),jm)*c_as
186 
187  IF (do_k) THEN
188  ijkt = top_of(ijk)
189  ijkb = bottom_of(ijk)
190 
191 ! Top face (i, j, k+1/2)
192  d_ft = avg_z_h(dif(ijk),dif(ijkt),k)*c_at
193 ! Bottom face (i, j, k-1/2)
194  d_fb = avg_z_h(dif(ijkb),dif(ijk),km)*c_ab
195  ENDIF
196 
197  RETURN
198  END SUBROUTINE get_phicell_diff_terms
199 
200 
201 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
202 ! C
203 ! Subroutine: CONV_DIF_PHI0 C
204 ! Purpose: Determine convection diffusion terms for Phi balance C
205 ! The off-diagonal coefficients calculated here must be positive. C
206 ! The center coefficient and the source vector are negative; C
207 ! See source_phi C
208 ! Implement FOUP discretization C
209 ! C
210 ! Author: M. Syamlal Date: 21-APR-97 C
211 ! C
212 ! C
213 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
214  SUBROUTINE conv_dif_phi0(PHI, DIF, UF, VF, WF, &
215  flux_e, flux_n, flux_t, m, a_m)
216 
217 ! Modules
218 !---------------------------------------------------------------------//
219  USE compar, only: ijkstart3, ijkend3
220 
221  USE functions, only: fluid_at
222  USE functions, only: ip_of, jp_of, kp_of
223  USE functions, only: im_of, jm_of, km_of
224  USE geometry, only: do_k
225  USE param
226  USE param1, only: zero
227  IMPLICIT NONE
228 
229 ! Dummy arguments
230 !---------------------------------------------------------------------//
231 ! Scalar
232  DOUBLE PRECISION, INTENT(IN) :: Phi(dimension_3)
233 ! Gamma -- diffusion coefficient
234  DOUBLE PRECISION, INTENT(IN) :: Dif(dimension_3)
235 ! Velocity components
236  DOUBLE PRECISION, INTENT(IN) :: Uf(dimension_3)
237  DOUBLE PRECISION, INTENT(IN) :: Vf(dimension_3)
238  DOUBLE PRECISION, INTENT(IN) :: Wf(dimension_3)
239 ! Mass flux components
240  DOUBLE PRECISION, INTENT(IN) :: Flux_E(dimension_3)
241  DOUBLE PRECISION, INTENT(IN) :: Flux_N(dimension_3)
242  DOUBLE PRECISION, INTENT(IN) :: Flux_T(dimension_3)
243 ! Phase index
244  INTEGER, INTENT(IN) :: M
245 ! Septadiagonal matrix A_m
246  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
247 
248 ! Local variables
249 !---------------------------------------------------------------------//
250 ! Indices
251  INTEGER :: IJK
252  INTEGER :: IPJK, IJPK, IJKP
253  INTEGER :: IMJK, IJMK, IJKM
254 ! Diffusion parameter
255  DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
256 !---------------------------------------------------------------------//
257 
258 
259 !!!$omp parallel do &
260 !!!$omp& private(IJK, IPJK, IJPK, IJKM, IMJK, IJMK, IJKM, &
261 !!!$omp& d_fe, df_w, df_n, df_s, df_t, df_b)
262  DO ijk = ijkstart3, ijkend3
263 
264  IF (fluid_at(ijk)) THEN
265 
266 ! Calculate convection-diffusion fluxes through each of the faces
267  CALL get_phicell_diff_terms(dif, d_fe, d_fw, d_fn, d_fs, &
268  d_ft, d_fb, ijk)
269 
270  ipjk = ip_of(ijk)
271  ijpk = jp_of(ijk)
272  imjk = im_of(ijk)
273  ijmk = jm_of(ijk)
274 
275 ! East face (i+1/2, j, k)
276  IF (uf(ijk) >= zero) THEN
277  a_m(ijk,east,m) = d_fe
278  a_m(ipjk,west,m) = d_fe + flux_e(ijk)
279  ELSE
280  a_m(ijk,east,m) = d_fe - flux_e(ijk)
281  a_m(ipjk,west,m) = d_fe
282  ENDIF
283 ! West face (i-1/2, j, k)
284  IF (.NOT.fluid_at(imjk)) THEN
285  IF (uf(imjk) >= zero) THEN
286  a_m(ijk,west,m) = d_fw + flux_e(imjk)
287  ELSE
288  a_m(ijk,west,m) = d_fw
289  ENDIF
290  ENDIF
291 
292 
293 ! North face (i, j+1/2, k)
294  IF (vf(ijk) >= zero) THEN
295  a_m(ijk,north,m) = d_fn
296  a_m(ijpk,south,m) = d_fn + flux_n(ijk)
297  ELSE
298  a_m(ijk,north,m) = d_fn - flux_n(ijk)
299  a_m(ijpk,south,m) = d_fn
300  ENDIF
301 ! South face (i, j-1/2, k)
302  IF (.NOT.fluid_at(ijmk)) THEN
303  IF (vf(ijmk) >= zero) THEN
304  a_m(ijk,south,m) = d_fs + flux_n(ijmk)
305  ELSE
306  a_m(ijk,south,m) = d_fs
307  ENDIF
308  ENDIF
309 
310 
311  IF (do_k) THEN
312  ijkp = kp_of(ijk)
313  ijkm = km_of(ijk)
314 
315 ! Top face (i, j, k+1/2)
316  IF (wf(ijk) >= zero) THEN
317  a_m(ijk,top,m) = d_ft
318  a_m(ijkp,bottom,m) = d_ft + flux_t(ijk)
319  ELSE
320  a_m(ijk,top,m) = d_ft - flux_t(ijk)
321  a_m(ijkp,bottom,m) = d_ft
322  ENDIF
323 
324 ! Bottom face (i, j, k-1/2)
325 
326  IF (.NOT.fluid_at(ijkm)) THEN
327  IF (wf(ijkm) >= zero) THEN
328  a_m(ijk,bottom,m) = d_fb + flux_t(ijkm)
329  ELSE
330  a_m(ijk,bottom,m) = d_fb
331  ENDIF
332  ENDIF
333  ENDIF
334 
335  ENDIF
336  ENDDO
337 
338  RETURN
339  END SUBROUTINE conv_dif_phi0
340 
341 
342 
343 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
344 ! C
345 ! Subroutine: CONV_DIF_PHI_DC C
346 ! Purpose: Use deferred correction method to solve the scalar C
347 ! transport equation. This method combines first order upwind and C
348 ! a user specified higher order method C
349 ! C
350 ! Author: C. GUENTHER Date: 1-ARP-99 C
351 ! C
352 ! C
353 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
354 
355  SUBROUTINE conv_dif_phi_dc(PHI, DISC, UF, VF, WF, &
356  flux_e, flux_n, flux_t, m, b_m)
357 
358 ! Modules
359 !---------------------------------------------------------------------//
360  USE compar, only: ijkstart3, ijkend3
361 
362  USE discretization, only: fpfoi_of
363 
364  USE function3, only: funijk3
365  USE functions, only: fluid_at
366  USE functions, only: ip_of, jp_of, kp_of
367  USE functions, only: im_of, jm_of, km_of
368 
369  USE geometry, only: do_k
370 
371  USE indices, only: i_of, j_of, k_of
372 
373  USE param, only: dimension_3, dimension_m, dimension_4
374  USE param1, only: zero, one
375 
376  USE run, only: fpfoi
377  USE sendrecv3, only: send_recv3
378 
379  USE xsi, only: calc_xsi
380  IMPLICIT NONE
381 
382 ! Dummy arguments
383 !---------------------------------------------------------------------//
384 ! Scalar
385  DOUBLE PRECISION, INTENT(IN) :: Phi(dimension_3)
386 ! Discretization index
387  INTEGER, INTENT(IN) :: Disc
388 ! Velocity components
389  DOUBLE PRECISION, INTENT(IN) :: Uf(dimension_3)
390  DOUBLE PRECISION, INTENT(IN) :: Vf(dimension_3)
391  DOUBLE PRECISION, INTENT(IN) :: Wf(dimension_3)
392 ! Mass flux components
393  DOUBLE PRECISION, INTENT(IN) :: Flux_E(dimension_3)
394  DOUBLE PRECISION, INTENT(IN) :: Flux_N(dimension_3)
395  DOUBLE PRECISION, INTENT(IN) :: Flux_T(dimension_3)
396 ! Phase index
397  INTEGER, INTENT(IN) :: M
398 ! Vector b_m
399  DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
400 
401 ! Dummy arguments
402 !---------------------------------------------------------------------//
403 ! Indices
404  INTEGER :: I, J, K, IJK
405  INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
406  INTEGER :: IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
407  INTEGER :: IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
408 ! indication for shear
409  INTEGER :: incr
410 ! Deferred correction contribution from high order method
411  DOUBLE PRECISION :: PHI_HO
412 ! low order approximation
413  DOUBLE PRECISION :: PHI_LO
414 ! deferred correction contribution from each face
415  DOUBLE PRECISION :: EAST_DC
416  DOUBLE PRECISION :: WEST_DC
417  DOUBLE PRECISION :: NORTH_DC
418  DOUBLE PRECISION :: SOUTH_DC
419  DOUBLE PRECISION :: TOP_DC
420  DOUBLE PRECISION :: BOTTOM_DC
421  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP4
422 
423  DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
424 
425 !---------------------------------------------------------------------//
426 
427  allocate(tmp4(dimension_4))
428 
429 ! Send recv the third ghost layer
430  IF (fpfoi) THEN
431  Do ijk = ijkstart3, ijkend3
432  i = i_of(ijk)
433  j = j_of(ijk)
434  k = k_of(ijk)
435  ijk4 = funijk3(i,j,k)
436  tmp4(ijk4) = phi(ijk)
437  ENDDO
438  CALL send_recv3(tmp4)
439  ENDIF
440 
441 ! shear indicator:
442  incr=0
443  CALL calc_xsi (disc, phi, uf, vf, wf, xsi_e, xsi_n, xsi_t, incr)
444 
445 !!!$omp parallel do &
446 !!!$omp& private(I, J, K, IJK, IPJK, IJPK, IJKP, &
447 !!!$omp& IMJK, IJMK, IJKM V_f, &
448 !!!$omp& PHI_HO, PHI_LO, EAST_DC, WEST_DC, NORTH_DC, &
449 !!!$omp& SOUTH_DC, TOP_DC, BOTTOM_DC)
450  DO ijk = ijkstart3, ijkend3
451 
452  IF (fluid_at(ijk)) THEN
453 
454  ipjk = ip_of(ijk)
455  imjk = im_of(ijk)
456  ijpk = jp_of(ijk)
457  ijmk = jm_of(ijk)
458  ijkp = kp_of(ijk)
459  ijkm = km_of(ijk)
460 
461 ! Third Ghost layer information
462  ippp = ip_of(ip_of(ipjk))
463  ippp4 = funijk3(i_of(ippp), j_of(ippp), k_of(ippp))
464  immm = im_of(im_of(imjk))
465  immm4 = funijk3(i_of(immm), j_of(immm), k_of(immm))
466  jppp = jp_of(jp_of(ijpk))
467  jppp4 = funijk3(i_of(jppp), j_of(jppp), k_of(jppp))
468  jmmm = jm_of(jm_of(ijmk))
469  jmmm4 = funijk3(i_of(jmmm), j_of(jmmm), k_of(jmmm))
470  kppp = kp_of(kp_of(ijkp))
471  kppp4 = funijk3(i_of(kppp), j_of(kppp), k_of(kppp))
472  kmmm = km_of(km_of(ijkm))
473  kmmm4 = funijk3(i_of(kmmm), j_of(kmmm), k_of(kmmm))
474 
475 
476 ! East face (i+1/2, j, k)
477  IF(uf(ijk)>= zero)THEN
478  phi_lo = phi(ijk)
479  IF (fpfoi) phi_ho = fpfoi_of(phi(ipjk), phi(ijk), &
480  phi(imjk), phi(im_of(imjk)))
481  ELSE
482  phi_lo = phi(ipjk)
483  IF (fpfoi) phi_ho = fpfoi_of(phi(ijk), phi(ipjk), &
484  phi(ip_of(ipjk)), tmp4(ippp4))
485  ENDIF
486  IF (.NOT. fpfoi) phi_ho = xsi_e(ijk)*phi(ipjk)+&
487  (1.0-xsi_e(ijk))*phi(ijk)
488  east_dc = flux_e(ijk)*(phi_lo - phi_ho)
489 
490 
491 ! North face (i, j+1/2, k)
492  IF(vf(ijk) >= zero)THEN
493  phi_lo = phi(ijk)
494  IF (fpfoi) phi_ho = fpfoi_of(phi(ijpk), phi(ijk), &
495  phi(ijmk), phi(jm_of(ijmk)))
496  ELSE
497  phi_lo = phi(ijpk)
498  IF (fpfoi) phi_ho = fpfoi_of(phi(ijk), phi(ijpk), &
499  phi(jp_of(ijpk)), tmp4(jppp4))
500  ENDIF
501  IF (.NOT. fpfoi) phi_ho = xsi_n(ijk)*phi(ijpk)+&
502  (1.0-xsi_n(ijk))*phi(ijk)
503  north_dc = flux_n(ijk)*(phi_lo - phi_ho)
504 
505 
506 ! Top face (i, j, k+1/2)
507  IF (do_k) THEN
508  ijkp = kp_of(ijk)
509  IF(wf(ijk) >= zero)THEN
510  phi_lo = phi(ijk)
511  IF (fpfoi) phi_ho = fpfoi_of(phi(ijkp), phi(ijk), &
512  phi(ijkm), phi(km_of(ijkm)))
513  ELSE
514  phi_lo = phi(ijkp)
515  IF (fpfoi) phi_ho = fpfoi_of(phi(ijk), phi(ijkp), &
516  phi(kp_of(ijkp)), tmp4(kppp4))
517  ENDIF
518  IF (.NOT. fpfoi) phi_ho = xsi_t(ijk)*phi(ijkp)+&
519  (1.0-xsi_t(ijk))*phi(ijk)
520  top_dc = flux_t(ijk)*(phi_lo - phi_ho)
521  ELSE
522  top_dc = zero
523  ENDIF
524 
525 
526 ! West face (i-1/2, j, k)
527  imjk = im_of(ijk)
528  IF(uf(imjk) >= zero)THEN
529  phi_lo = phi(imjk)
530  IF (fpfoi) phi_ho = fpfoi_of(phi(ijk), phi(imjk), &
531  phi(im_of(imjk)), tmp4(immm4))
532  ELSE
533  phi_lo = phi(ijk)
534  IF (fpfoi) phi_ho = fpfoi_of(phi(imjk), phi(ijk), &
535  phi(ipjk), phi(ip_of(ipjk)))
536  ENDIF
537  IF (.NOT. fpfoi) phi_ho = xsi_e(imjk)*phi(ijk)+&
538  (one-xsi_e(imjk))*phi(imjk)
539  west_dc = flux_e(imjk)*(phi_lo - phi_ho)
540 
541 
542 ! South face (i, j-1/2, k)
543  ijmk = jm_of(ijk)
544  IF(vf(ijmk) >= zero)THEN
545  phi_lo = phi(ijmk)
546  IF (fpfoi) phi_ho = fpfoi_of(phi(ijk), phi(ijmk), &
547  phi(jm_of(ijmk)), tmp4(jmmm4))
548  ELSE
549  phi_lo = phi(ijk)
550  IF (fpfoi) phi_ho = fpfoi_of(phi(ijmk), phi(ijk), &
551  phi(ijpk), phi(jp_of(ijpk)))
552  ENDIF
553  IF (.NOT. fpfoi) phi_ho = xsi_n(ijmk)*phi(ijk)+&
554  (one-xsi_n(ijmk))*phi(ijmk)
555  south_dc = flux_n(ijmk)*(phi_lo - phi_ho)
556 
557 
558 ! Bottom face (i, j, k-1/2)
559  IF (do_k) THEN
560  ijkm = km_of(ijk)
561  IF(wf(ijkm) >= zero)THEN
562  phi_lo = phi(ijkm)
563  IF (fpfoi) phi_ho = fpfoi_of(phi(ijk), phi(ijkm), &
564  phi(km_of(ijkm)), tmp4(kmmm4))
565  ELSE
566  phi_lo = phi(ijk)
567  IF (fpfoi) phi_ho = fpfoi_of(phi(ijkm), phi(ijk), &
568  phi(ijkp), phi(kp_of(ijkp)))
569  ENDIF
570  IF (.NOT. fpfoi) phi_ho = xsi_t(ijkm)*phi(ijk)+&
571  (1.0-xsi_t(ijkm))*phi(ijkm)
572  bottom_dc = flux_t(ijkm)*(phi_lo - phi_ho)
573  ELSE
574  bottom_dc = zero
575  ENDIF
576 
577 
578 ! CONTRIBUTION DUE TO DEFERRED CORRECTION
579  b_m(ijk,m) = b_m(ijk,m)+west_dc-east_dc+south_dc-&
580  north_dc+bottom_dc-top_dc
581 
582  ENDIF ! end if fluid_at
583  ENDDO ! end do ijk
584 
585  DEALLOCATE(tmp4)
586 
587  RETURN
588  END SUBROUTINE conv_dif_phi_dc
589 
590 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
591 ! C
592 ! Subroutine: CONV_DIF_Phi1 C
593 ! Purpose: Determine convection diffusion terms for scalar transport C
594 ! equations. The off-diagonal coefficients calculated here must be C
595 ! positive. The center coefficient and the source vector are C
596 ! negative; C
597 ! See source_Phi C
598 ! C
599 ! Author: M. Syamlal Date: 21-APR-97 C
600 ! C
601 ! C
602 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
603  SUBROUTINE conv_dif_phi1(PHI, DIF, DISC, UF, VF, WF, &
604  flux_e, flux_n, flux_t, m, a_m)
605 
606 ! Modules
607 !---------------------------------------------------------------------//
608  USE compar, only: ijkstart3, ijkend3
609 
610  USE functions, only: fluid_at
611  USE functions, only: ip_of, jp_of, kp_of
612  USE functions, only: im_of, jm_of, km_of
613  USE geometry, only: do_k
614  USE param
615  USE param1, only: one
616  USE xsi, only: calc_xsi
617  IMPLICIT NONE
618 
619 ! Dummy arguments
620 !---------------------------------------------------------------------//
621 ! Scalar
622  DOUBLE PRECISION, INTENT(IN) :: Phi(dimension_3)
623 ! Gamma -- diffusion coefficient
624  DOUBLE PRECISION, INTENT(IN) :: Dif(dimension_3)
625 ! Discretization index
626  INTEGER, INTENT(IN) :: Disc
627 ! Velocity components
628  DOUBLE PRECISION, INTENT(IN) :: Uf(dimension_3)
629  DOUBLE PRECISION, INTENT(IN) :: Vf(dimension_3)
630  DOUBLE PRECISION, INTENT(IN) :: Wf(dimension_3)
631 ! Mass flux components
632  DOUBLE PRECISION, INTENT(IN) :: Flux_E(dimension_3)
633  DOUBLE PRECISION, INTENT(IN) :: Flux_N(dimension_3)
634  DOUBLE PRECISION, INTENT(IN) :: Flux_T(dimension_3)
635 ! Phase index
636  INTEGER, INTENT(IN) :: M
637 ! Septadiagonal matrix A_m
638  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
639 
640 ! Local variables
641 !---------------------------------------------------------------------//
642 ! Indices
643  INTEGER :: IJK
644  INTEGER :: IPJK, IJPK, IJKP
645  INTEGER :: IMJK, IJMK, IJKM
646 ! indicator for shear
647  INTEGER :: incr
648 ! Diffusion parameter
649  DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
650 
651  DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
652 
653 !---------------------------------------------------------------------//
654 
655 ! shear indicator
656  incr=0
657  CALL calc_xsi (disc, phi, uf, vf, wf, xsi_e, xsi_n, xsi_t, incr)
658 
659 !!!$omp parallel do &
660 !!!$omp& private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM, &
661 !!!$omp& d_fe, d_fw, d_fn, d_fs, d_ft, d_fb)
662  DO ijk = ijkstart3, ijkend3
663 
664  IF (fluid_at(ijk)) THEN
665 ! Calculate convection-diffusion fluxes through each of the faces
666  CALL get_phicell_diff_terms(dif, d_fe, d_fw, d_fn, d_fs, &
667  d_ft, d_fb, ijk)
668 
669  ipjk = ip_of(ijk)
670  ijpk = jp_of(ijk)
671  imjk = im_of(ijk)
672  ijmk = jm_of(ijk)
673 
674 ! East face (i+1/2, j, k)
675  a_m(ijk,east,m) = d_fe - xsi_e(ijk)*flux_e(ijk)
676  a_m(ipjk,west,m) = d_fe + (one - xsi_e(ijk))*flux_e(ijk)
677 ! West face (i-1/2, j, k)
678  IF (.NOT.fluid_at(imjk)) THEN
679  a_m(ijk,west,m) = d_fw + (one - xsi_e(imjk))*flux_e(imjk)
680  ENDIF
681 
682 
683 ! North face (i, j+1/2, k)
684  a_m(ijk,north,m) = d_fn - xsi_n(ijk)*flux_n(ijk)
685  a_m(ijpk,south,m) = d_fn + (one - xsi_n(ijk))*flux_n(ijk)
686 ! South face (i, j-1/2, k)
687  IF (.NOT.fluid_at(ijmk)) THEN
688  a_m(ijk,south,m) = d_fs + (one - xsi_n(ijmk))*flux_n(ijmk)
689  ENDIF
690 
691 
692  IF (do_k) THEN
693  ijkp = kp_of(ijk)
694  ijkm = km_of(ijk)
695 ! Top face (i, j, k+1/2)
696  a_m(ijk,top,m) = d_ft - xsi_t(ijk)*flux_t(ijk)
697  a_m(ijkp,bottom,m)=d_ft + (one-xsi_t(ijk))*flux_t(ijk)
698 ! Bottom face (i, j, k-1/2)
699  IF (.NOT.fluid_at(ijkm)) THEN
700  a_m(ijk,bottom,m) = d_fb + (one - xsi_t(ijkm))*flux_t(ijkm)
701  ENDIF
702  ENDIF
703 
704  ENDIF
705  ENDDO
706 
707  RETURN
708  END SUBROUTINE conv_dif_phi1
709 
710 
711 
712 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
713 ! C
714 ! Subroutine: DIF_phi_IS C
715 ! Purpose: Remove diffusive fluxes across internal surfaces. C
716 ! (Make user defined internal surfaces non-conducting) C
717 ! C
718 ! Author: M. Syamlal Date: 30-APR-97 C
719 ! C
720 ! C
721 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
722  SUBROUTINE dif_phi_is(DIF, A_M, M)
723 
724 ! Modules
725 !---------------------------------------------------------------------//
726  USE param
727  USE geometry, only: do_k
728  USE geometry, only: ody_n, odx_e, odz_t, ox
729  USE geometry, only: axz, axy, ayz
730 
731  USE is, only: is_defined, is_plane
732  USE is, only: is_i_w, is_i_e, is_j_s, is_j_n, is_k_t, is_k_b
733 
734  USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
735 
736  USE functions, only: funijk, ip_of, jp_of, kp_of
737  USE functions, only: east_of, north_of, top_of
738 
739  USE compar, only: dead_cell_at
740  USE compar, only: istart2, jstart2, kstart2
741  USE compar, only: iend2, jend2, kend2
742  IMPLICIT NONE
743 
744 ! Dummy arguments
745 !---------------------------------------------------------------------//
746 ! Gamma -- diffusion coefficient
747  DOUBLE PRECISION, INTENT(IN) :: Dif(dimension_3)
748 ! Septadiagonal matrix A_m
749  DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
750 ! Solids phase
751  INTEGER, INTENT(IN) :: M
752 
753 ! Local variables
754 !---------------------------------------------------------------------//
755 ! Internal surface
756  INTEGER :: L
757 ! Indices
758  INTEGER :: I1, I2, J1, J2, K1, K2
759  INTEGER :: I, J, K, IJK
760  INTEGER :: IJKE, IJKN, IJKT, IPJK, IJPK, IJKP
761 ! Diffusion parameter
762  DOUBLE PRECISION :: D_f
763 !---------------------------------------------------------------------//
764 
765 
766  DO l = 1, dimension_is
767  IF (is_defined(l)) THEN
768  i1 = is_i_w(l)
769  i2 = is_i_e(l)
770  j1 = is_j_s(l)
771  j2 = is_j_n(l)
772  k1 = is_k_b(l)
773  k2 = is_k_t(l)
774 
775 ! Limit I1, I2 and all to local processor first ghost layer
776  IF(i1.LE.iend2) i1 = max(i1, istart2)
777  IF(j1.LE.jend2) j1 = max(j1, jstart2)
778  IF(k1.LE.kend2) k1 = max(k1, kstart2)
779  IF(i2.GE.istart2) i2 = min(i2, iend2)
780  IF(j2.GE.jstart2) j2 = min(j2, jend2)
781  IF(k2.GE.kstart2) k2 = min(k2, kend2)
782 
783  IF (is_plane(l) == 'E') THEN
784  DO k = k1, k2
785  DO j = j1, j2
786  DO i = i1, i2
787  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
788  ijk = funijk(i,j,k)
789  ijke = east_of(ijk)
790  ipjk = ip_of(ijk)
791 
792  d_f = avg_x_h(dif(ijk),dif(ijke),i)*odx_e(i)*ayz(ijk)
793  a_m(ijk,east,m) = a_m(ijk,east,m) - d_f
794  a_m(ipjk,west,m) = a_m(ipjk,west,m) - d_f
795  ENDDO
796  ENDDO
797  ENDDO
798 
799  ELSEIF(is_plane(l) == 'N') THEN
800  DO k = k1, k2
801  DO j = j1, j2
802  DO i = i1, i2
803  IF (dead_cell_at(i,j,k)) cycle ! skip dead cells
804  ijk = funijk(i,j,k)
805  ijkn = north_of(ijk)
806  ijpk = jp_of(ijk)
807 
808  d_f = avg_y_h(dif(ijk),dif(ijkn),j)*ody_n(j)*axz(ijk)
809  a_m(ijk,north,m) = a_m(ijk,north,m) - d_f
810  a_m(ijpk,south,m) = a_m(ijpk,south,m) - d_f
811  ENDDO
812  ENDDO
813  ENDDO
814 
815  ELSEIF(is_plane(l) == 'T') THEN
816  IF (do_k) THEN
817  DO k = k1, k2
818  DO j = j1, j2
819  DO i = i1, i2
820  ijkt = top_of(ijk)
821  ijkp = kp_of(ijk)
822 
823  d_f = avg_z_h(dif(ijk),dif(ijkt),k)*&
824  ox(i)*odz_t(k)*axy(ijk)
825  a_m(ijk,top,m) = a_m(ijk,top,m) - d_f
826  a_m(ijkp,bottom,m) = a_m(ijkp,bottom,m) - d_f
827  ENDDO
828  ENDDO
829  ENDDO
830  ENDIF ! end if do_k
831  ENDIF ! endif/else is_plane
832  ENDIF ! endif is_defined
833  ENDDO ! end do dimension_is
834 
835  RETURN
836  END SUBROUTINE dif_phi_is
compar::jend2
integer jend2
Definition: compar_mod.f:80
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
conv_dif_phi1
subroutine conv_dif_phi1(PHI, DIF, DISC, UF, VF, WF, Flux_E, Flux_N, Flux_T, M, A_M)
Definition: conv_dif_phi.f:605
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
param::dimension_is
integer, parameter dimension_is
Definition: param_mod.f:63
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
is::is_i_w
integer, dimension(dimension_is) is_i_w
Definition: is_mod.f:45
compar::istart2
integer istart2
Definition: compar_mod.f:80
indices::im1
integer, dimension(:), allocatable im1
Definition: indices_mod.f:50
compar::iend2
integer iend2
Definition: compar_mod.f:80
geometry::dy
double precision, dimension(0:dim_j) dy
Definition: geometry_mod.f:70
xsi::calc_xsi
subroutine calc_xsi(DISCR, PHI, U, V, W, XSI_E, XSI_N, XSI_T, incr)
Definition: xsi_mod.f:23
conv_dif_phi0
subroutine conv_dif_phi0(PHI, DIF, UF, VF, WF, Flux_E, Flux_N, Flux_T, M, A_M)
Definition: conv_dif_phi.f:216
param::east
integer east
Definition: param_mod.f:29
geometry::dz
double precision, dimension(0:dim_k) dz
Definition: geometry_mod.f:72
geometry::ayz
double precision, dimension(:), allocatable ayz
Definition: geometry_mod.f:206
compar::kend2
integer kend2
Definition: compar_mod.f:80
is::is_plane
character, dimension(dimension_is) is_plane
Definition: is_mod.f:80
is
Definition: is_mod.f:11
compar::kstart2
integer kstart2
Definition: compar_mod.f:80
xsi
Definition: xsi_mod.f:15
dif_phi_is
subroutine dif_phi_is(DIF, A_M, M)
Definition: conv_dif_phi.f:723
conv_dif_phi_dc
subroutine conv_dif_phi_dc(PHI, DISC, UF, VF, WF, Flux_E, Flux_N, Flux_T, M, B_M)
Definition: conv_dif_phi.f:357
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
discretization::fpfoi_of
double precision function fpfoi_of(PHI_D, PHI_C, PHI_U, PHI_UU)
Definition: discretization_mod.f:366
geometry::ody_n
double precision, dimension(:), allocatable ody_n
Definition: geometry_mod.f:123
fun_avg
Definition: fun_avg_mod.f:1
compar::dead_cell_at
logical, dimension(:,:,:), allocatable dead_cell_at
Definition: compar_mod.f:127
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
geometry::odx_e
double precision, dimension(:), allocatable odx_e
Definition: geometry_mod.f:121
indices::jm1
integer, dimension(:), allocatable jm1
Definition: indices_mod.f:51
is::is_k_b
integer, dimension(dimension_is) is_k_b
Definition: is_mod.f:61
compar::jstart2
integer jstart2
Definition: compar_mod.f:80
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
param::dimension_4
integer dimension_4
Definition: param_mod.f:16
param::north
integer north
Definition: param_mod.f:37
geometry::dx
double precision, dimension(0:dim_i) dx
Definition: geometry_mod.f:68
sendrecv3::send_recv3
Definition: sendrecv3_mod.f:75
param::south
integer south
Definition: param_mod.f:41
cutcell::cut_treatment_at
logical, dimension(:), allocatable cut_treatment_at
Definition: cutcell_mod.f:349
run
Definition: run_mod.f:13
geometry::axz
double precision, dimension(:), allocatable axz
Definition: geometry_mod.f:208
param
Definition: param_mod.f:2
run::def_cor
logical def_cor
Definition: run_mod.f:103
is::is_defined
logical, dimension(dimension_is) is_defined
Definition: is_mod.f:73
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
is::is_j_s
integer, dimension(dimension_is) is_j_s
Definition: is_mod.f:53
cutcell::cut_cell_at
logical, dimension(:), allocatable cut_cell_at
Definition: cutcell_mod.f:355
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
conv_dif_phi
subroutine conv_dif_phi(PHI, DIF, DISC, UF, VF, WF, Flux_E, Flux_N, Flux_T, M, A_M, B_M)
Definition: conv_dif_phi.f:23
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
param::west
integer west
Definition: param_mod.f:33
param::top
integer top
Definition: param_mod.f:45
get_phicell_diff_terms
subroutine get_phicell_diff_terms(Dif, D_FE, D_FW, D_FN, D_FS, D_FT, D_FB, IJK)
Definition: conv_dif_phi.f:89
is::is_j_n
integer, dimension(dimension_is) is_j_n
Definition: is_mod.f:57
is::is_i_e
integer, dimension(dimension_is) is_i_e
Definition: is_mod.f:49
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
function3::funijk3
integer function funijk3(LI3, LJ3, LK3)
Definition: function3_mod.f:49
run::fpfoi
logical fpfoi
Definition: run_mod.f:106
param::bottom
integer bottom
Definition: param_mod.f:49
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
is::is_k_t
integer, dimension(dimension_is) is_k_t
Definition: is_mod.f:65