MFIX  2016-1
stiff_chem_dbg_mod.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2 ! !
3 ! Module name: STIFF_CHEM_DEBUG !
4 ! !
5 ! Purpose: !
6 ! !
7 ! Author: J.Musser Date: !
8 ! !
9 ! Comments: !
10 ! !
11 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
13 
14  PRIVATE
15 
16 ! Variable Access:
17 !---------------------------------------------------------------------//
18  PUBLIC :: ode_debug_level
19 
20 ! Subroutine Access:
21 !---------------------------------------------------------------------//
22  PUBLIC :: allocate_stiff_chem_dbg
23 
24  PUBLIC :: update_ode_old
25  PUBLIC :: reset_ode
26 
27  PUBLIC :: check_ode_data
28  PUBLIC :: write_ode_log
29 
30 ! Routine used to compare to values.
31  LOGICAL, external :: compare
32 
33 ! Static variables/parameters.
34 !---------------------------------------------------------------------//
35 ! Debug Level: (Messages)
36 ! 0 - None
37 ! 1 - Limited
38 ! 2 - Aggressive
39  INTEGER :: ode_debug_level = 1
40 
41 ! File unit for ODE Error Log.
42  INTEGER, parameter :: oel_unit = 6589
43 
44 ! Original field variables. Used to reset after failed integration.
45  DOUBLE PRECISION, allocatable :: ode_vars_o(:)
46 
47  LOGICAL :: print_err_header
48 
49  contains
50 
51 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
52 ! !
53 ! Module name: ALLOCATE_STIFF_CHEM_STATS !
54 ! !
55 ! Purpose: !
56 ! !
57 ! Author: J.Musser Date: !
58 ! !
59 ! Comments: !
60 ! !
61 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
62  SUBROUTINE allocate_stiff_chem_dbg(lODE_DIMN)
63 
64  implicit none
65 
66 ! The number ODEs (maximum).
67  INTEGER, intent(in) :: lODE_DIMN
68 
69  allocate( ode_vars_o(lode_dimn) )
70 
71  RETURN
72  END SUBROUTINE allocate_stiff_chem_dbg
73 
74 
75 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
76 ! Module name: ODE_UPDATE_OLD !
77 ! !
78 ! Purpose: !
79 ! !
80 ! Author: J.Musser Date: !
81 ! !
82 ! Comments: !
83 ! !
84 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
85  SUBROUTINE update_ode_old(lODE_DIMN, lODE_VARS)
86 
87  implicit none
88 
89 ! The number ODEs (maximum).
90  INTEGER, intent(in) :: lODE_DIMN
91 
92 ! MFIX field variables mapped into the ODE format.
93  DOUBLE PRECISION, intent(in) :: lODE_VARS(lode_dimn)
94 
95  ode_vars_o = lode_vars
96 
97  RETURN
98  END SUBROUTINE update_ode_old
99 
100 
101 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
102 ! Module name: ODE_RESET !
103 ! !
104 ! Purpose: !
105 ! !
106 ! Author: J.Musser Date: !
107 ! !
108 ! Comments: !
109 ! !
110 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
111  SUBROUTINE reset_ode(lDIMN, lVARS, lAtps)
113  use compar, only: mype
114 
115  use stiff_chem_stats, only: failedcount
116 
117  implicit none
118 
119 ! The number ODEs (maximum).
120  INTEGER, intent(in) :: lDIMN
121 
122 ! MFIX field variables mapped into the ODE format.
123  DOUBLE PRECISION, intent(inout) :: lVARS(ldimn)
124 
125  INTEGER, intent(in) :: lAtps
126 
127  lvars = ode_vars_o
128 
129  IF(latps == 3) failedcount(mype) = failedcount(mype) + 1
130 
131  RETURN
132  END SUBROUTINE reset_ode
133 
134 
135 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
136 ! !
137 ! Module name: CHECK_ODE_DATA !
138 ! !
139 ! Purpose: !
140 ! !
141 ! Author: J.Musser Date: !
142 ! !
143 ! Comments: !
144 ! !
145 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
146  SUBROUTINE check_ode_data(lnD, lNEQ, loD, VARS, lUNLIMITED, &
147  lstate, lerr)
149  use param1, only : zero
150  use physprop, only : nmax, mmax
151  use toleranc, only : tmax, tmin
152 
153  implicit none
154 
155  INTEGER, intent(in) :: lnD
156 
157  INTEGER, intent(in) :: lNEQ(lnd)
158 
159 ! The number ODEs (maximum).
160  INTEGER, intent(in) :: loD
161 
162 ! MFIX field variables mapped into the ODE format.
163  DOUBLE PRECISION, intent(in) :: VARS(lod)
164 ! Flag for unlimited steps.
165  LOGICAL, intent(in) :: lUNLIMITED
166 ! Sate value returned from ODEPACK
167  INTEGER, intent(in) :: lState
168 ! Error Flag
169  INTEGER, intent(out) :: lErr
170 
171 
172 ! Loop indicies:
173  INTEGER :: M ! phase
174  INTEGER :: N ! species
175  INTEGER :: Node ! other
176 
177 
178  lerr = 0
179 
180 ! Check ODE_VAR for NaNs.
181 !---------------------------------------------------------------------//
182  DO node=1,lod
183  IF(vars(node).NE.vars(node)) THEN
184  lerr = -8
185  return
186  ENDIF
187  ENDDO
188 
189 
190 ! Check ODEPACK State
191 !---------------------------------------------------------------------//
192  IF(lstate == -1) THEN
193  lerr = lstate
194  if(lunlimited) return
195  ELSEIF(lstate == 2) THEN
196  lerr = 0
197  ELSE
198  lerr = lstate
199  ENDIF
200 
201 
202 ! Initialize.
203  node = 1
204 
205 ! Check variables for physical values.
206 !---------------------------------------------------------------------//
207 
208 ! Gas phase density.
209  IF(vars(node) <= zero) THEN
210  lerr = -(100 + node)
211  return
212  ENDIF
213  node = node + 1
214 
215 
216 ! Gas phase temperature.
217  IF(vars(node) > tmax .OR. vars(node) < tmin) THEN
218  lerr = -(100 + node)
219  return
220  ENDIF
221  node = node + 1
222 
223 ! Gas phase species mass fractions.
224  DO n=1,nmax(0)
225  IF(vars(node) > 1.009d0) THEN
226  lerr = -(100 + node)
227  return
228  ENDIF
229  IF(vars(node) < 0.0d0) THEN
230  IF(abs(vars(node)) > 1.0d-5) THEN
231  lerr = -(100 + node)
232  return
233  ENDIF
234  ENDIF
235  node = node + 1
236  ENDDO
237 
238 ! Solids temperature.
239  DO m = 1, mmax
240  IF(vars(node) > tmax .OR. vars(node) < tmin) THEN
241  lerr = -(100 + node)
242  return
243  ENDIF
244  node = node + 1
245  ENDDO
246 
247 
248  DO m = 1, mmax
249  IF(lneq(2+m) == 1) THEN
250 
251 ! Solids volume fraction.
252  IF(vars(node) <= zero) THEN
253  lerr = -(100 + node)
254  return
255  ENDIF
256  node = node + 1
257 
258 ! Solids phase species mass fractions.
259  DO n=1,nmax(m)
260 
261  IF(vars(node) > 1.009d0) THEN
262  lerr = -(100 + node)
263  return
264  ENDIF
265  IF(vars(node) < 0.0d0) THEN
266  IF(abs(vars(node)) > 1.0d-5) THEN
267  lerr = -(100 + node)
268  return
269  ENDIF
270  ENDIF
271  node = node + 1
272 
273  ENDDO
274  ENDIF
275  ENDDO
276 
277  RETURN
278  END SUBROUTINE check_ode_data
279 
280 
281 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
282 ! !
283 ! Module name: ODE_ErrorLog !
284 ! !
285 ! Purpose: !
286 ! !
287 ! Author: J.Musser Date: !
288 ! !
289 ! Comments: !
290 ! !
291 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
292  SUBROUTINE write_ode_log(lErr, lnD, lNEQ, loD, lVARS)
294  use compar, only : mype
295  use run, only : time
296  use physprop, only : mmax
297 
298  use indices
299 
300  implicit none
301 
302  INTEGER, intent(in) :: lErr
303 
304  INTEGER, intent(in) :: lnD
305  INTEGER, intent(in) :: lNEQ(lnd)
306 
307 ! The number ODEs (maximum).
308  INTEGER, intent(in) :: loD
309  DOUBLE PRECISION, intent(in) :: lVARS(lod)
310 
311  INTEGER :: lc1
312  INTEGER :: IJK
313 
314  DOUBLE PRECISION :: ddt_lVARS(lod)
315 
316  CHARACTER(LEN=255) :: lFile
317 
318  LOGICAL :: lExist
319 
320  ijk = lneq(2)
321 
322  lfile = ''; write(lfile,"('StiffChem_',I2.2,'.log')") mype
323 
324  inquire(file=lfile,exist=lexist)
325  if(lexist) then
326  open(oel_unit,file=trim(adjustl(lfile)), position='append',convert='big_endian')
327  else
328  open(oel_unit,file=trim(adjustl(lfile)), status='new',convert='big_endian')
329  endif
330 
331 
332  if(print_err_header) then
333  WRITE(oel_unit,9000) time
334  print_err_header = .false.
335  endif
336 
337 
338  WRITE(oel_unit,9001) ijk, i_of(ijk), j_of(ijk), mype
339  WRITE(oel_unit,9002) lerr
340  WRITE(oel_unit,9003) lneq(1)
341 
342  DO lc1=1,mmax
343  IF(lneq(lc1+2) == 1) THEN
344  WRITE(oel_unit,9004) lc1
345  ELSEIF(lneq(lc1+2) == 0) THEN
346  WRITE(oel_unit,9005) lc1
347  ELSE
348  WRITE(oel_unit,9006) lc1, lneq(lc1+2)
349  ENDIF
350  ENDDO
351 
352  WRITE(oel_unit,"(/5x,'Field Variables:')")
353  WRITE(oel_unit,9007)
354  DO lc1=1, lod
355  WRITE(oel_unit,9008) lc1, ode_vars_o(lc1), lvars(lc1), &
356  (lvars(lc1) - ode_vars_o(lc1))
357  ENDDO
358 
359  CALL stiff_chem_rrates(lneq, time, lvars, ddt_lvars)
360 
361  WRITE(oel_unit,"(/5x,'Derivatives:')")
362  WRITE(oel_unit,9009)
363  DO lc1=1, lod
364  WRITE(oel_unit,9010) lc1, ddt_lvars(lc1)
365  ENDDO
366 
367 
368  RETURN
369 
370  9000 FORMAT(//3x,'Time: ',g15.8)
371  9001 FORMAT(//5x,'Fluid Cell: ',i4,4x,'(',i3,',',i3,')',10x,&
372  'Process: ',i4)
373  9002 FORMAT( 5x,'Error Flag/ODEPACK Status: ',i4,/)
374 
375  9003 FORMAT( 5x,'Number of ODEs Solved: ', i4,/)
376  9004 FORMAT( 5x,'Solving solids phase ',i2,&
377  ' bulk density and species equations.')
378 
379  9005 FORMAT( 5x,'NOT solving solids phase ',i2,&
380  ' bulk density and species equations.')
381 
382  9006 FORMAT( 5x,'Unknown flag for solids phase ',i2,' Flag = ',i4)
383 
384 
385  9007 FORMAT( 5x,'Var',6x,'Incoming',13x,'Returned',13x,'Difference')
386  9008 FORMAT( 5x,i3,3(3x,g18.8))
387 
388 
389  9009 FORMAT( 5x,'Var',6x,' YDOT ')
390  9010 FORMAT( 5x,i3,3x,g18.8)
391 
392  END SUBROUTINE write_ode_log
393 
394 
395  END MODULE stiff_chem_dbg
subroutine, public write_ode_log(lErr, lnD, lNEQ, loD, lVARS)
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
double precision, parameter tmax
Definition: toleranc_mod.f:31
subroutine, public check_ode_data(lnD, lNEQ, loD, VARS, lUNLIMITED, lState, lErr)
double precision, parameter tmin
Definition: toleranc_mod.f:34
subroutine, public update_ode_old(lODE_DIMN, lODE_VARS)
integer, dimension(:), allocatable, public failedcount
integer, public ode_debug_level
integer mmax
Definition: physprop_mod.f:19
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
Definition: run_mod.f:13
integer, dimension(0:dim_m) nmax
Definition: physprop_mod.f:119
integer mype
Definition: compar_mod.f:24
subroutine stiff_chem_rrates(lNEQ, lTime, Y, YDOT)
subroutine, public reset_ode(lDIMN, lVARS, lAtps)
double precision time
Definition: run_mod.f:45
double precision, parameter zero
Definition: param1_mod.f:27
subroutine, public allocate_stiff_chem_dbg(lODE_DIMN)