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