MFIX  2016-1
discretization_mod.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: DISCRETIZATION C
4 ! Purpose: A collection of functions for higher-order discretization C
5 ! C
6 ! Author: M. Syamlal Date: 18-MAR-97 C
7 ! C
8 ! C
9 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10 
11  MODULE discretization
12  IMPLICIT NONE
13 
14  CONTAINS
15 
16 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
17 ! !
18 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
19  DOUBLE PRECISION FUNCTION superbee (PHI_C)
20  USE param1, only: one, half, zero
21  IMPLICIT NONE
22  DOUBLE PRECISION, INTENT(IN) :: PHI_C
23 
24  DOUBLE PRECISION :: TH
25 
26  IF (phi_c>=zero .AND. phi_c<one) THEN !monotonic region
27  th = phi_c/(one - phi_c)
28  superbee = half*max(zero,min(one,2.d0*th),min(2.d0,th))
29  ELSE IF (phi_c == one) THEN
30  superbee = one
31  ELSE !first order upwinding
32  superbee = zero
33  ENDIF
34  RETURN
35  END FUNCTION superbee
36 
37 
38 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
39 ! !
40 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
41  DOUBLE PRECISION FUNCTION smart (PHI_C)
42  USE param1, only: zero, half, one
43  IMPLICIT NONE
44  DOUBLE PRECISION, INTENT(IN) :: PHI_C
45 
46  DOUBLE PRECISION :: TH
47 
48  IF (phi_c>=zero .AND. phi_c<one) THEN !monotonic region
49  th = phi_c/(one - phi_c)
50  smart = half*max(zero,min(4.0d0*th,0.75d0 + 0.25d0*th,2.0d0))
51  ELSE IF (phi_c == one) THEN
52  smart = one
53  ELSE !first order upwinding
54  smart = zero
55  ENDIF
56  RETURN
57  END FUNCTION smart
58 
59 
60 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
61 ! !
62 ! calculate DWF from Chi_SMART scheme !
63 ! !
64 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
65  DOUBLE PRECISION FUNCTION chi_smart (PHI_C, Chi)
66  USE param1, only: zero, one
67  IMPLICIT NONE
68  DOUBLE PRECISION, INTENT(IN) :: PHI_C
69  DOUBLE PRECISION, INTENT(IN) :: Chi
70 
71  if(phi_c < one)then
72  chi_smart = chi * (3.d0/8.d0 - phi_c/4.d0) / (one-phi_c)
73  elseif(chi == zero)then ! insures that all species equations uses the same chi
74  chi_smart = zero
75  else
76  chi_smart = one
77  endif
78  RETURN
79  END FUNCTION chi_smart
80 
81 
82 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
83 ! !
84 ! calculate CHI for SMART scheme C
85 ! !
86 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
87  DOUBLE PRECISION FUNCTION chi4smart (PHI_C, PHIU, PHIC, PHID)
88  USE param1, only: zero, one, large_number
89  IMPLICIT NONE
90  DOUBLE PRECISION, INTENT(IN) :: PHI_C
91  DOUBLE PRECISION, INTENT(IN) :: PHIU
92  DOUBLE PRECISION, INTENT(IN) :: PHIC
93  DOUBLE PRECISION, INTENT(IN) :: PHID
94 
95  IF(phi_c > zero .AND. phi_c <= (1.d0/6.d0))THEN
96  chi4smart = 16.d0 * phi_c/(3.d0 - 2.d0 * phi_c)
97  ELSEIF(phi_c > (1.d0/6.d0) .AND. phi_c <= (5.d0/6.d0))THEN
98  chi4smart = one
99  ELSEIF(phi_c > (5.d0/6.d0) .AND. phi_c <= one)THEN
100  chi4smart = 8.d0*(one-phi_c)/(3.d0 - 2.d0 * phi_c)
101  ELSEIF( (phiu < 1d-15) .AND. (phic < 1d-15) .AND. (phid < 1d-15) )THEN
102 ! if a species Xg is less than machine precision, do not use its chi.
103 ! this will avoid calculating chi = 0 for a vanishing species.
104  chi4smart = large_number
105  ELSE
106  chi4smart = zero
107  ENDIF
108  RETURN
109  END FUNCTION chi4smart
110 
111 
112 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
113 ! !
114 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
115  DOUBLE PRECISION FUNCTION ultra_quick (PHI_C, CF)
116  USE param1, only: zero, half, one
117  IMPLICIT NONE
118  DOUBLE PRECISION, INTENT(IN) :: PHI_C
119  DOUBLE PRECISION, INTENT(IN) :: CF
120 
121  DOUBLE PRECISION, PARAMETER :: FIVEOSIX = 5.d0/6.d0
122  DOUBLE PRECISION :: TH, OCF
123 
124  ocf = max(one,one/max(1.d-2,cf))
125  IF (phi_c > one) THEN
126  ultra_quick = half
127  ELSE IF (phi_c > fiveosix) THEN
128  ultra_quick = one
129  ELSE IF (phi_c > 3.d0/(8.d0*ocf - 6.d0)) THEN
130  th = phi_c/(one - phi_c)
131  ultra_quick = half - 0.125d0*(one - th)
132  ELSE IF (phi_c > zero) THEN
133  th = phi_c/(one - phi_c)
134  ultra_quick = (ocf - one)*th
135  ELSE
136  th = phi_c/(one - phi_c)
137  ultra_quick = half*th
138  ENDIF
139  RETURN
140  END FUNCTION ultra_quick
141 
142 
143 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
144 ! !
145 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
146  DOUBLE PRECISION FUNCTION quickest (PHI_C, CF, ODXC, ODXUC, ODXCD)
147  USE param1, only: zero, half, one, small_number
148  IMPLICIT NONE
149  DOUBLE PRECISION, INTENT(IN) :: PHI_C
150  DOUBLE PRECISION, INTENT(IN) :: CF
151  DOUBLE PRECISION, INTENT(IN) :: ODXC
152  DOUBLE PRECISION, INTENT(IN) :: ODXUC
153  DOUBLE PRECISION, INTENT(IN) :: ODXCD
154 
155  DOUBLE PRECISION :: FCF, TH
156 
157  IF (phi_c>zero .AND. phi_c<one) THEN !monotonic region
158  fcf = -(one - cf*cf)/3.d0
159  th = phi_c/(one - phi_c + small_number)
160  quickest = half*(one - cf) + fcf*(odxc/odxcd - odxc*odxuc*th/odxcd**2)
161  IF(phi_c<cf)quickest=min(quickest,(one/cf-one)*phi_c/(one-phi_c))
162  quickest = max(zero,min(one,quickest))
163  ELSE !first order upwinding
164  quickest = zero
165  ENDIF
166  RETURN
167  END FUNCTION quickest
168 
169 
170 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
171 ! !
172 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
173  DOUBLE PRECISION FUNCTION muscl (PHI_C)
174  USE param1, only: zero, half, one
175  IMPLICIT NONE
176  DOUBLE PRECISION, INTENT(IN) :: PHI_C
177 
178  DOUBLE PRECISION :: TH
179 
180  IF (phi_c>=zero .AND. phi_c<one) THEN !monotonic region
181  th = phi_c/(one - phi_c)
182  muscl = half*max(zero,min(2.0d0*th,(one + th)/2.0d0,2.0d0))
183  ELSE IF (phi_c == one) THEN
184  muscl = one
185  ELSE !first order upwinding
186  muscl = zero
187  ENDIF
188  RETURN
189  END FUNCTION muscl
190 
191 
192 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
193 ! !
194 ! calculate DWF from Chi_MUSCL scheme !
195 ! !
196 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
197  DOUBLE PRECISION FUNCTION chi_muscl (PHI_C, Chi)
198  USE param1, only: zero, one
199  IMPLICIT NONE
200  DOUBLE PRECISION, INTENT(IN) :: PHI_C
201  DOUBLE PRECISION, INTENT(IN) :: Chi
202 
203  if(phi_c < one)then
204  chi_muscl = chi /(4.d0 * (one - phi_c))
205  elseif(chi == zero)then ! insures that all species equations uses the same chi
206  chi_muscl = zero
207  else
208  chi_muscl = one
209  endif
210  RETURN
211  END FUNCTION chi_muscl
212 
213 
214 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
215 ! !
216 ! calculate CHI for MUSCL scheme !
217 ! !
218 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
219  DOUBLE PRECISION FUNCTION chi4muscl (PHI_C, PHIU, PHIC, PHID)
220  USE param1, only: zero, one, large_number
221  IMPLICIT NONE
222  DOUBLE PRECISION, INTENT(IN) :: PHI_C
223  DOUBLE PRECISION, INTENT(IN) :: PHIU
224  DOUBLE PRECISION, INTENT(IN) :: PHIC
225  DOUBLE PRECISION, INTENT(IN) :: PHID
226 
227  IF(phi_c > zero .AND. phi_c <= 0.25d0)THEN
228  chi4muscl = 4.d0 * phi_c
229  ELSEIF(phi_c > 0.25d0 .AND. phi_c <= 0.75d0)THEN
230  chi4muscl = one
231  ELSEIF(phi_c > 0.75d0 .AND. phi_c <= one)THEN
232  chi4muscl = 4.d0*(one - phi_c)
233  ELSEIF( (phiu < 1d-15) .AND. (phic < 1d-15) .AND. (phid < 1d-15) )THEN
234 ! if a species Xg is less than machine precision, do not use its chi.
235 ! this will avoid calculating chi = 0 for a vanishing species.
236  chi4muscl = large_number
237  ELSE
238  chi4muscl = zero
239  ENDIF
240  RETURN
241  END FUNCTION chi4muscl
242 
243 
244 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
245 ! !
246 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
247  DOUBLE PRECISION FUNCTION vanleer (PHI_C)
248  USE param1, only: zero, one
249  IMPLICIT NONE
250  DOUBLE PRECISION, INTENT(IN) :: PHI_C
251 
252  IF (phi_c>=zero .AND. phi_c<=one) THEN !monotonic region
253  vanleer = phi_c
254  ELSE !first order upwinding
255  vanleer = zero
256  ENDIF
257  RETURN
258  END FUNCTION vanleer
259 
260 
261 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
262 ! !
263 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
264  DOUBLE PRECISION FUNCTION minmod (PHI_C)
265  USE param1, only: zero, half, one
266  IMPLICIT NONE
267  DOUBLE PRECISION, INTENT(IN) :: PHI_C
268 
269  IF (phi_c>=zero .AND. phi_c<=one) THEN !monotonic region
270  IF (phi_c >= half) THEN !central differencing
271  minmod = half
272  ELSE !second order upwinding
273  minmod = half*phi_c/(one - phi_c)
274  ENDIF
275  ELSE !first order upwinding
276  minmod = zero
277  ENDIF
278 
279  RETURN
280  END FUNCTION minmod
281 
282 
283 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
284 ! !
285 ! Central scheme for MMS cases. !
286 ! !
287 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
288  DOUBLE PRECISION FUNCTION central_scheme (PHI_C)
289  USE param1, only: half
290  IMPLICIT NONE
291  DOUBLE PRECISION, INTENT(IN) :: PHI_C
292 
293  central_scheme = half
294  RETURN
295  END FUNCTION central_scheme
296 
297 
298 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
299 ! !
300 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
301  DOUBLE PRECISION FUNCTION umist (PHI_C)
302  USE param1, only: zero, half, one
303  IMPLICIT NONE
304  DOUBLE PRECISION, INTENT(IN) :: PHI_C
305 
306  DOUBLE PRECISION :: TH
307 
308  IF (phi_c>=zero .AND. phi_c<one) THEN !monotonic region
309  th = phi_c/(one - phi_c)
310  umist = half*max(zero,min(2.0d0*th,0.75d0 + 0.25d0*th,2.0d0))
311  ELSE IF (phi_c == one) THEN
312  umist = one
313  ELSE !first order upwinding
314  umist = zero
315  ENDIF
316  RETURN
317  END FUNCTION umist
318 
319 
320 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
321 ! !
322 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
323  DOUBLE PRECISION FUNCTION xsi (V, DWF)
324  USE param1, only: zero, one
325  IMPLICIT NONE
326  DOUBLE PRECISION, INTENT(IN) :: V
327  DOUBLE PRECISION, INTENT(IN) :: DWF
328 
329  IF (v >= zero) THEN
330  xsi = dwf
331  ELSE
332  xsi = one - dwf
333  ENDIF
334  RETURN
335  END FUNCTION xsi
336 
337 
338 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
339 ! !
340 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
341  DOUBLE PRECISION FUNCTION phi_c_of (PHI_U, PHI_C, PHI_D)
342  USE param1, only: zero
343  USE toleranc, only: compare
344  IMPLICIT NONE
345  DOUBLE PRECISION, INTENT(IN) :: PHI_U
346  DOUBLE PRECISION, INTENT(IN) :: PHI_C
347  DOUBLE PRECISION, INTENT(IN) :: PHI_D
348 
349  DOUBLE PRECISION :: DEN
350 
351  IF (compare(phi_d,phi_u)) THEN
352  phi_c_of = zero
353  ELSE
354  den = phi_d - phi_u
355  phi_c_of = (phi_c - phi_u)/den
356  ENDIF
357  RETURN
358  END FUNCTION phi_c_of
359 
360 
361 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
362 ! !
363 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
364  DOUBLE PRECISION FUNCTION fpfoi_of (PHI_D, PHI_C, &
365  phi_u, phi_uu)
366  IMPLICIT NONE
367  DOUBLE PRECISION, INTENT(IN) :: PHI_D
368  DOUBLE PRECISION, INTENT(IN) :: PHI_C
369  DOUBLE PRECISION, INTENT(IN) :: PHI_U
370  DOUBLE PRECISION, INTENT(IN) :: PHI_UU
371 
372  fpfoi_of = phi_c + (5.0d0/16.0d0)*(phi_d - phi_c) + &
373  (1.0d0/4.0d0)*(phi_c - phi_u) - &
374  (1.0d0/16.0d0)*(phi_u - phi_uu)
375 
376 ! LIMIT THE HIGH ORDER INTERPOLATION
377  fpfoi_of = univ_limiter_of(fpfoi_of, phi_d, phi_c, phi_u)
378  RETURN
379  END FUNCTION fpfoi_of
380 
381 
382 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
383 ! !
384 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
385  DOUBLE PRECISION FUNCTION univ_limiter_of (PHI_TEMP, PHI_D,&
386  phi_c, phi_u)
387 
388  USE param1, only: zero
389  USE run, only: c_fac
390  IMPLICIT NONE
391  DOUBLE PRECISION, INTENT(IN) :: PHI_TEMP
392  DOUBLE PRECISION, INTENT(IN) :: PHI_D
393  DOUBLE PRECISION, INTENT(IN) :: PHI_C
394  DOUBLE PRECISION, INTENT(IN) :: PHI_U
395 
396  DOUBLE PRECISION :: PHI_REF, DEL, CURV
397 
398  del = phi_d - phi_u
399  curv = phi_d - 2.0*phi_c + phi_u
400  IF (abs(curv) >= abs(del)) THEN ! NON-MONOTONIC
401  univ_limiter_of = phi_c
402  ELSE
403  phi_ref = phi_u + (phi_c-phi_u)/c_fac
404  IF (del > zero) THEN
405  IF (phi_temp < phi_c)univ_limiter_of = phi_c
406  IF (phi_temp > min(phi_ref,phi_d)) &
407  univ_limiter_of = min(phi_ref,phi_d)
408  IF (phi_temp >= phi_c .AND. phi_temp <= min(phi_ref,phi_d)) &
409  univ_limiter_of = phi_temp
410  ELSE
411  IF (phi_temp > phi_c)univ_limiter_of = phi_c
412  IF (phi_temp < max(phi_ref,phi_d)) &
413  univ_limiter_of = max(phi_ref,phi_d)
414  IF (phi_temp >= max(phi_ref,phi_d) .AND. phi_temp <= phi_c) &
415  univ_limiter_of = phi_temp
416  ENDIF
417  ENDIF
418  RETURN
419  END FUNCTION univ_limiter_of
420 
421  END MODULE discretization
toleranc::compare
logical function compare(V1, V2)
Definition: toleranc_mod.f:94
param1::one
double precision, parameter one
Definition: param1_mod.f:29
discretization::chi_smart
double precision function chi_smart(PHI_C, Chi)
Definition: discretization_mod.f:66
discretization::superbee
double precision function superbee(PHI_C)
Definition: discretization_mod.f:20
discretization::quickest
double precision function quickest(PHI_C, CF, ODXC, ODXUC, ODXCD)
Definition: discretization_mod.f:147
xsi
Definition: xsi_mod.f:15
discretization::muscl
double precision function muscl(PHI_C)
Definition: discretization_mod.f:174
discretization::univ_limiter_of
double precision function univ_limiter_of(PHI_TEMP, PHI_D, PHI_C, PHI_U)
Definition: discretization_mod.f:387
discretization::phi_c_of
double precision function phi_c_of(PHI_U, PHI_C, PHI_D)
Definition: discretization_mod.f:342
discretization::fpfoi_of
double precision function fpfoi_of(PHI_D, PHI_C, PHI_U, PHI_UU)
Definition: discretization_mod.f:366
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
discretization::ultra_quick
double precision function ultra_quick(PHI_C, CF)
Definition: discretization_mod.f:116
discretization::umist
double precision function umist(PHI_C)
Definition: discretization_mod.f:302
discretization::minmod
double precision function minmod(PHI_C)
Definition: discretization_mod.f:265
discretization::vanleer
double precision function vanleer(PHI_C)
Definition: discretization_mod.f:248
discretization::central_scheme
double precision function central_scheme(PHI_C)
Definition: discretization_mod.f:289
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
param1::large_number
double precision, parameter large_number
Definition: param1_mod.f:23
discretization::chi4smart
double precision function chi4smart(PHI_C, PHIU, PHIC, PHID)
Definition: discretization_mod.f:88
discretization::chi4muscl
double precision function chi4muscl(PHI_C, PHIU, PHIC, PHID)
Definition: discretization_mod.f:220
discretization::smart
double precision function smart(PHI_C)
Definition: discretization_mod.f:42
discretization::chi_muscl
double precision function chi_muscl(PHI_C, Chi)
Definition: discretization_mod.f:198
run::c_fac
double precision c_fac
Definition: run_mod.f:237
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27