module_cbmz_rodas3_solver.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! CBMZ module: see module_cbmz.F for information and terms of use
8 !**********************************************************************************
9 module module_cbmz_rodas3_solver
10
11 contains
12
13
14 !-----------------------------------------------------------------------
15
16 subroutine rodas3_ff_x2( n, t, tnext, hmin, hmax, hstart, &
17 y, abstol, reltol, yposlimit, yneglimit, &
18 yfixed, rconst, &
19 lu_nonzero_v, lu_crow_v, lu_diag_v, lu_icol_v, &
20 info, iok, lunerr, &
21 dydtsubr, yjacsubr, decompsubr, solvesubr )
22 ! include 'sparse_s.h'
23
24 ! stiffly accurate rosenbrock 3(2), with
25 ! stiffly accurate embedded formula for error control.
26 !
27 ! all the arguments aggree with the kpp syntax.
28 !
29 ! input arguments:
30 ! n = number of dependent variables (i.e., time-varying species)
31 ! [t, tnext] = the integration interval
32 ! hmin, hmax = lower and upper bounds for the selected step-size.
33 ! note that for step = hmin the current computed
34 ! solution is unconditionally accepted by the error
35 ! control mechanism.
36 ! hstart = initial guess step-size
37 ! y = vector of (n) concentrations, contains the
38 ! initial values on input
39 ! abstol, reltol = (n) dimensional vectors of
40 ! componentwise absolute and relative tolerances.
41 ! yposlimit = vector of (n) upper-limit positive concentrations
42 ! yneglimit = vector of (n) upper-limit negative concentrations
43 ! yfixed = vector of (*) concentrations of fixed/non-reacting species
44 ! rconst = vector of (*) concentrations of fixed/non-reacting species
45 !
46 ! lu_nonzero_v = number of non-zero entries in the "augmented jacobian"
47 ! (i.e., the jacobian with all diagonal elements filled in)
48 ! lu_crow_v = vector of (n) pointers to the crow (?) elements
49 ! lu_diag_v = vector of (n) pointers to the diagonal elements
50 ! of the sparsely organized jacobian.
51 ! jac(lu_diag_v(i)) is the i,i diagonal element
52 ! lu_icol_v = vector of (lu_nonzero_v) pointers to the icol (?) elements
53 ! info(1) = 1 for autonomous system
54 ! = 0 for nonautonomous system
55 ! iok = completion code (output)
56 ! +1 = successful integration
57 ! +2 = successful integration, but some substeps were h=hmin
58 ! and error > abstol,reltol
59 ! -1 = failure -- lu decomposition failed with minimum stepsize h=hmin
60 ! -1001 --> -1999 = failure -- with minimum stepsize h=hmin,
61 ! species i = -(iok+1000) was NaN
62 ! -2001 --> -2999 = failure -- with minimum stepsize h=hmin,
63 ! species i = -(iok+2000) was > yneglimit
64 ! -3001 --> -3999 = failure -- with minimum stepsize h=hmin,
65 ! species i = -(iok+3000) was < yposlimit
66 ! -5001 --> -5999 = failure -- with minimum stepsize h=hmin,
67 ! species i = -(iok+5000) had abs(kn(i)) > ylimit_solvesubr
68 ! before call to solvesubr, where kn is either k1,k2,k3,k4
69 !
70 ! dydtsubr = name of routine of derivatives. kpp syntax.
71 ! see the header below.
72 ! yjacsubr = name of routine that computes the jacobian, in
73 ! sparse format. kpp syntax. see the header below.
74 ! decompsubr = name of routine that does sparse lu decomposition
75 ! solvesubr = name of routine that does sparse lu backsolve
76 !
77 ! output arguments:
78 ! y = the values of concentrations at tend.
79 ! t = equals tend on output.
80 ! info(2) = # of dydtsubr calls.
81 ! info(3) = # of yjacsubr calls.
82 ! info(4) = # of accepted steps.
83 ! info(5) = # of rejected steps.
84 ! info(6) = # of steps that were accepted but had
85 ! (hold.eq.hmin) .and. (err.gt.1)
86 ! which means stepsize=minimum and error>tolerance
87 !
88 ! adrian sandu, march 1996
89 ! the center for global and regional environmental research
90
91 use module_peg_util, only: peg_message, peg_error_fatal
92
93 implicit none
94
95 integer nvar_maxd, lu_nonzero_v_maxd
96 parameter (nvar_maxd=99)
97 parameter (lu_nonzero_v_maxd=999)
98
99 ! common block variables
100
101 ! subr parameters
102 integer n, info(6), iok, lunerr, lu_nonzero_v
103 integer lu_crow_v(n), lu_diag_v(n), lu_icol_v(lu_nonzero_v)
104 real t, tnext, hmin, hmax, hstart
105 real y(n), abstol(n), reltol(n), yposlimit(n), yneglimit(n)
106 real yfixed(*), rconst(*)
107 external dydtsubr, yjacsubr, decompsubr, solvesubr
108
109 ! local variables
110 logical isreject, autonom
111 integer nfcn, njac, naccept, nreject, nnocnvg, i, j
112 integer ier
113
114 real k1(nvar_maxd), k2(nvar_maxd)
115 real k3(nvar_maxd), k4(nvar_maxd)
116 real f1(nvar_maxd), ynew(nvar_maxd)
117 real jac(lu_nonzero_v_maxd)
118 real ghinv, uround
119 real tin, tplus, h, hold, hlowest
120 real err, factor, facmax
121 real dround, c43, tau, x1, x2, ytol
122 real ylimit_solvesubr
123
124 character*80 errmsg
125
126 ylimit_solvesubr = 1.0e18
127
128 ! check n and lu_nonzero_v
129 if (n .gt. nvar_maxd) then
130 call peg_message( lunerr, '*** rodas3 dimensioning problem' )
131 write(errmsg,9050) 'n, nvar_maxd = ', n, nvar_maxd
132 call peg_message( lunerr, errmsg )
133 call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
134 else if (lu_nonzero_v .gt. lu_nonzero_v_maxd) then
135 call peg_message( lunerr, '*** rodas3 dimensioning problem' )
136 write(errmsg,9050) 'lu_nonvero_v, lu_nonzero_v_maxd = ', &
137 lu_nonzero_v, lu_nonzero_v_maxd
138 call peg_message( lunerr, errmsg )
139 call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
140 end if
141 9050 format( a, 2(1x,i6) )
142
143 ! initialization
144 uround = 1.e-7
145 dround = sqrt(uround)
146
147 ! check hmin and hmax
148 hlowest = dround
149 if (info(1) .eq. 1) hlowest = 1.0e-7
150 if (hmin .lt. hlowest) then
151 call peg_message( lunerr, '*** rodas3 -- hmin is too small' )
152 write(errmsg,9060) 'hmin and minimum allowed value = ', &
153 hmin, hlowest
154 call peg_message( lunerr, errmsg )
155 call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
156 else if (hmin .ge. hmax) then
157 call peg_message( lunerr, '*** rodas3 -- hmin >= hmax' )
158 write(errmsg,9060) 'hmin, hmax = ', hmin, hmax
159 call peg_message( lunerr, errmsg )
160 call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
161 end if
162 9060 format( a, 1p, 2e14.4 )
163
164 ! initialization of counters, etc.
165 autonom = info(1) .eq. 1
166 !##checkthis##
167 c43 = - 8.e0/3.e0
168 ! h = 60.e0 ! hmin
169 h = max( hstart, hmin )
170 tplus = t
171 tin = t
172 isreject = .false.
173 naccept = 0
174 nreject = 0
175 nnocnvg = 0
176 nfcn = 0
177 njac = 0
178
179
180 ! === starting the time loop ===
181 10 continue
182 tplus = t + h
183 if ( tplus .gt. tnext ) then
184 h = tnext - t
185 tplus = tnext
186 end if
187
188 call yjacsubr( n, t, y, jac, yfixed, rconst )
189 njac = njac+1
190 ghinv = -2.0e0/h
191 do 20 j=1,n
192 jac(lu_diag_v(j)) = jac(lu_diag_v(j)) + ghinv
193 20 continue
194 call decompsubr( n, jac, ier, lu_crow_v, lu_diag_v, lu_icol_v )
195
196 if (ier.ne.0) then
197 if ( h.gt.hmin) then
198 h = 5.0e-1*h
199 go to 10
200 else
201 ! print *,'ier <> 0, h=',h
202 ! stop
203 iok = -1
204 goto 200
205 end if
206 end if
207
208 call dydtsubr( n, t, y, f1, yfixed, rconst )
209
210 ! ====== nonautonomous case ===============
211 if (.not. autonom) then
212 ! tau = sign(dround*max( 1.0e-6, abs(t) ), t)
213 tau = dround*max( 1.0e-6, abs(t) )
214 call dydtsubr( n, t+tau, y, k2, yfixed, rconst )
215 nfcn=nfcn+1
216 do 30 j = 1,n
217 k3(j) = ( k2(j)-f1(j) )/tau
218 30 continue
219
220 ! ----- stage 1 (nonautonomous) -----
221 x1 = 0.5*h
222 do 40 j = 1,n
223 k1(j) = f1(j) + x1*k3(j)
224 if (abs(k1(j)) .gt. ylimit_solvesubr) then
225 iok = -(5000+j)
226 goto 135
227 end if
228 40 continue
229 call solvesubr( jac, k1 )
230
231 ! ----- stage 2 (nonautonomous) -----
232 x1 = 4.e0/h
233 x2 = 1.5e0*h
234 do 50 j = 1,n
235 k2(j) = f1(j) - x1*k1(j) + x2*k3(j)
236 if (abs(k2(j)) .gt. ylimit_solvesubr) then
237 iok = -(5000+j)
238 goto 135
239 end if
240 50 continue
241 call solvesubr( jac, k2 )
242
243 ! ====== autonomous case ===============
244 else
245 ! ----- stage 1 (autonomous) -----
246 do 60 j = 1,n
247 k1(j) = f1(j)
248 if (abs(k1(j)) .gt. ylimit_solvesubr) then
249 iok = -(5000+j)
250 goto 135
251 end if
252 60 continue
253 call solvesubr( jac, k1 )
254
255 ! ----- stage 2 (autonomous) -----
256 x1 = 4.e0/h
257 do 70 j = 1,n
258 k2(j) = f1(j) - x1*k1(j)
259 if (abs(k2(j)) .gt. ylimit_solvesubr) then
260 iok = -(5000+j)
261 goto 135
262 end if
263 70 continue
264 call solvesubr( jac, k2 )
265 end if
266
267 ! ----- stage 3 -----
268 do 80 j = 1,n
269 ynew(j) = y(j) - 2.0e0*k1(j)
270 80 continue
271 call dydtsubr( n, t+h, ynew, f1, yfixed, rconst )
272 nfcn=nfcn+1
273 do 90 j = 1,n
274 k3(j) = f1(j) + ( -k1(j) + k2(j) )/h
275 if (abs(k3(j)) .gt. ylimit_solvesubr) then
276 iok = -(5000+j)
277 goto 135
278 end if
279 90 continue
280 call solvesubr( jac, k3 )
281
282 ! ----- stage 4 -----
283 do 100 j = 1,n
284 ynew(j) = y(j) - 2.0e0*k1(j) - k3(j)
285 100 continue
286 call dydtsubr( n, t+h, ynew, f1, yfixed, rconst )
287 nfcn=nfcn+1
288 do 110 j = 1,n
289 k4(j) = f1(j) + ( -k1(j) + k2(j) - c43*k3(j) )/h
290 if (abs(k4(j)) .gt. ylimit_solvesubr) then
291 iok = -(5000+j)
292 goto 135
293 end if
294 110 continue
295 call solvesubr( jac, k4 )
296
297 ! ---- the solution ---
298
299 do 120 j = 1,n
300 ynew(j) = y(j) - 2.0e0*k1(j) - k3(j) - k4(j)
301 120 continue
302
303
304 ! ====== error estimation ========
305
306 err=0.e0
307 do 130 i=1,n
308 ytol = abstol(i) + reltol(i)*abs(ynew(i))
309 err = err + ( k4(i)/ytol )**2
310 130 continue
311 err = max( uround, sqrt( err/n ) )
312
313 ! ======= choose the stepsize ===============================
314
315 factor = 0.9/err**(1.e0/3.e0)
316 if (isreject) then
317 facmax=1.0
318 else
319 facmax=10.0
320 end if
321 factor = max( 1.0e-1, min(factor,facmax) )
322 hold = h
323 h = min( hmax, max(hmin,factor*h) )
324
325 ! ======= check for nan, too big, too small =================
326 ! if any of these conditions occur, either
327 ! exit if hold <= hmin
328 ! reduce step size and retry if hold > hmin
329
330 iok = 1
331 do i = 1, n
332 if (y(i) .ne. y(i)) then
333 iok = -(1000+i)
334 goto 135
335 else if (y(i) .lt. yneglimit(i)) then
336 iok = -(2000+i)
337 goto 135
338 else if (y(i) .gt. yposlimit(i)) then
339 iok = -(3000+i)
340 goto 135
341 end if
342 end do
343 135 if (iok .lt. 0) then
344 if (hold .le. hmin) then
345 goto 200
346 else
347 isreject = .true.
348 nreject = nreject+1
349 h = max(hmin, 0.5*hold)
350 if (t.eq.tin) h = max(hmin, 1.0e-1*hold)
351 goto 10
352 end if
353 end if
354
355
356 ! ======= rejected/accepted step ============================
357
358 if ( (err.gt.1).and.(hold.gt.hmin) ) then
359 isreject = .true.
360 nreject = nreject+1
361 if (t.eq.tin) h = max(hmin, 1.0e-1*hold)
362 else
363 isreject = .false.
364 do 140 i=1,n
365 y(i) = ynew(i)
366 140 continue
367 t = tplus
368 if (err.gt.1) then
369 nnocnvg = nnocnvg+1
370 else
371 naccept = naccept+1
372 end if
373 end if
374
375 ! ======= end of the time loop ===============================
376 if ( t .lt. tnext ) go to 10
377
378 iok = 1
379 if (nnocnvg .gt. 0) iok = 2
380
381
382 ! ======= output information =================================
383 200 info(2) = nfcn
384 info(3) = njac
385 info(4) = naccept
386 info(5) = nreject
387 info(6) = nnocnvg
388
389 return
390 end subroutine rodas3_ff_x2
391
392
393 end module module_cbmz_rodas3_solver