File: /nfs/home/0/users/jenkins/mfix.git/model/check_data/check_odepack_stiff_chem.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Module name: CHECK_DATA_CHEM                                        !
4     !  Author: J.Musser                                   Date: 02-Aug-13  !
5     !                                                                      !
6     !  Purpose: Check the chemical rxns namelist variables.                !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE CHECK_ODEPACK_STIFF_CHEM
10     
11     ! Global Variables:
12     !---------------------------------------------------------------------//
13     ! File unit of log file.
14           use funits, only: UNIT_LOG
15     ! Dimension of IJK arrays.
16           use param, only: DIMENSION_3
17     ! Double precision zero.
18           use param1, only: ZERO
19     ! Double precision value for undefined variables.
20           use param1, only: UNDEFINED
21     ! Constant gas phase density
22           use physprop, only : RO_G0
23     ! Runtime logical for solving energy equations.
24           use run, only: ENERGY_EQ
25     ! Run time logical for solving species equations.
26           use run, only: SPECIES_EQ
27     ! Net rate of gas phase production/consumption
28           use rxns, only: SUM_R_g
29     ! Net rate of solids phase production/consumption
30           use rxns, only: SUM_R_s
31     ! Run time logical for using stiff chemistry solver
32           use stiff_chem, only: STIFF_CHEMISTRY
33     ! Run time logicals for identifying cells owned by myPE
34           use stiff_chem, only: notOwner
35     
36     ! Full access to the following modules:
37     !---------------------------------------------------------------------//
38           use compar
39           use geometry
40           use indices
41     
42           use error_manager
43           use functions
44     
45           implicit none
46     
47     ! Local Variables:
48     !---------------------------------------------------------------------//
49           INTEGER :: I, J, K, IJK
50     
51           CALL INIT_ERR_MSG('CHECK_ODEPACK_STIFF_CHEM')
52     
53     ! Verify that there is sufficient run complexity to use the stiff solver
54           IF(STIFF_CHEMISTRY) THEN
55     
56     ! Energy equations must be solved.
57              IF(.NOT.ENERGY_EQ) THEN
58                 WRITE(ERR_MSG,1004)'ENERGY_EQ = .FALSE.'
59                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
60              ENDIF
61     
62     ! Must be compressible
63              IF(RO_G0 /= UNDEFINED) THEN
64                 WRITE(ERR_MSG,1004)'RO_G0 /= UNDEFINED'
65                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
66              ENDIF
67     
68              IF(.NOT.SPECIES_EQ(0)) THEN
69                 WRITE(ERR_MSG,1004)'SPECIES_EQ(0) = .FALSE.'
70                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
71              ENDIF
72     
73     ! The stiff chemistry solver only needs to loop over physical cells
74     ! owned by a process (e.g., not ghost cells). To avoid having a
75     ! triple do loop, this array is populated to identify the cells that
76     ! are not owned.
77              ALLOCATE( notOwner(DIMENSION_3) ); notOwner = .TRUE.
78              do k=kstart, kend
79              do j=jstart, jend
80              do i=istart, iend
81                 ijk = funijk(i,j,k)
82                 notOwner(IJK) = .FALSE.
83              enddo
84              enddo
85              enddo
86     
87     ! Initialize ODEPACK operating parameters.
88              CALL ODEPACK_INIT
89     
90     ! Clear the interphase mass transfer terms as the stiff solver
91     ! does no use them.
92              IF(allocated(SUM_R_g)) SUM_R_g = ZERO
93              IF(allocated(SUM_R_s)) SUM_R_s = ZERO
94     
95           ENDIF
96     
97           CALL FINL_ERR_MSG
98     
99     
100     
101      1004 FORMAT('Error 1004: ',                                           &
102              'Invalid parameters for stiff chemistry solver!',//           &
103              ' The following criteria must be satisfied:',/                &
104              '   > Solving the energy equations.',/                        &
105              '   > Compressible gas phase.',/                              &
106              '   > Solving gas phase species equations.',//                &
107              ' >>> Invalid Parameter: ',A,//                               &
108              'Check the user documentation and correct the mfix.dat file.')
109     
110           RETURN
111           END SUBROUTINE CHECK_ODEPACK_STIFF_CHEM
112     
113     
114     
115     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
116     !                                                                         C
117     !     Module name: MCHEM_ODEPACK_INIT                                     C
118     !     Purpose: controlling values for ODEAPCK(reference to ODEPACK manual)C
119     !                                                                         C
120     !     Author: Nan Xie                                   Date: 02-Aug-04   C
121     !     Reviewer:                                         Date:             C
122     !                                                                         C
123     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
124           SUBROUTINE ODEPACK_INIT
125     
126     ! Global Variables:
127     !---------------------------------------------------------------------//
128           use physprop, only: MMAX
129           use physprop, only: NMAX
130     
131           use run, only: SPECIES_EQ
132     
133           use stiff_chem, only: NEQ_DIMN
134           use stiff_chem, only: ODE_DIMN_g
135           use stiff_chem, only: ODE_DIMN_all
136     
137           use stiff_chem, only: ODE_ITOL
138           use stiff_chem, only: ODE_RTOL
139           use stiff_chem, only: ODE_ATOL
140     
141           use stiff_chem, only: ODE_LRW
142           use stiff_chem, only: ODE_LIW
143           use stiff_chem, only: ODE_JT
144     
145           use stiff_chem, only: STIFF_CHEM_MAX_STEPS
146     
147           use stiff_chem_dbg, only: ALLOCATE_STIFF_CHEM_DBG
148           use stiff_chem_stats, only: ALLOCATE_STIFF_CHEM_STATS
149     
150     ! Double precision value for undefined variables.
151           use param1, only: UNDEFINED_I
152     
153           implicit none
154     
155     
156     ! Local Variables:
157     !---------------------------------------------------------------------//
158           INTEGER :: LRN, LRS
159           INTEGER :: M, N_VAR
160     
161           LOGICAL :: LIMIT_MAX_STEPS
162     
163     ! Number of ODEs (maximum)
164           NEQ_DIMN = 2 + MMAX
165     
166     ! Dimension of ODEs solved if only the gas phase is present:
167     ! Gas density, temperature, and species
168           ODE_DIMN_g = 2 + NMAX(0)
169     ! Solids temperature excluding for all phases.
170           ODE_DIMN_g = ODE_DIMN_g + MMAX
171     
172     ! Calculate the total number of ODEs that are solve.
173           ODE_DIMN_all = ODE_DIMN_g
174           DO M=1, MMAX
175     ! Solids bulk density and species.
176              IF(SPECIES_EQ(M)) ODE_DIMN_all = ODE_DIMN_all + (1 + NMAX(M))
177           ENDDO
178     
179     ! Indicates type of Error control.
180           ODE_ITOL = 2 ! :: EWT(i) = RTOL * ABS(Y(i)) * ATOL(i)
181     
182     ! Relative error tolerance paramter.
183           ODE_RTOL(1) = 1.0D-5
184     
185     ! Absolue error tolerance parameter.
186           IF(.NOT.(allocated(ODE_ATOL))) allocate(ODE_ATOL(ODE_DIMN_all))
187           ODE_ATOL(:) = 1.0D-6  ! All Equations
188     
189     ! Declared length of RWORK.
190           LRN = 20 + 16*ODE_DIMN_all
191           LRS = 22 + 9*ODE_DIMN_all + (ODE_DIMN_all**2)
192           ODE_LRW = max(LRN, LRS)
193     
194     ! Declared length of IWORK.
195           ODE_LIW = 20 + ODE_DIMN_all
196     
197     ! Jacobian type indicator.
198           ODE_JT = 2 ! Internally generated.
199     
200     
201           IF(STIFF_CHEM_MAX_STEPS == UNDEFINED_I) THEN
202              STIFF_CHEM_MAX_STEPS = 500000
203              LIMIT_MAX_STEPS = .FALSE.
204           ELSE
205              LIMIT_MAX_STEPS = .TRUE.
206           ENDIF
207     
208     
209           CALL ALLOCATE_STIFF_CHEM_DBG(ODE_DIMN_all)
210           CALL ALLOCATE_STIFF_CHEM_STATS
211     
212           return
213           END SUBROUTINE ODEPACK_INIT
214