Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2016-2019 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module net_burn_support
21 : use const_def, only: dp
22 : use math_lib
23 : use utils_lib, only: is_bad, mesa_error
24 : use mtx_def, only: lapack
25 :
26 : implicit none
27 :
28 : integer, parameter :: qcol_imax=13, kmaxx = 7, stifbs_imax = kmaxx+1
29 : logical, parameter :: dbg = .false.
30 :
31 : contains
32 :
33 2 : subroutine netint( &
34 2 : start,stptry,stpmin,max_steps,stopp,y, &
35 2 : eps,species,nvar,nok,nbad,nstp,odescal,dens_dfdy,dmat, &
36 : derivs,jakob,burner_finish_substep,ierr)
37 :
38 : ! input:
39 : ! start = beginning integration point
40 : ! stptry = suggested first step size
41 : ! stpmin = minimum allowable step size
42 : ! stopp = ending integration point
43 : ! bc = initial conditions, array of physical dimension yphys
44 : ! eps = desired fraction error during the integration
45 : ! odescal = error scaling factor
46 : ! derivs = name of the routine that contains the odes
47 : ! jakob = name of the routine that contains the jacobian of the odes
48 : ! steper = name of the routine that will take a single step
49 :
50 : ! nvar >= species. abundances come first in args,
51 : ! optionally followed by anything else such as lnT.
52 :
53 : ! output:
54 : ! nok = number of successful steps taken
55 : ! nbad = number of bad steps taken, bad but retried and then successful
56 : ! nstp = total number of steps taken
57 :
58 : real(dp) :: dens_dfdy(:,:),dmat(:,:)
59 : integer, intent(out) :: ierr
60 :
61 : interface
62 : include 'burner_derivs.inc'
63 : end interface
64 : interface
65 : include 'burner_jakob.inc'
66 : end interface
67 : interface
68 : include 'burner_finish_substep.inc'
69 : end interface
70 :
71 : integer :: species,nvar,nok,nbad,nstp,max_steps
72 : real(dp) :: start,stptry,stpmin,stopp,y(:),eps,odescal
73 :
74 88 : real(dp) :: yscal(nvar),dydx(nvar),cons,x,h,hdid,hnext,xx
75 : real(dp), parameter :: tiny=1.0d-15
76 :
77 176 : real(dp) :: y0(nvar),a(stifbs_imax),alf(kmaxx,kmaxx),epsold,xnew,scale,red
78 : integer :: i,kmax,kopt,nseq(stifbs_imax),nvold
79 : logical :: first
80 :
81 : include 'formats'
82 :
83 : ! initialize
84 2 : first = .true.
85 2 : epsold = -1d0
86 2 : nvold = -1
87 2 : nseq = [ 2, 6, 10, 14, 22, 34, 50, 70 ]
88 2 : x = start
89 2 : h = sign(stptry,stopp-start)
90 2 : nok = 0
91 2 : nbad = 0
92 : ierr = 0
93 :
94 44 : do i=1,nvar
95 44 : y0(i) = y(i)
96 : end do
97 :
98 : ! take at most max_steps steps
99 54 : do nstp=1,max_steps
100 :
101 : ! positive definite abundance fractions
102 1188 : do i=1,species
103 : if (dbg) then
104 : if (is_bad(y(i))) then
105 : write(*,*) 'bad y for nstp', nstp, i, y(i)
106 : call mesa_error(__FILE__,__LINE__,'netint')
107 : end if
108 : end if
109 1188 : y(i) = min(1.0d0, max(y(i),1.0d-30))
110 : end do
111 :
112 54 : call burner_finish_substep(nstp, x, y, ierr)
113 54 : if (ierr /= 0) return
114 :
115 : ! get the right hand sides and form scaling vector
116 54 : call derivs(x,y,dydx,nvar,ierr)
117 54 : if (ierr /= 0) then
118 : return
119 : write(*,*) 'derivs failed in netint'
120 : end if
121 :
122 1188 : do i=1,nvar
123 1188 : yscal(i) = max(odescal,abs(y(i)))
124 : end do
125 :
126 : ! if the step can overshoot the stop point cut it
127 54 : if ((x+h-stopp)*(x+h-start) > 0.0d0) h = stopp - x
128 :
129 : ! do an integration step
130 : call stifbs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext, &
131 : a,alf,epsold,first,kmax,kopt,nseq,nvold,xnew,scale,red, &
132 54 : dens_dfdy,dmat,derivs,jakob,nstp,ierr)
133 54 : if (ierr /= 0) then
134 0 : write(*,*) 'stifbs ierr'
135 0 : return
136 : end if
137 :
138 : ! tally if we did or did not do the step
139 54 : if (hdid==h) then
140 46 : nok = nok+1
141 : else
142 8 : nbad = nbad+1
143 : end if
144 :
145 : ! this is the normal exit point
146 54 : if ( (nstp == max_steps) .or. &
147 : (x-stopp)*(stopp-start)>= 0.0d0) then
148 2 : call burner_finish_substep(nstp, x, y, ierr)
149 2 : return
150 : end if
151 :
152 : ! normal timestep choice
153 52 : h = hnext
154 :
155 : ! die
156 52 : if (abs(h)<stpmin) then
157 0 : write(*,*) 'netint failed: abs(h).lt.stpmin', abs(h), stpmin
158 0 : ierr = -1
159 0 : return
160 :
161 : write(6,210) 'nstp=',nstp
162 : write(6,210) 'nok nbad',nok,nbad
163 : write(6,220) 'attempted time step =',stptry
164 : write(6,220) 'current time step =',h
165 : write(6,220) 'input composition:'
166 : write(6,230) (y0(i), i=1,nvar)
167 : write(6,220) 'current composition:'
168 : write(6,230) (y(i), i=1,nvar)
169 :
170 : !call mesa_error(__FILE__,__LINE__,'h < stpmin in netint')
171 :
172 : 210 format(1x,a,4i6)
173 : 220 format(1x,a,1p3e24.16)
174 : 230 format(1x,1p3e24.16)
175 :
176 : end if
177 :
178 :
179 : ! back for another iteration or death
180 : end do
181 0 : ierr = -1
182 0 : return
183 : call mesa_error(__FILE__,__LINE__,'more than max_steps steps required in netint')
184 : end subroutine netint
185 :
186 :
187 54 : subroutine stifbs(y,dydx,nvar,x,htry,eps,yscal,hdid,hnext, &
188 : a,alf,epsold,first,kmax,kopt,nseq,nvold,xnew,scale,red, &
189 54 : dens_dfdy,dmat,derivs,jakob,nstp,ierr)
190 : !
191 : ! for dense analytic jacobians, lu decomposition linear algebra
192 : !
193 : ! semi-implicit extrapolation step for integrating stiff ode's with monitoring
194 : ! of local truncation error to adjust stepsize. inputs are the dependent
195 : ! variable vector y(1:nvar) and its derivative dydx(1:nvar) at the starting of the
196 : ! independent variable x. also input are the stepsize to be attempted htry,
197 : ! the required accuracy eps, and the vector yscal against which the error is
198 : ! scaled. on output, y and x are replaced by their new values, hdid is the
199 : ! stepsize actually accomplished, and hnext is the estimated next stepsize.
200 : ! dervs is a user supplied function that computes the right hand side of
201 : ! the equations.
202 : !
203 : ! declare
204 : real(dp) :: a(stifbs_imax),alf(kmaxx,kmaxx),epsold,xnew,scale,red
205 : integer :: kmax,kopt,nseq(stifbs_imax),nvold
206 : logical :: first
207 : integer, intent(out) :: ierr
208 :
209 :
210 : real(dp) :: dens_dfdy(:,:),dmat(:,:)
211 :
212 : interface
213 : include 'burner_derivs.inc'
214 : end interface
215 : interface
216 : include 'burner_jakob.inc'
217 : end interface
218 :
219 : logical :: reduct
220 : integer :: nvar
221 : integer :: i,iq,j,k,kk,km,i_errmax
222 : real(dp) :: y(:),dydx(:),x,htry,eps,yscal(:),hdid,hnext, &
223 486 : eps1,errmax,fact,h,work,wrkmin,xest,err(kmaxx), &
224 3456 : yerr(nvar),ysav(nvar),yseq(nvar),safe1,safe2,dum, &
225 16254 : redmax,redmin,tiny,scalmx,qcol(nvar,qcol_imax),x_pzextr(qcol_imax)
226 : parameter (safe1 = 0.25d0, safe2 = 0.7d0, redmax=1.0d-5, &
227 : redmin = 0.7d0, tiny = 1.0d-30, &
228 : scalmx = 0.1d0)
229 :
230 : integer :: nstp, ierr2
231 :
232 54 : ierr = 0
233 :
234 : ! a new tolerance or a new number, so reinitialize
235 54 : if (eps /= epsold .or. nvar /= nvold) then
236 2 : hnext = -1.0d29
237 2 : xnew = -1.0d29
238 2 : eps1 = safe1 * eps
239 :
240 : ! compute the work coefficients a_k
241 2 : a(1) = nseq(1) + 1
242 16 : do k=1,kmaxx
243 16 : a(k+1) = a(k) + nseq(k+1)
244 : end do
245 :
246 : ! compute alf(k,q)
247 14 : do iq=2,kmaxx
248 56 : do k=1,iq-1
249 : alf(k,iq) = pow(eps1,((a(k+1) - a(iq+1)) / &
250 54 : ((a(iq+1) - a(1) + 1.0d0) * (2*k + 1))))
251 : end do
252 : end do
253 2 : epsold = eps
254 2 : nvold = nvar
255 :
256 : ! add cost of jacobians to work coefficients
257 2 : a(1) = nvar + a(1)
258 16 : do k=1,kmaxx
259 16 : a(k+1) = a(k) + nseq(k+1)
260 : end do
261 :
262 : ! determine optimal row number for convergence
263 12 : do kopt=2,kmaxx-1
264 12 : if (a(kopt+1) > a(kopt)*alf(kopt-1,kopt)) GOTO 01
265 : end do
266 2 : 01 kmax = kopt
267 : end if
268 :
269 : ! save the starting values
270 54 : h = htry
271 1188 : do i=1,nvar
272 1188 : ysav(i) = y(i)
273 : end do
274 :
275 : ! get the dense jacobian in dens_dfdy
276 54 : call jakob(x,y,dens_dfdy,nvar,ierr)
277 54 : if (ierr /= 0) then
278 : if (dbg) write(*,*) 'jakob failed in stifbs'
279 : return
280 : end if
281 :
282 : ! a new stepsize or a new integration, re-establish the order window
283 54 : if (h /= hnext .or. x /= xnew) then
284 4 : first = .true.
285 4 : kopt = kmax
286 : end if
287 : reduct = .false.
288 :
289 : ! evaluate the sequence of semi implicit midpoint rules
290 264 : 02 do k=1,kmax
291 :
292 : ! write(6,119) 'xnew x and h',xnew,x,h
293 : ! 119 format(1x,a,' ',1p3e12.4)
294 :
295 264 : xnew = x + h
296 :
297 : ! write(6,119) 'xnew x and h',xnew,x,h
298 : ! read(5,*)
299 :
300 264 : if (xnew == x) then
301 0 : ierr = -1
302 : if (dbg) write(*,*) 'ierr: stepsize too small in routine stiffbs'
303 0 : return
304 : call mesa_error(__FILE__,__LINE__,'stepsize too small in routine stiffbs')
305 : end if
306 :
307 : call simpr( &
308 : ysav,dydx,nvar,x,h,nseq(k),yseq,dens_dfdy,dmat, &
309 264 : derivs,ierr)
310 264 : if (ierr /= 0) then
311 :
312 0 : h = h * 0.1d0
313 0 : i_errmax = 0
314 0 : reduct = .true.
315 0 : ierr = 0
316 : if (dbg) write(*,*) 'ierr: simpr failed in stifbs'
317 0 : GOTO 02
318 :
319 : write(*,*) 'simpr failed in stifbs'
320 : return
321 :
322 : end if
323 264 : xest = (h/nseq(k))*(h/nseq(k))
324 264 : call net_pzextr(k,xest,yseq,y,yerr,nvar,qcol,x_pzextr)
325 :
326 : ! compute normalized error estimate
327 264 : if (k /= 1) then
328 198 : errmax = tiny
329 198 : i_errmax = 0
330 4356 : do i=1,nvar
331 : ! errmax = max(errmax,abs(yerr(i)/yscal(i)))
332 4158 : dum = abs(yerr(i)/yscal(i))
333 4356 : if (dum >= errmax) then
334 732 : errmax = dum
335 732 : i_errmax = i
336 : end if
337 : end do
338 :
339 198 : errmax = errmax/eps
340 198 : km = k - 1
341 198 : err(km) = pow(errmax/safe1,1.0d0/(2*km+1))
342 : end if
343 :
344 : ! in order window
345 264 : if (k /= 1 .and. (k >= kopt-1 .or. first)) then
346 :
347 : ! converged
348 114 : if (errmax < 1.0d0) GOTO 04
349 :
350 : ! possible step size reductions
351 60 : if (k == kmax .or. k == kopt + 1) then
352 12 : red = safe2/err(km)
353 12 : GOTO 03
354 48 : else if (k == kopt) then
355 8 : if (alf(kopt-1,kopt) < err(km)) then
356 0 : red = 1.0d0/err(km)
357 0 : GOTO 03
358 : end if
359 40 : else if (kopt == kmax) then
360 12 : if (alf(km,kmax-1) < err(km)) then
361 0 : red = alf(km,kmax-1) * safe2/err(km)
362 0 : GOTO 03
363 : end if
364 28 : else if (alf(km,kopt) < err(km)) then
365 0 : red = alf(km,kopt-1)/err(km)
366 0 : GOTO 03
367 : end if
368 : end if
369 : end do
370 :
371 : ! reduce stepsize by at least redmin and at most redmax
372 12 : 03 red = min(red,redmin)
373 12 : red = max(red,redmax)
374 12 : h = h * red
375 12 : i = i_errmax
376 : if (dbg) write(*,*) 'reduce stepsize', i, errmax, yerr(i), yscal(i), red, h
377 12 : i_errmax = 0
378 12 : reduct = .true.
379 66 : GOTO 02
380 :
381 :
382 : ! successful step; get optimal row for convergence and corresponding stepsize
383 54 : 04 x = xnew
384 54 : hdid = h
385 54 : first = .false.
386 54 : wrkmin = 1.0d35
387 188 : do kk=1,km
388 134 : fact = max(err(kk),scalmx)
389 134 : work = fact * a(kk+1)
390 188 : if (work < wrkmin) then
391 134 : scale = fact
392 134 : wrkmin = work
393 134 : kopt = kk + 1
394 : end if
395 : end do
396 : !
397 : ! check for possible order increase, but not if stepsize was just reduced
398 54 : hnext = h/scale
399 54 : if (kopt >= k .and. kopt /= kmax .and. .not.reduct) then
400 44 : fact = max(scale/alf(kopt-1,kopt),scalmx)
401 44 : if (a(kopt+1)*fact <= wrkmin) then
402 20 : hnext = h/fact
403 20 : kopt = kopt + 1
404 : end if
405 : end if
406 : return
407 : end subroutine stifbs
408 :
409 :
410 264 : subroutine simpr( &
411 264 : y,dydx,nvar,xs,htot,nstep,yout,dens_dfdy,dmat, &
412 : derivs,ierr)
413 : !
414 : real(dp) :: dens_dfdy(:,:),dmat(:,:)
415 :
416 : interface
417 : include 'burner_derivs.inc'
418 : end interface
419 :
420 : integer, intent(out) :: ierr
421 :
422 : integer :: nvar,nstep
423 : integer :: i,j,nn,ii
424 : real(dp) :: y(:),dydx(:),xs,htot, &
425 11352 : yout(:),h,x,del(nvar),ytemp(nvar)
426 :
427 : !..for the linear algebra
428 264 : integer, target :: indx_a(nvar)
429 : integer, pointer :: indx(:)
430 :
431 : include 'formats'
432 :
433 264 : indx => indx_a
434 :
435 : ! stepsize this trip, and make the a matrix
436 264 : h = htot/nstep
437 5808 : do j=1,nvar
438 122232 : do i=1,nvar
439 121968 : dmat(i,j) = -h * dens_dfdy(i,j)
440 : end do
441 : end do
442 5808 : do i=1,nvar
443 5808 : dmat(i,i) = 1.0d0 + dmat(i,i)
444 : end do
445 :
446 : !..factor the matrix
447 264 : call my_getf2(nvar, dmat, nvar, indx, ierr)
448 264 : if (ierr /= 0) then
449 : if (dbg) write(*,*) 'my_getf2 failed in simpr'
450 : return
451 : end if
452 :
453 : ! use yout as temporary storage; the first step
454 5808 : do i=1,nvar
455 5544 : yout(i) = h * dydx(i)
456 264 : if (dbg) then
457 : if (is_bad(yout(i))) then
458 : write(*,*) 'bad yout in simpr nstep i yout', nstep, i, yout(i), dydx(i)
459 : call mesa_error(__FILE__,__LINE__,'simpr')
460 : end if
461 : end if
462 : end do
463 :
464 264 : call my_getrs1(nvar, dmat, nvar, indx, yout, nvar, ierr)
465 264 : if (ierr /= 0) then
466 : if (dbg) write(*,*) 'my_getrs1 failed in simpr'
467 : return
468 : end if
469 :
470 5808 : do i=1,nvar
471 5544 : del(i) = yout(i)
472 5544 : ytemp(i) = y(i) + del(i)
473 264 : if (dbg) then
474 : if (is_bad(ytemp(i))) then
475 :
476 : do j=1,nvar
477 : do ii=1,nvar
478 : if (dens_dfdy(ii,j) /= 0) write(*,3) 'dens_dfdy(ii,j)', ii, j, dens_dfdy(ii,j)
479 : end do
480 : end do
481 :
482 : do ii=1,nvar
483 : if (dydx(ii) /= 0) write(*,2) 'dydx(ii)', ii, dydx(ii)
484 : end do
485 :
486 : do j=1,nvar
487 : do ii=1,nvar
488 : if (dmat(ii,j) /= 0) write(*,3) 'dmat(ii,j)', ii, j, dmat(ii,j)
489 : end do
490 : end do
491 :
492 : do ii=1,nvar
493 : if (yout(ii) /= 0) write(*,2) 'yout(ii)', ii, yout(ii)
494 : end do
495 :
496 : write(*,*) 'first step: bad ytemp in simpr nstep i ytemp', nstep, i, ytemp(i), del(i), y(i)
497 : call mesa_error(__FILE__,__LINE__,'simpr')
498 :
499 : end if
500 :
501 : end if
502 : end do
503 :
504 264 : x = xs + h
505 264 : call derivs(x,ytemp,yout,nvar,ierr)
506 264 : if (ierr /= 0) then
507 : if (dbg) write(*,*) 'init derivs failed in simpr'
508 : return
509 : end if
510 :
511 : ! use yout as temporary storage; general step
512 :
513 3008 : do nn=2,nstep
514 60368 : do i=1,nvar
515 57624 : yout(i) = h*yout(i) - del(i)
516 2744 : if (dbg) then
517 : if (is_bad(yout(i))) then
518 : write(*,*) 'bad yout in simpr nn i yout', nn, i, yout(i)
519 : call mesa_error(__FILE__,__LINE__,'simpr')
520 : end if
521 : end if
522 : end do
523 2744 : call my_getrs1(nvar, dmat, nvar, indx, yout, nvar, ierr)
524 2744 : if (ierr /= 0) then
525 : if (dbg) write(*,*) 'my_getrs1 failed in simpr'
526 : return
527 : end if
528 60368 : do i=1,nvar
529 57624 : del(i) = del(i) + 2.0d0 * yout(i)
530 57624 : ytemp(i) = ytemp(i) + del(i)
531 2744 : if (dbg) then
532 : if (is_bad(ytemp(i))) then
533 : write(*,*) 'general step: bad ytemp in simpr nn i yout', nn, i, ytemp(i)
534 : call mesa_error(__FILE__,__LINE__,'simpr')
535 : end if
536 : end if
537 : end do
538 :
539 2744 : x = x + h
540 2744 : call derivs(x,ytemp,yout,nvar,ierr)
541 3008 : if (ierr /= 0) then
542 : if (dbg) write(*,*) 'derivs failed in simpr general step'
543 : return
544 : end if
545 : end do
546 :
547 : ! take the last step
548 5808 : do i=1,nvar
549 5544 : yout(i) = h * yout(i) - del(i)
550 264 : if (dbg) then
551 : if (is_bad(yout(i))) then
552 : write(*,*) 'bad yout in simpr last step: nstep i yout', nstep, i, yout(i)
553 : call mesa_error(__FILE__,__LINE__,'simpr')
554 : end if
555 : end if
556 : end do
557 264 : call my_getrs1(nvar, dmat, nvar, indx, yout, nvar, ierr)
558 264 : if (ierr /= 0) then
559 0 : write(*,*) 'my_getrs1 failed in simpr'
560 0 : return
561 : end if
562 :
563 5808 : do i=1,nvar
564 5544 : yout(i) = ytemp(i) + yout(i)
565 264 : if (dbg) then
566 : if (is_bad(yout(i))) then
567 : write(*,*) 'bad yout in simpr result: nstep i yout', nstep, i, yout(i)
568 : call mesa_error(__FILE__,__LINE__,'simpr')
569 : end if
570 : end if
571 : end do
572 :
573 : return
574 : end subroutine simpr
575 :
576 :
577 264 : subroutine net_pzextr(iest,xest,yest,yz,dy,nvar,qcol,x)
578 : ! use polynomial extrapolation to evaluate nvar functions at x=0 by fitting
579 : ! a polynomial to a sequence of estimates with progressively smaller values
580 : ! x=xest, and corresponding function vectors yest(1:nvar). the call is number
581 : ! iest in the sequence of calls. extrapolated function values are output as
582 : ! yz(1:nvar), and their estimated error is output as dy(1:nvar)
583 :
584 : ! declare the pass
585 : integer :: iest,nvar
586 : real(dp) :: xest,dy(:),yest(:),yz(:),qcol(:,:),x(:)
587 :
588 : ! locals; qcol and x must be "saved" between successive calls
589 :
590 : integer :: j,k1
591 5808 : real(dp) :: delta,f1,f2,q,d(nvar)
592 :
593 :
594 : ! sanity checks
595 :
596 264 : if (iest > qcol_imax) call mesa_error(__FILE__,__LINE__,'iest > qcol_imax in net_pzextr')
597 :
598 : ! save current independent variables
599 264 : x(iest) = xest
600 5808 : do j=1,nvar
601 5544 : dy(j) = yest(j)
602 5808 : yz(j) = yest(j)
603 : end do
604 :
605 : ! store first estimate in first column
606 264 : if (iest == 1) then
607 1452 : do j=1,nvar
608 1452 : qcol(j,1) = yest(j)
609 : end do
610 : else
611 4356 : do j=1,nvar
612 4356 : d(j) = yest(j)
613 : end do
614 696 : do k1=1,iest-1
615 498 : delta = 1.0d0/(x(iest-k1) - xest)
616 498 : f1 = xest * delta
617 498 : f2 = x(iest-k1) * delta
618 :
619 : ! propagate tableu 1 diagonal more
620 11154 : do j=1,nvar
621 10458 : q = qcol(j,k1)
622 10458 : qcol(j,k1) = dy(j)
623 10458 : delta = d(j) - q
624 10458 : dy(j) = f1*delta
625 10458 : d(j) = f2*delta
626 10956 : yz(j) = yz(j) + dy(j)
627 : end do
628 : end do
629 4356 : do j=1,nvar
630 4356 : qcol(j,iest) = dy(j)
631 : end do
632 : end if
633 264 : return
634 : end subroutine net_pzextr
635 :
636 :
637 : include 'mtx_solve_routines.inc'
638 :
639 :
640 : end module net_burn_support
641 :
|
