discretization

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

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
422