File: N:\mfix\model\chem\stiff_chem_mod.f

1           MODULE STIFF_CHEM
2     
3     ! External Routines.
4     !---------------------------------------------------------------------//
5     ! Routine used to calculate the reaction rates and populate the
6     ! fluid variable ODEs for ODEPACK.
7           external STIFF_CHEM_RRATES
8     ! Routine used to compare to values.
9           LOGICAL, external :: COMPARE
10     ! Routine used to calculate species enthalpies.
11           DOUBLE PRECISION, external :: CALC_H0
12     
13     
14     ! Runtime Flags:
15     !---------------------------------------------------------------------//
16     ! Flag to invoke stiff chemistry solver.
17           LOGICAL :: STIFF_CHEMISTRY
18     ! Flag to invoke variable solids density.
19     !      LOGICAL :: VARIABLE_DENSITY
20     ! Flag indicating if cell IJK is own by myPE.
21           LOGICAL, dimension(:), allocatable :: notOwner
22     
23     ! ODEPACK Controlling parameters:
24     !---------------------------------------------------------------------//
25     ! Dimension of ODEs solved in stiff solver.
26           INTEGER :: ODE_DIMN_all
27     ! Dimension of ODEs solved in stiff solver for gas phase only.
28           INTEGER :: ODE_DIMN_g
29     
30     ! Dimension of ODEs solved in stiff solver for gas phase only.
31           INTEGER :: NEQ_DIMN
32     
33     ! Indicates type of Error control.
34           INTEGER :: ODE_ITOL
35     ! Relative error tolerance paramter.
36           DOUBLE PRECISION, DIMENSION(1) :: ODE_RTOL
37     ! Absolute error tolerance parameter. (Dimension (ODE_DIMN))
38           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: ODE_ATOL
39     ! Declared length of RWORK.
40           INTEGER :: ODE_LRW
41     ! Declared length of IWORK.
42           INTEGER :: ODE_LIW
43     ! Jacobian type indicator.
44           INTEGER :: ODE_JT
45     ! The maximum number of steps ODEPACK may use to integrate.
46           INTEGER :: STIFF_CHEM_MAX_STEPS
47     ! Flag indicating that the max number of steps is unlimited.
48           LOGICAL :: UNLIMITED_STEPS
49     
50     ! Explicit interface for ODEPACK
51     !---------------------------------------------------------------------//
52           INTERFACE
53              SUBROUTINE DLSODA (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, &
54                 ITASK,ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, JT)
55                 external F
56                 INTEGER :: ITOL, ITASK, ISTATE, IOPT, LRW, LIW, JT
57                 INTEGER, dimension(2) :: NEQ
58                 INTEGER, dimension(LIW) :: IWORK
59                 DOUBLE PRECISION :: T, TOUT
60                 DOUBLE PRECISION :: JAC
61                 DOUBLE PRECISION, dimension(1) :: RTOL
62                 DOUBLE PRECISION, dimension(LRW) :: RWORK
63                 DOUBLE PRECISION, dimension(NEQ(1)) :: Y, ATOL
64              END SUBROUTINE DLSODA
65           END INTERFACE
66     
67           contains
68     
69     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
70     !                                                                         C
71     !     Module name: MCHEM_TIME_MARCH                                       C
72     !     Purpose: Called in time_march.f to do rxns calcs                    C
73     !                                                                         C
74     !     Author: Nan Xie                                   Date: 02-Aug-04   C
75     !     Reviewer:                                         Date:             C
76     !                                                                         C
77     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
78           SUBROUTINE STIFF_CHEM_SOLVER(ODE_DT, iErr)
79     
80     ! External Module Procedures:
81     !---------------------------------------------------------------------//
82           use stiff_chem_maps, only : mapODEtoMFIX
83           use stiff_chem_maps, only : mapMFIXtoODE
84     
85     ! Global Variables:
86     !---------------------------------------------------------------------//
87           use output,   only : FULL_LOG
88           use param1,   only : zero
89           use run,      only : TIME
90     
91           use mpi_utility
92     
93           use stiff_chem_dbg
94           use stiff_chem_stats
95           use functions
96     
97           implicit none
98     
99     ! Passed Variables:
100     !----------------------------------------------------------------------!
101     ! Time integral length.
102           DOUBLE PRECISION, intent(IN) :: ODE_DT
103     ! Error Flag
104           INTEGER, intent(OUT) :: iErr
105     
106     ! Local Variables:
107     !----------------------------------------------------------------------!
108     ! Error flag -> Integration failed in one or more fluid cells.
109           LOGICAL :: lErr_l   ! local
110     
111     ! Fluid Cell Index
112           INTEGER :: IJK
113     ! The maximum number of ODEs to solve.
114           DOUBLE PRECISION, dimension(ODE_DIMN_all) :: ODE_VARS
115     
116     ! (1) :: Number of ODEs
117     ! (2) :: Fluid cell index (IJK) passed into ODEPACK
118     ! (:) :: Flag for solving solid phase M density and species equations.
119           INTEGER, dimension(NEQ_DIMN) :: lNEQ
120     ! Start time for integration
121           DOUBLE PRECISION :: lT
122     ! Stop time for integration
123           DOUBLE PRECISION :: lTOUT
124     ! Indicates type of Error control.
125           INTEGER :: lITOL
126     ! Relative error tolerance paramter.
127           DOUBLE PRECISION :: lRTOL(1)
128     ! Absolute error tolerance parameter. (Dimension (ODE_DIMN))
129           DOUBLE PRECISION :: lATOL(ODE_DIMN_all)
130     ! Index specifying the ODEPACK task.
131           INTEGER :: lITASK
132     ! Specifies the state of ODEPACK
133           INTEGER :: lISTATE
134     ! Flag indicating optional inputs are used.
135           INTEGER :: lIOPT
136     ! Array for REAL* work
137           DOUBLE PRECISION :: RWORK(ODE_LRW)
138     ! Declared length of RWORK.
139           INTEGER :: lLRW
140     ! Array for Integer work
141           INTEGER :: IWORK(ODE_LIW)
142     ! Declared length of IWORK.
143           INTEGER :: lLIW
144     ! Jacobain Routine (not used)
145           DOUBLE PRECISION :: lJAC
146     ! Jacobian type indicator.
147           INTEGER :: lJT
148     
149     ! The number of attempts of a specific fluid cell.
150           INTEGER :: lAtps
151     
152           LOGICAL :: lReset
153           LOGICAL :: lIncpt
154     
155           lErr_l = .FALSE.
156     
157           CALL INIT_STIFF_CHEM_STATS
158     
159           IJK_LP: DO IJK = IJKSTART3, IJKEND3
160              IF(notOwner(IJK)) cycle IJK_LP
161              IF(FLUID_AT(IJK)) THEN
162     
163                 lAtps = 0
164                 lReset = .FALSE.
165                 lIncpt = .FALSE.
166     
167     ! Forced restset of tolerance values.
168                 lRTOL = ODE_RTOL
169                 lATOL = ODE_ATOL
170     
171     ! Increment the attempts counter.
172       50        lAtps = lAtps + 1
173     
174     ! Forced restset to initial values.
175                 lT    = 0.0d0
176                 lTOUT = ODE_DT
177                 lITOL = ODE_ITOL
178                 lLRW  = ODE_LRW
179                 lLIW  = ODE_LIW
180                 lJT   = ODE_JT
181     
182     ! Fixed parameters
183                 lITASK  = 1
184                 lISTATE = 1
185                 lIOPT   = 1
186     
187     ! Calculate the number of ODEs to solve.
188                 CALL CALC_ODE_COEFF(lNEQ, IJK)
189     
190     ! Clear the work arrays.
191                 IWORK = 0
192                 RWORK = ZERO
193     
194     ! The maximum number of internal steps ODEPACK may use to integrate over
195     ! the time interval. The default value is 500.
196                 IWORK(6) = STIFF_CHEM_MAX_STEPS
197     
198                 IF(CALC_REACTIONS(IJK)) THEN
199     
200     ! Map MFIX variables to ODE variables.
201                    CALL mapMFIXtoODE(NEQ_DIMN, lNEQ, ODE_DIMN_all, ODE_VARS)
202     
203     ! Store a copy of the original field variables. This allows for these
204     ! values to be 'reset' in the event that the stiff solver fails.
205                    CALL UPDATE_ODE_OLD(ODE_DIMN_all, ODE_VARS)
206     
207     ! Clear the error flag.
208      100           iErr = 0
209     
210     ! Integrate flow field variables to incorporate reactions.
211                    CALL DLSODA(STIFF_CHEM_RRATES, lNEQ, ODE_VARS, lT,      &
212                       lTOUT, lITOL, lRTOL, lATOL, lITASK, lISTATE, lIOPT,  &
213                       RWORK, lLRW, IWORK, lLIW, lJAC, lJT)
214     
215     ! Verify that the results are well defined.
216                    CALL CHECK_ODE_DATA(NEQ_DIMN, lNEQ, ODE_DIMN_all,       &
217                       ODE_VARS, UNLIMITED_STEPS, lISTATE, iErr)
218     
219     ! Successfully Integrated ODEs.
220                    IF(iErr == 0) THEN
221                       lReset = .FALSE.
222     ! Additional integration steps are needed (lT < lTOUT).
223                    ELSEIF(iErr == -1) THEN
224     ! Reste the state flag and keep integrating.
225                       IF(UNLIMITED_STEPS) THEN
226                          lISTATE = 2
227                          goto 100
228                       ELSE
229                          lReset = .FALSE.
230                          lIncpt = .TRUE.
231                       ENDIF
232     
233     ! Too much accuracy was requested.
234                    ELSEIF(iErr == -2) THEN
235                       IF(lAtps < 3) THEN
236     ! Write to the error log file.
237                          IF(ODE_DEBUG_LEVEL >= 2) CALL WRITE_ODE_LOG(iErr, &
238                             NEQ_DIMN, lNEQ, ODE_DIMN_all, ODE_VARS)
239     ! Reset the ODE variable array.
240                          CALL RESET_ODE(ODE_DIMN_all, ODE_VARS, lAtps)
241     ! Reset the field variables.
242                          CALL mapODEtoMFIX(NEQ_DIMN, lNEQ,               &
243                             ODE_DIMN_all, ODE_VARS)
244     ! Loosen the convergence criteria and try again.
245                          lRTOL = ODE_RTOL*10.0d0
246                          lATOL = ODE_ATOL*10.0d0
247                          goto 50
248                       ELSE
249     ! Write to the error log file.
250                          IF(ODE_DEBUG_LEVEL >= 1) CALL WRITE_ODE_LOG(iErr, &
251                             NEQ_DIMN, lNEQ, ODE_DIMN_all, ODE_VARS)
252     ! Set the flag to reset the field variables to the initial values.
253                          lReset = .TRUE.
254                       ENDIF
255     ! All other errors.
256                    ELSE
257     ! Tighten the convergence criteria and try again.
258                       IF(lAtps < 3) THEN
259     ! Write to the error log file.
260                          IF(ODE_DEBUG_LEVEL >= 2) CALL WRITE_ODE_LOG(iErr, &
261                             NEQ_DIMN, lNEQ, ODE_DIMN_all, ODE_VARS)
262     ! Reset the ODE variable array.
263                          CALL RESET_ODE(ODE_DIMN_all, ODE_VARS, lAtps)
264     ! Rest the filed variables to their original values.
265                          CALL mapODEtoMFIX(NEQ_DIMN, lNEQ,                &
266                             ODE_DIMN_all, ODE_VARS)
267     ! Reduce the tolerances and try again.
268                          lRTOL = ODE_RTOL/(10.0d0**lAtps)
269                          lATOL = ODE_ATOL/(10.0d0**lAtps)
270                          goto 50
271                       ELSE
272     ! Write to the error log file.
273                          IF(ODE_DEBUG_LEVEL >= 1) CALL WRITE_ODE_LOG(iErr, &
274                             NEQ_DIMN, lNEQ, ODE_DIMN_all, ODE_VARS)
275     ! Set the flag to reset the field variables to the initial values.
276                          lReset = .TRUE.
277                       ENDIF
278     
279                    ENDIF ! IF(iErr == 0)
280     
281     ! Reset the field variables.
282                    if(lReset) CALL RESET_ODE(ODE_DIMN_all, ODE_VARS, lAtps)
283     ! Store the results in the field variables.
284                    CALL mapODEtoMFIX(NEQ_DIMN, lNEQ, ODE_DIMN_all, ODE_VARS)
285     ! Collect solver stats.
286                    IF(FULL_LOG) CALL UPDATE_STIFF_CHEM_STATS(lNEQ, &
287                       NEQ_DIMN, IWORK(11), ODE_DIMN_all, lAtps, lIncpt)
288     
289     
290                 ENDIF  ! EndIF CALC_REACTIONS
291              ENDIF  ! IF(CALC_REACTIONS(IJK))
292           END DO IJK_LP ! End Loop over fluiod Cells, IJK
293     
294     !      gErr_l = .FALSE.
295     !      CALL GLOBAL_ALL_OR(lErr_l, gErr_l)
296     !      IF(gErr_l) CALL WRITE_VTU_FILE
297     
298     
299           CALL FINALIZE_STIFF_SOLVER()
300     
301           IF(FULL_LOG) CALL WRITE_STIFF_CHEM_STATS()
302     
303           iErr = 0
304     
305           RETURN
306           END SUBROUTINE STIFF_CHEM_SOLVER
307     
308     
309     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
310     !                                                                      !
311     !  Module name: CALC_REACTIONS                                         !
312     !                                                                      !
313     !  Purpose:                                                            !
314     !                                                                      !
315     !  Author: J.Musser                                   Date: 07-Feb-13  !
316     !                                                                      !
317     !  Comments:                                                           !
318     !                                                                      !
319     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
320           LOGICAL FUNCTION CALC_REACTIONS(IJK)
321     
322           use rxns,   only : NO_OF_RXNS
323     
324           implicit none
325     
326           INTEGER, intent(in) :: IJK
327     
328           DOUBLE PRECISION :: RATES(NO_OF_RXNS)
329     
330           DOUBLE PRECISION, parameter :: rLimit = 1.0d-8
331     
332     ! Initialize
333           RATES = 0.0d0
334     
335     ! Calculate user defined reaction rates.
336           CALL USR_RATES(IJK, RATES)
337     
338     ! If there is little to no reaction in the cell, then set the ODE
339     ! Time to zero to avoid calling the stiff solver.
340           CALC_REACTIONS = .TRUE.
341           if(maxval(RATES) < rLimit) CALC_REACTIONS = .FALSE.
342     
343     
344           RETURN
345           END FUNCTION CALC_REACTIONS
346     
347     
348     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
349     !                                                                      !
350     !  Module name: CALC_DIMN_ODE                                          !
351     !                                                                      !
352     !  Purpose:                                                            !
353     !                                                                      !
354     !  Author: J.Musser                                   Date: 07-Feb-13  !
355     !                                                                      !
356     !  Comments:                                                           !
357     !                                                                      !
358     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
359           SUBROUTINE CALC_ODE_COEFF(lNEQ, IJK)
360     
361           use fldvar, only : EP_S
362           use physprop, only : MMAX
363           use run, only : SPECIES_EQ
364     
365           implicit none
366     
367           INTEGER, intent(in)  :: IJK
368           INTEGER, intent(out) :: lNEQ(NEQ_DIMN)
369     
370           INTEGER :: M
371     
372           LOGICAL :: USE_SOLIDS_ODEs
373     
374     ! Initialize.
375           USE_SOLIDS_ODEs = .FALSE.
376           lNEQ(2) = IJK
377           lNEQ(3:) = 0
378     
379     ! If there is little to no solids in the cell, then set the ODE
380     ! dimension to gas phase only.
381           DO M=1, MMAX
382              IF(SPECIES_EQ(M)) THEN
383                 IF(EP_s(IJK,M) > 1.0d-6) THEN
384                    USE_SOLIDS_ODEs = .TRUE.
385                    lNEQ(2+M) = 1
386                 ENDIF
387              ENDIF
388           ENDDO
389     
390           IF(USE_SOLIDS_ODEs)THEN
391              lNEQ(1) = ODE_DIMN_all
392           ELSE
393              lNEQ(1) = ODE_DIMN_g
394              lNEQ(3:) = 0
395           ENDIF
396     
397     
398           RETURN
399           END SUBROUTINE CALC_ODE_COEFF
400     
401     
402     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
403     !                                                                      !
404     !  Module name: FINALIZE_STIFF_SOLVER                                  !
405     !                                                                      !
406     !  Purpose:                                                            !
407     !                                                                      !
408     !  Author: J.Musser                                   Date: 07-Feb-13  !
409     !                                                                      !
410     !  Comments:                                                           !
411     !                                                                      !
412     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
413           SUBROUTINE FINALIZE_STIFF_SOLVER
414     
415     ! Gas phase volume fraction
416           use fldvar, only : EP_g
417     ! Gas phase denisty
418           use fldvar, only: RO_g, ROP_g
419     ! Gas phase pressure
420           use fldvar, only: P_g
421     ! Gas phase temperature
422           use fldvar, only: T_g
423     ! Gas phase species mass fraction
424           use fldvar, only: X_g
425     ! Molecular weight of each gas phase species
426           use physprop, only : MW_g
427     ! Gas phase mixture molecular weight.
428           use physprop, only : MW_MIX_g
429     
430     ! Solids phase bulk density
431           use fldvar, only: ROP_S
432     ! Solids phase temperature
433           use fldvar, only: T_s
434     ! Solids phase species mass fractions
435           use fldvar, only: X_s
436     
437     ! Number of solids phases
438           use physprop, only: MMAX
439     ! Number of species in each phase.
440           use physprop, only: NMAX
441     
442     ! Double precision:  ONE = 1.0d0
443           use param1, only: ONE
444     ! Universal gas constant
445           use constant, only : GAS_CONST
446     
447           use compar
448           use mpi_utility
449           use sendrecv
450           use functions
451           use utilities
452     
453           implicit none
454     
455     ! Local loop indicies.
456           INTEGER :: IJK  ! Fluid Cell index.
457           INTEGER :: M    ! Solids phase index
458           INTEGER :: NN    ! Species index
459     
460           CALL send_recv(EP_G,2)
461           CALL send_recv(RO_G,2)
462           CALL send_recv(ROP_G,2)
463           CALL send_recv(T_G,2)
464     
465           DO NN=1,NMAX(0)
466              CALL send_recv(X_G(:,NN),2)
467              CALL BOUND_X (X_G(1,NN), IJKMAX2)
468           ENDDO
469     
470           DO M = 1, MMAX
471     ! Solids temperature.
472              CALL send_recv(T_S(:,M),2)
473     ! Solids volume fraction. (Constant Solids Density)
474              CALL send_recv(ROP_S(:,M),2)
475     ! Solids phase species mass fractions.
476              DO NN=1,NMAX(M)
477                 CALL send_recv(X_S(:,M,NN),2)
478                 CALL BOUND_X (X_S(1,M,NN), IJKMAX2)
479              ENDDO
480           ENDDO
481     
482           DO IJK = ijkStart3, ijkEnd3
483              IF(.NOT.FLUID_AT(IJK)) CYCLE
484     ! Calculate the mixture molecular weight.
485              MW_MIX_G(IJK) = sum(X_G(IJK,1:NMAX(0))/MW_g(1:NMAX(0)))
486              MW_MIX_G(IJK) = ONE/MW_MIX_G(IJK)
487     ! Calculate the gas phase pressure.
488              P_G(IJK) = (RO_G(IJK)*GAS_CONST*T_G(IJK))/MW_MIX_G(IJK)
489           ENDDO
490     
491           RETURN
492           END SUBROUTINE FINALIZE_STIFF_SOLVER
493     
494           END MODULE STIFF_CHEM
495