diff --git a/net/private/net_approx21.f90 b/net/private/net_approx21.f90 index 2a0f5e15a9..96b672174a 100644 --- a/net/private/net_approx21.f90 +++ b/net/private/net_approx21.f90 @@ -34,6 +34,15 @@ module net_approx21 integer, parameter :: species_21 = 21, species_co56 = 22 + integer, parameter :: num_reduced_flows = 7, & + i_fe52_fe54_n = 1, & + i_fe54_fe56_n = 2, & + i_fe54_fe56_ap = 3, & + i_fe54_ni56 = 4, & + i_fe52_fe54 = 5, & + i_fe52_ni56 = 6, & + i_he4_breakup = 7 + integer :: iso_cid(species_co56) ! these are corresponding chem ids for the isos ! e.g., iso_cid(ife52) is = the iso number for fe52 as defined in mesa/chem @@ -923,11 +932,287 @@ end subroutine turn_off_reaction end subroutine approx21_special_reactions + ! Keep reduced forward and reverse flows paired through cancellation. + pure subroutine eval_approx21_reduced_flows( & + y, rate, temp, reduced_dydt, reduced_dfdy, & + drate, reduced_drate) + real(qp), intent(in) :: y(:), rate(:), temp + real(qp), intent(out) :: reduced_dydt(size(y)) + real(qp), intent(out), optional :: reduced_dfdy(size(y),size(y)) + real(qp), intent(in), optional :: drate(:) + real(qp), intent(out), optional :: reduced_drate(size(y)) + + integer :: j + real(qp) :: denom, denom_drate, denom_dyneut, denom_dyprot + real(qp) :: numer, numer_drate + real(qp) :: flows(num_reduced_flows) + real(qp) :: dflows_dy(num_reduced_flows,size(y)) + real(qp) :: dflows_drate(num_reduced_flows) + real(qp) :: forward_product, reverse_product + real(qp), parameter :: tiny_denom = 1d-50 + + reduced_dydt = 0d0 + flows = 0d0 + dflows_dy = 0d0 + dflows_drate = 0d0 + if (present(reduced_dfdy)) reduced_dfdy = 0d0 + if (present(reduced_drate)) reduced_drate = 0d0 + + if (.not. (temp > 1.5d9)) return + + ! fe52(n,g)fe53(n,g)fe54 + denom = rate(ir53gn) + y(ineut)*rate(ir53ng) + if (denom > tiny_denom) then + reverse_product = rate(ir54gn)*rate(ir53gn) + forward_product = rate(ir52ng)*rate(ir53ng) + numer = y(ife54)*reverse_product - & + y(ife52)*y(ineut)*y(ineut)*forward_product + flows(i_fe52_fe54_n) = numer/denom + dflows_dy(i_fe52_fe54_n,ife52) = & + -y(ineut)*y(ineut)*forward_product/denom + dflows_dy(i_fe52_fe54_n,ife54) = reverse_product/denom + dflows_dy(i_fe52_fe54_n,ineut) = ( & + -2d0*y(ife52)*y(ineut)*forward_product - & + flows(i_fe52_fe54_n)*rate(ir53ng))/denom + if (present(drate)) then + denom_drate = drate(ir53gn) + y(ineut)*drate(ir53ng) + numer_drate = y(ife54)*( & + drate(ir54gn)*rate(ir53gn) + rate(ir54gn)*drate(ir53gn)) - & + y(ife52)*y(ineut)*y(ineut)*( & + drate(ir52ng)*rate(ir53ng) + rate(ir52ng)*drate(ir53ng)) + dflows_drate(i_fe52_fe54_n) = (numer_drate - & + flows(i_fe52_fe54_n)*denom_drate)/denom + end if + end if + + ! fe54(n,g)fe55(n,g)fe56 + denom = rate(ir55gn) + y(ineut)*rate(ir55ng) + if (denom > tiny_denom) then + reverse_product = rate(ir56gn)*rate(ir55gn) + forward_product = rate(ir54ng)*rate(ir55ng) + numer = y(ife56)*reverse_product - & + y(ife54)*y(ineut)*y(ineut)*forward_product + flows(i_fe54_fe56_n) = numer/denom + dflows_dy(i_fe54_fe56_n,ife54) = & + -y(ineut)*y(ineut)*forward_product/denom + dflows_dy(i_fe54_fe56_n,ife56) = reverse_product/denom + dflows_dy(i_fe54_fe56_n,ineut) = ( & + -2d0*y(ife54)*y(ineut)*forward_product - & + flows(i_fe54_fe56_n)*rate(ir55ng))/denom + if (present(drate)) then + denom_drate = drate(ir55gn) + y(ineut)*drate(ir55ng) + numer_drate = y(ife56)*( & + drate(ir56gn)*rate(ir55gn) + rate(ir56gn)*drate(ir55gn)) - & + y(ife54)*y(ineut)*y(ineut)*( & + drate(ir54ng)*rate(ir55ng) + rate(ir54ng)*drate(ir55ng)) + dflows_drate(i_fe54_fe56_n) = (numer_drate - & + flows(i_fe54_fe56_n)*denom_drate)/denom + end if + end if + + ! fe54(a,p)co57(g,p)fe56 + denom = rate(irco57gp) + y(iprot)*rate(irco57pa) + if (denom > tiny_denom) then + reverse_product = rate(irfe56pg)*rate(irco57pa) + forward_product = rate(irfe54ap)*rate(irco57gp) + numer = y(ife56)*y(iprot)*y(iprot)*reverse_product - & + y(ife54)*y(ihe4)*forward_product + flows(i_fe54_fe56_ap) = numer/denom + dflows_dy(i_fe54_fe56_ap,ihe4) = -y(ife54)*forward_product/denom + dflows_dy(i_fe54_fe56_ap,ife54) = -y(ihe4)*forward_product/denom + dflows_dy(i_fe54_fe56_ap,ife56) = & + y(iprot)*y(iprot)*reverse_product/denom + dflows_dy(i_fe54_fe56_ap,iprot) = ( & + 2d0*y(ife56)*y(iprot)*reverse_product - & + flows(i_fe54_fe56_ap)*rate(irco57pa))/denom + if (present(drate)) then + denom_drate = drate(irco57gp) + y(iprot)*drate(irco57pa) + numer_drate = y(ife56)*y(iprot)*y(iprot)*( & + drate(irfe56pg)*rate(irco57pa) + & + rate(irfe56pg)*drate(irco57pa)) - & + y(ife54)*y(ihe4)*( & + drate(irfe54ap)*rate(irco57gp) + & + rate(irfe54ap)*drate(irco57gp)) + dflows_drate(i_fe54_fe56_ap) = (numer_drate - & + flows(i_fe54_fe56_ap)*denom_drate)/denom + end if + end if + + ! Links through the eliminated co55 intermediate. + denom = rate(ircogp) + y(iprot)*(rate(ircopg) + rate(ircopa)) + if (denom > tiny_denom) then + denom_dyprot = rate(ircopg) + rate(ircopa) + + reverse_product = rate(irnigp)*rate(ircogp) + forward_product = rate(irfepg)*rate(ircopg) + numer = y(ini56)*reverse_product - & + y(ife54)*y(iprot)*y(iprot)*forward_product + flows(i_fe54_ni56) = numer/denom + dflows_dy(i_fe54_ni56,ife54) = & + -y(iprot)*y(iprot)*forward_product/denom + dflows_dy(i_fe54_ni56,ini56) = reverse_product/denom + dflows_dy(i_fe54_ni56,iprot) = ( & + -2d0*y(ife54)*y(iprot)*forward_product - & + flows(i_fe54_ni56)*denom_dyprot)/denom + + reverse_product = rate(irfepg)*rate(ircopa) + forward_product = rate(irfeap)*rate(ircogp) + numer = y(ife54)*y(iprot)*y(iprot)*reverse_product - & + y(ife52)*y(ihe4)*forward_product + flows(i_fe52_fe54) = numer/denom + dflows_dy(i_fe52_fe54,ihe4) = -y(ife52)*forward_product/denom + dflows_dy(i_fe52_fe54,ife52) = -y(ihe4)*forward_product/denom + dflows_dy(i_fe52_fe54,ife54) = & + y(iprot)*y(iprot)*reverse_product/denom + dflows_dy(i_fe52_fe54,iprot) = ( & + 2d0*y(ife54)*y(iprot)*reverse_product - & + flows(i_fe52_fe54)*denom_dyprot)/denom + + reverse_product = rate(irnigp)*rate(ircopa) + forward_product = rate(irfeap)*rate(ircopg) + numer = y(ini56)*y(iprot)*reverse_product - & + y(ife52)*y(ihe4)*y(iprot)*forward_product + flows(i_fe52_ni56) = numer/denom + dflows_dy(i_fe52_ni56,ihe4) = & + -y(ife52)*y(iprot)*forward_product/denom + dflows_dy(i_fe52_ni56,ife52) = & + -y(ihe4)*y(iprot)*forward_product/denom + dflows_dy(i_fe52_ni56,ini56) = y(iprot)*reverse_product/denom + dflows_dy(i_fe52_ni56,iprot) = ( & + y(ini56)*reverse_product - y(ife52)*y(ihe4)*forward_product - & + flows(i_fe52_ni56)*denom_dyprot)/denom + + if (present(drate)) then + denom_drate = drate(ircogp) + & + y(iprot)*(drate(ircopg) + drate(ircopa)) + + numer_drate = y(ini56)*( & + drate(irnigp)*rate(ircogp) + rate(irnigp)*drate(ircogp)) - & + y(ife54)*y(iprot)*y(iprot)*( & + drate(irfepg)*rate(ircopg) + rate(irfepg)*drate(ircopg)) + dflows_drate(i_fe54_ni56) = (numer_drate - & + flows(i_fe54_ni56)*denom_drate)/denom + + numer_drate = y(ife54)*y(iprot)*y(iprot)*( & + drate(irfepg)*rate(ircopa) + rate(irfepg)*drate(ircopa)) - & + y(ife52)*y(ihe4)*( & + drate(irfeap)*rate(ircogp) + rate(irfeap)*drate(ircogp)) + dflows_drate(i_fe52_fe54) = (numer_drate - & + flows(i_fe52_fe54)*denom_drate)/denom + + numer_drate = y(ini56)*y(iprot)*( & + drate(irnigp)*rate(ircopa) + rate(irnigp)*drate(ircopa)) - & + y(ife52)*y(ihe4)*y(iprot)*( & + drate(irfeap)*rate(ircopg) + rate(irfeap)*drate(ircopg)) + dflows_drate(i_fe52_ni56) = (numer_drate - & + flows(i_fe52_ni56)*denom_drate)/denom + end if + end if + + ! he4 photodisintegration and light-particle reconstruction. + denom = rate(irhegp)*rate(irdgn) + & + y(ineut)*rate(irheng)*rate(irdgn) + & + y(ineut)*y(iprot)*rate(irheng)*rate(irdpg) + if (denom > tiny_denom) then + reverse_product = rate(irhegn)*rate(irhegp)*rate(irdgn) + forward_product = rate(irheng)*rate(irdpg)*rate(irhng) + numer = y(ihe4)*reverse_product - & + y(ineut)*y(ineut)*y(iprot)*y(iprot)*forward_product + flows(i_he4_breakup) = numer/denom + denom_dyneut = rate(irheng)*rate(irdgn) + & + y(iprot)*rate(irheng)*rate(irdpg) + denom_dyprot = y(ineut)*rate(irheng)*rate(irdpg) + dflows_dy(i_he4_breakup,ihe4) = reverse_product/denom + dflows_dy(i_he4_breakup,ineut) = ( & + -2d0*y(ineut)*y(iprot)*y(iprot)*forward_product - & + flows(i_he4_breakup)*denom_dyneut)/denom + dflows_dy(i_he4_breakup,iprot) = ( & + -2d0*y(ineut)*y(ineut)*y(iprot)*forward_product - & + flows(i_he4_breakup)*denom_dyprot)/denom + if (present(drate)) then + denom_drate = & + drate(irhegp)*rate(irdgn) + rate(irhegp)*drate(irdgn) + & + y(ineut)*( & + drate(irheng)*rate(irdgn) + rate(irheng)*drate(irdgn)) + & + y(ineut)*y(iprot)*( & + drate(irheng)*rate(irdpg) + rate(irheng)*drate(irdpg)) + numer_drate = y(ihe4)*( & + drate(irhegn)*rate(irhegp)*rate(irdgn) + & + rate(irhegn)*drate(irhegp)*rate(irdgn) + & + rate(irhegn)*rate(irhegp)*drate(irdgn)) - & + y(ineut)*y(ineut)*y(iprot)*y(iprot)*( & + drate(irheng)*rate(irdpg)*rate(irhng) + & + rate(irheng)*drate(irdpg)*rate(irhng) + & + rate(irheng)*rate(irdpg)*drate(irhng)) + dflows_drate(i_he4_breakup) = (numer_drate - & + flows(i_he4_breakup)*denom_drate)/denom + end if + end if + + reduced_dydt(ihe4) = flows(i_fe52_fe54) + flows(i_fe52_ni56) - & + flows(i_he4_breakup) + flows(i_fe54_fe56_ap) + reduced_dydt(ife52) = flows(i_fe52_fe54_n) + & + flows(i_fe52_fe54) + flows(i_fe52_ni56) + reduced_dydt(ife54) = -flows(i_fe52_fe54_n) + & + flows(i_fe54_fe56_n) + flows(i_fe54_fe56_ap) + & + flows(i_fe54_ni56) - flows(i_fe52_fe54) + reduced_dydt(ife56) = -flows(i_fe54_fe56_n) - flows(i_fe54_fe56_ap) + reduced_dydt(ini56) = -flows(i_fe54_ni56) - flows(i_fe52_ni56) + reduced_dydt(ineut) = 2d0*( & + flows(i_fe52_fe54_n) + flows(i_fe54_fe56_n) + flows(i_he4_breakup)) + reduced_dydt(iprot) = 2d0*( & + flows(i_fe54_ni56) - flows(i_fe52_fe54) + & + flows(i_he4_breakup) - flows(i_fe54_fe56_ap)) + + if (present(reduced_dfdy)) then + do j = 1, size(y) + reduced_dfdy(ihe4,j) = dflows_dy(i_fe52_fe54,j) + & + dflows_dy(i_fe52_ni56,j) - dflows_dy(i_he4_breakup,j) + & + dflows_dy(i_fe54_fe56_ap,j) + reduced_dfdy(ife52,j) = dflows_dy(i_fe52_fe54_n,j) + & + dflows_dy(i_fe52_fe54,j) + dflows_dy(i_fe52_ni56,j) + reduced_dfdy(ife54,j) = -dflows_dy(i_fe52_fe54_n,j) + & + dflows_dy(i_fe54_fe56_n,j) + dflows_dy(i_fe54_fe56_ap,j) + & + dflows_dy(i_fe54_ni56,j) - dflows_dy(i_fe52_fe54,j) + reduced_dfdy(ife56,j) = -dflows_dy(i_fe54_fe56_n,j) - & + dflows_dy(i_fe54_fe56_ap,j) + reduced_dfdy(ini56,j) = -dflows_dy(i_fe54_ni56,j) - & + dflows_dy(i_fe52_ni56,j) + reduced_dfdy(ineut,j) = 2d0*(dflows_dy(i_fe52_fe54_n,j) + & + dflows_dy(i_fe54_fe56_n,j) + dflows_dy(i_he4_breakup,j)) + reduced_dfdy(iprot,j) = 2d0*(dflows_dy(i_fe54_ni56,j) - & + dflows_dy(i_fe52_fe54,j) + dflows_dy(i_he4_breakup,j) - & + dflows_dy(i_fe54_fe56_ap,j)) + end do + end if + + if (present(reduced_drate)) then + reduced_drate(ihe4) = dflows_drate(i_fe52_fe54) + & + dflows_drate(i_fe52_ni56) - dflows_drate(i_he4_breakup) + & + dflows_drate(i_fe54_fe56_ap) + reduced_drate(ife52) = dflows_drate(i_fe52_fe54_n) + & + dflows_drate(i_fe52_fe54) + dflows_drate(i_fe52_ni56) + reduced_drate(ife54) = -dflows_drate(i_fe52_fe54_n) + & + dflows_drate(i_fe54_fe56_n) + dflows_drate(i_fe54_fe56_ap) + & + dflows_drate(i_fe54_ni56) - dflows_drate(i_fe52_fe54) + reduced_drate(ife56) = -dflows_drate(i_fe54_fe56_n) - & + dflows_drate(i_fe54_fe56_ap) + reduced_drate(ini56) = -dflows_drate(i_fe54_ni56) - & + dflows_drate(i_fe52_ni56) + reduced_drate(ineut) = 2d0*(dflows_drate(i_fe52_fe54_n) + & + dflows_drate(i_fe54_fe56_n) + dflows_drate(i_he4_breakup)) + reduced_drate(iprot) = 2d0*(dflows_drate(i_fe54_ni56) - & + dflows_drate(i_fe52_fe54) + dflows_drate(i_he4_breakup) - & + dflows_drate(i_fe54_fe56_ap)) + end if + end subroutine eval_approx21_reduced_flows + + subroutine approx21_dydt( & - y, rate, ratdum, dydt, deriva, & + y_in, rate_in, ratdum_in, dydt, deriva, & fe56ec_fake_factor_in, min_T, fe56ec_n_neut, temp, den, plus_co56, ierr) logical, intent(in) :: deriva ! false for dydt, true for partials wrt T, Rho - real(dp), dimension(:), intent(in) :: y, rate, ratdum + real(dp), dimension(:), intent(in) :: y_in, rate_in, ratdum_in integer, intent(in) :: fe56ec_n_neut real(dp), dimension(:), intent(out) :: dydt real(dp), intent(in) :: fe56ec_fake_factor_in, temp, den @@ -941,6 +1226,8 @@ subroutine approx21_dydt( & real(qp) :: a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,& a11,a12,a13,a14,a15,a16,a17,a18,a19,a20 real(qp) :: qray(species(plus_co56)) + real(qp) :: y(size(y_in)), rate(size(rate_in)), ratdum(size(ratdum_in)) + real(qp) :: reduced_dydt(species(plus_co56)), reduced_drate(species(plus_co56)) logical :: okay @@ -948,6 +1235,10 @@ subroutine approx21_dydt( & ierr = 0 + y = real(y_in,kind=qp) + rate = real(rate_in,kind=qp) + ratdum = real(ratdum_in,kind=qp) + ! Turn on special fe56ec rate above some temperature fe56ec_fake_factor = 0d0 if(.not.deriva) then @@ -957,6 +1248,15 @@ subroutine approx21_dydt( & dydt(1:species(plus_co56)) = 0.0d0 qray(1:species(plus_co56)) = 0.0_qp + if (deriva) then + call eval_approx21_reduced_flows( & + y, ratdum, real(temp,kind=qp), reduced_dydt, & + drate=rate, reduced_drate=reduced_drate) + else + call eval_approx21_reduced_flows( & + y, ratdum, real(temp,kind=qp), reduced_dydt) + end if + ! hydrogen reactions a1 = -1.5d0 * y(ih1) * y(ih1) * rate(irpp) a2 = y(ihe3) * y(ihe3) * rate(ir33) @@ -1086,20 +1386,6 @@ subroutine approx21_dydt( & qray(ihe4) = qray(ihe4) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 end if - - ! photodisintegration reactions - a1 = y(ife54) * y(iprot) * y(iprot) * rate(ir5f54) - a2 = -y(ife52) * y(ihe4) * rate(ir6f54) - a3 = -y(ife52) * y(ihe4) * y(iprot) * rate(ir7f54) - a4 = y(ini56) * y(iprot) * rate(ir8f54) - a5 = -y(ihe4) * rate(iralf1) - a6 = y(ineut)*y(ineut) * y(iprot)*y(iprot) * rate(iralf2) - a7 = y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3) - a8 = -y(ife54) * y(ihe4) * rate(irfe56_aux4) - - qray(ihe4) = qray(ihe4) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 - - ! ppchain a1 = 0.5d0 * y(ihe3) * y(ihe3) * rate(ir33) a2 = y(ihe3) * y(ihe4) * rate(irhe3ag) @@ -1431,32 +1717,6 @@ subroutine approx21_dydt( & qray(ife52) = qray(ife52) + a1 + a2 + a3 + a4 end if - a1 = y(ife54) * rate(ir1f54) - a2 = -y(ife52) * y(ineut) * y(ineut) * rate(ir2f54) - a3 = y(ife54) * y(iprot) * y(iprot) * rate(ir5f54) - a4 = -y(ife52) * y(ihe4) * rate(ir6f54) - a5 = -y(ife52) * y(ihe4) * y(iprot) * rate(ir7f54) - a6 = y(ini56) * y(iprot) * rate(ir8f54) - - qray(ife52) = qray(ife52) + a1 + a2 + a3 + a4 + a5 + a6 - - - ! fe54 reactions - a1 = -y(ife54) * rate(ir1f54) - a2 = y(ife52) * y(ineut) * y(ineut) * rate(ir2f54) - a3 = -y(ife54) * y(iprot) * y(iprot) * rate(ir3f54) - a4 = y(ini56) * rate(ir4f54) - a5 = -y(ife54) * y(iprot) * y(iprot) * rate(ir5f54) - a6 = y(ife52) * y(ihe4) * rate(ir6f54) - a7 = y(ife56) * rate(irfe56_aux1) - a8 = -y(ife54) * y(ineut) * y(ineut) * rate(irfe56_aux2) - a9 = y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3) - a10 = -y(ife54) * y(ihe4) * rate(irfe56_aux4) - - qray(ife54) = qray(ife54) + & - a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10 - - ! fe56 reactions if (plus_co56) then a1 = y(ico56) * rate(irco56ec) @@ -1464,12 +1724,8 @@ subroutine approx21_dydt( & a1 = y(ini56) * rate(irn56ec) end if a2 = -y(ife56) * fe56ec_fake_factor * rate(irn56ec) - a3 = -y(ife56) * rate(irfe56_aux1) - a4 = y(ife54) * y(ineut) * y(ineut) * rate(irfe56_aux2) - a5 = -y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3) - a6 = y(ife54) * y(ihe4) * rate(irfe56_aux4) - qray(ife56) = qray(ife56) + a1 + a2 + a3 + a4 + a5 + a6 + qray(ife56) = qray(ife56) + a1 + a2 if (plus_co56) then ! co56 reactions @@ -1486,47 +1742,31 @@ subroutine approx21_dydt( & qray(ini56) = qray(ini56) + a1 + a2 + a3 - a1 = y(ife54) * y(iprot) * y(iprot) * rate(ir3f54) - a2 = -y(ini56) * rate(ir4f54) - a3 = y(ife52) * y(ihe4)* y(iprot) * rate(ir7f54) - a4 = -y(ini56) * y(iprot) * rate(ir8f54) - - qray(ini56) = qray(ini56) + a1 + a2 + a3 + a4 - ! neutrons - a1 = 2.0d0 * y(ife54) * rate(ir1f54) - a2 = -2.0d0 * y(ife52) * y(ineut) * y(ineut) * rate(ir2f54) - a3 = 2.0d0 * y(ihe4) * rate(iralf1) - a4 = -2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * rate(iralf2) a5 = y(iprot) * rate(irpen) a6 = -y(ineut) * rate(irnep) - a7 = 2.0d0 * y(ife56) * rate(irfe56_aux1) - a8 = -2.0d0 * y(ife54) * y(ineut) * y(ineut) * rate(irfe56_aux2) a9 = -fe56ec_n_neut * y(ife56) * fe56ec_fake_factor * rate(irn56ec) - qray(ineut) = qray(ineut) + a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + qray(ineut) = qray(ineut) + a5 + a6 + a9 ! photodisintegration protons - a1 = -2.0d0 * y(ife54) * y(iprot) * y(iprot) * rate(ir3f54) - a2 = 2.0d0 * y(ini56) * rate(ir4f54) - a3 = -2.0d0 * y(ife54) * y(iprot) * y(iprot) * rate(ir5f54) - a4 = 2.0d0 * y(ife52) * y(ihe4) * rate(ir6f54) - a5 = 2.0d0 * y(ihe4) * rate(iralf1) - a6 = -2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * rate(iralf2) a7 = -y(iprot) * rate(irpen) a8 = y(ineut) * rate(irnep) - a9 = -2.0d0 * y(ife56) * y(iprot) * y(iprot) * rate(irfe56_aux3) - a10 = 2.0d0 * y(ife54) * y(ihe4) * rate(irfe56_aux4) - qray(iprot) = qray(iprot) + & - a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8 + a9 + a10 + qray(iprot) = qray(iprot) + a7 + a8 + + if (deriva) then + qray = qray + reduced_drate + else + qray = qray + reduced_dydt + end if ! now set the real(dp) return argument dydt okay = .true. do i=1,species(plus_co56) dydt(i) = qray(i) if (is_bad(dydt(i))) then - write(*,*) 'dydt(i)', i, dydt(i), y(i) + write(*,*) 'dydt(i)', i, dydt(i), y_in(i) okay = .false. end if end do @@ -1536,7 +1776,7 @@ subroutine approx21_dydt( & write(*,*) 'log10(rho) = ',safe_log10(den) do i=1,num_reactions(plus_co56) - write(*,*) trim(ratnam(i)), i, rate(i) + write(*,*) trim(ratnam(i)), i, rate_in(i) end do call mesa_error(__FILE__,__LINE__,'approx21_dydt') end if @@ -2045,12 +2285,20 @@ subroutine approx21_dfdy( & real(dp) :: abar,zbar,ye,taud,taut, b1, & snuda,snudz,enuc,velx,posx,zz real(dp) :: fe56ec_fake_factor + real(qp) :: y_qp(size(y)), rate_qp(size(ratdum)) + real(qp) :: reduced_dydt_qp(size(y)), reduced_dfdy_qp(size(y),size(y)) ierr = 0 ! Turn on special fe56ec rate above some temperature fe56ec_fake_factor=eval_fe56ec_fake_factor(fe56ec_fake_factor_in,min_T,btemp) + y_qp = real(y,kind=qp) + rate_qp = real(ratdum,kind=qp) + call eval_approx21_reduced_flows( & + y_qp, rate_qp, real(btemp,kind=qp), reduced_dydt_qp, & + reduced_dfdy=reduced_dfdy_qp) + ! NOTE: use of quad precision for dfdy doesn't make a difference. dfdy(1:species(plus_co56),1:species(plus_co56)) = 0.0d0 @@ -2120,13 +2368,6 @@ subroutine approx21_dfdy( & - y(iti44) * ratdum(irtiap) * (1.0d0-ratdum(irw1)) & - y(icr48) * ratdum(ircrap) * (1.0d0-ratdum(irx1)) - dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) & - - y(ife52) * ratdum(ir6f54) & - - y(ife52) * y(iprot) * ratdum(ir7f54) & - - ratdum(iralf1) & - - y(ife54) * ratdum(irfe56_aux4) - - dfdy(ihe4,ihe4) = dfdy(ihe4,ihe4) & + y(ihe3) * ratdum(irhe3ag) & + y(ihe3) * y(ihe4) * dratdumdy1(irhe3ag) & @@ -2189,39 +2430,12 @@ subroutine approx21_dfdy( & dfdy(ihe4,ife52) = ratdum(irfega) & - y(ihe4) * ratdum(irfeag) & - + ratdum(irx1) * ratdum(irfegp) & - - y(ihe4) * ratdum(ir6f54) & - - y(ihe4) * y(iprot) * ratdum(ir7f54) - - dfdy(ihe4,ife54) = y(iprot) * y(iprot) * ratdum(ir5f54) & - - y(ihe4) * ratdum(irfe56_aux4) - - dfdy(ihe4,ife56) = y(iprot) * y(iprot) * ratdum(irfe56_aux3) - - dfdy(ihe4,ini56) = ratdum(irniga) & - + y(iprot) * ratdum(ir8f54) + + ratdum(irx1) * ratdum(irfegp) - - dfdy(ihe4,ineut) = -y(ihe4) * dratdumdy1(iralf1) & - + 2.0d0 * y(ineut) * y(iprot)*y(iprot) * ratdum(iralf2) & - + y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy1(iralf2) + dfdy(ihe4,ini56) = ratdum(irniga) include 'formats' - dfdy(ihe4,iprot) = 2.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) & - + y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir5f54) & - - y(ihe4) * y(ife52) * dratdumdy1(ir6f54) & - - y(ife52) * y(ihe4) * ratdum(ir7f54) & - - y(ife52) * y(ihe4) * y(iprot) * dratdumdy1(ir7f54) & - + y(ini56) * ratdum(ir8f54) & - + y(ini56) * y(iprot) * dratdumdy1(ir8f54) & - - y(ihe4) * dratdumdy2(iralf1) & - + 2.0d0 * y(ineut)*y(ineut) * y(iprot) * ratdum(iralf2) & - + y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy2(iralf2) & - + 2.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) & - + y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) & - - y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4) - ! c12 jacobian elements @@ -2453,92 +2667,19 @@ subroutine approx21_dfdy( & ! fe52 jacobian elements dfdy(ife52,ihe4) = y(icr48) * ratdum(ircrag) & - y(ife52) * ratdum(irfeag) & - + y(icr48) * ratdum(ircrap) * (1.0d0-ratdum(irx1)) & - - y(ife52) * ratdum(ir6f54) & - - y(ife52) * y(iprot) * ratdum(ir7f54) + + y(icr48) * ratdum(ircrap) * (1.0d0-ratdum(irx1)) dfdy(ife52,icr48) = y(ihe4) * ratdum(ircrag) & + y(ihe4) * ratdum(ircrap) * (1.0d0-ratdum(irx1)) dfdy(ife52,ife52) = - y(ihe4) * ratdum(irfeag) & - ratdum(irfega) & - - ratdum(irx1) * ratdum(irfegp) & - - y(ineut) * y(ineut) * ratdum(ir2f54) & - - y(ihe4) * ratdum(ir6f54) & - - y(ihe4) * y(iprot) * ratdum(ir7f54) - - dfdy(ife52,ife54) = ratdum(ir1f54) + & - y(iprot) * y(iprot) * ratdum(ir5f54) - - dfdy(ife52,ini56) = ratdum(irniga) & - + y(iprot) * ratdum(ir8f54) - - dfdy(ife52,ineut) = & - y(ife54) * dratdumdy1(ir1f54) & - - 2.0d0 * y(ife52) * y(ineut) * ratdum(ir2f54) & - - y(ife52) * y(ineut) * y(ineut) * dratdumdy1(ir2f54) - - dfdy(ife52,iprot) = 2.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) & - + y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir5f54) & - - y(ihe4) * y(ife52) * dratdumdy1(ir6f54) & - - y(ife52) * y(ihe4) * ratdum(ir7f54) & - - y(ife52) * y(ihe4) * y(iprot) * dratdumdy1(ir7f54) & - + y(ini56) * ratdum(ir8f54) & - + y(ini56) * y(iprot) * dratdumdy1(ir8f54) - - - ! fe54 jacobian elements - dfdy(ife54,ihe4) = y(ife52) * ratdum(ir6f54) & - - y(ife54) * ratdum(irfe56_aux4) - - dfdy(ife54,ife52) = & - y(ineut) * y(ineut) * ratdum(ir2f54) + & - y(ihe4) * ratdum(ir6f54) - - dfdy(ife54,ife54) = & - - ratdum(ir1f54) & - - y(ineut) * y(ineut) * ratdum(irfe56_aux2) & - - y(iprot) * y(iprot) * ratdum(ir3f54) & - - y(iprot) * y(iprot) * ratdum(ir5f54) & - - y(ihe4) * ratdum(irfe56_aux4) - - dfdy(ife54,ife56) = & - ratdum(irfe56_aux1) + & - y(iprot) * y(iprot) * ratdum(irfe56_aux3) - - dfdy(ife54,ini56) = ratdum(ir4f54) - - dfdy(ife54,ineut) = & - - y(ife54) * dratdumdy1(ir1f54) & - + 2.0d0 * y(ife52) * y(ineut) * ratdum(ir2f54) & - + y(ife52) * y(ineut) * y(ineut) * dratdumdy1(ir2f54) & - + y(ife56) * dratdumdy1(irfe56_aux1) & - - 2.0d0 * y(ife54) * y(ineut) * ratdum(irfe56_aux2) & - - y(ife54) * y(ineut) * y(ineut) * dratdumdy1(irfe56_aux2) - - dfdy(ife54,iprot) = -2.0d0 * y(ife54) * y(iprot) * ratdum(ir3f54) & - - y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir3f54) & - + y(ini56) * dratdumdy1(ir4f54) & - - 2.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) & - - y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir5f54) & - + y(ihe4) * y(ife52) * dratdumdy1(ir6f54) & - + 2.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) & - + y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) & - - y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4) + - ratdum(irx1) * ratdum(irfegp) + dfdy(ife52,ini56) = ratdum(irniga) ! fe56 jacobian elements - - dfdy(ife56,ihe4) = y(ife54) * ratdum(irfe56_aux4) - - - dfdy(ife56,ife54) = & - y(ineut) * y(ineut) * ratdum(irfe56_aux2) + & - y(ihe4) * ratdum(irfe56_aux4) - - dfdy(ife56,ife56) = - fe56ec_fake_factor * ratdum(irn56ec) & - - ratdum(irfe56_aux1) & - - y(iprot) * y(iprot) * ratdum(irfe56_aux3) + dfdy(ife56,ife56) = -fe56ec_fake_factor * ratdum(irn56ec) if (plus_co56) then dfdy(ife56,ico56) = ratdum(irco56ec) @@ -2546,17 +2687,6 @@ subroutine approx21_dfdy( & dfdy(ife56,ini56) = ratdum(irn56ec) end if - - dfdy(ife56,ineut) = & - -y(ife56) * dratdumdy1(irfe56_aux1) & - + 2.0d0 * y(ife54) * y(ineut) * ratdum(irfe56_aux2) & - + y(ife54) * y(ineut) * y(ineut) * dratdumdy1(irfe56_aux2) - - - dfdy(ife56,iprot) = -2.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) & - - y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) & - + y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4) - if (plus_co56) then ! co56 jacobian elements dfdy(ico56,ini56) = ratdum(irn56ec) @@ -2565,89 +2695,26 @@ subroutine approx21_dfdy( & ! ni56 jacobian elements - dfdy(ini56,ihe4) = y(ife52) * ratdum(irfeag) & - + y(ife52) * y(iprot) * ratdum(ir7f54) + dfdy(ini56,ihe4) = y(ife52) * ratdum(irfeag) - dfdy(ini56,ife52) = y(ihe4) * ratdum(irfeag) & - + y(ihe4)* y(iprot) * ratdum(ir7f54) - - dfdy(ini56,ife54) = y(iprot) * y(iprot) * ratdum(ir3f54) + dfdy(ini56,ife52) = y(ihe4) * ratdum(irfeag) dfdy(ini56,ini56) = -ratdum(irniga) & - - ratdum(ir4f54) & - - y(iprot) * ratdum(ir8f54) & - ratdum(irn56ec) - dfdy(ini56,iprot) = 2.0d0 * y(ife54) * y(iprot) * ratdum(ir3f54) & - + y(ife54) * y(iprot) * y(iprot) * dratdumdy1(ir3f54) & - - y(ini56) * dratdumdy1(ir4f54) & - + y(ife52) * y(ihe4)* ratdum(ir7f54) & - + y(ife52) * y(ihe4)* y(iprot) * dratdumdy1(ir7f54) & - - y(ini56) * ratdum(ir8f54) & - - y(ini56) * y(iprot) * dratdumdy1(ir8f54) - - ! photodisintegration neutrons jacobian elements - dfdy(ineut,ihe4) = 2.0d0 * ratdum(iralf1) - - dfdy(ineut,ife52) = -2.0d0 * y(ineut) * y(ineut) * ratdum(ir2f54) - - dfdy(ineut,ife54) = 2.0d0 * ratdum(ir1f54) & - - 2.0d0 * y(ineut) * y(ineut) * ratdum(irfe56_aux2) - - dfdy(ineut,ife56) = 2.0d0 * ratdum(irfe56_aux1) & - - fe56ec_n_neut * fe56ec_fake_factor * ratdum(irn56ec) - - dfdy(ineut,ineut) = & - 2.0d0 * y(ife54) * dratdumdy1(ir1f54) & - - 4.0d0 * y(ife52) * y(ineut) * ratdum(ir2f54) & - - 2.0d0 * y(ife52) * y(ineut) * y(ineut) * dratdumdy1(ir2f54) & - + 2.0d0 * y(ihe4) * dratdumdy1(iralf1) & - - 4.0d0 * y(ineut) * y(iprot)*y(iprot) * ratdum(iralf2) & - - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy1(iralf2) & - - ratdum(irnep) & - + 2.0d0 * y(ife56) * dratdumdy1(irfe56_aux1) & - - 4.0d0 * y(ife54) * y(ineut) * ratdum(irfe56_aux2) & - - 2.0d0 * y(ife54) * y(ineut) * y(ineut) * dratdumdy1(irfe56_aux2) - - dfdy(ineut,iprot) = 2.0d0 * y(ihe4) * dratdumdy2(iralf1) & - - 4.0d0 * y(ineut)*y(ineut) * y(iprot) * ratdum(iralf2) & - - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy2(iralf2) & - + ratdum(irpen) + dfdy(ineut,ife56) = & + -fe56ec_n_neut * fe56ec_fake_factor * ratdum(irn56ec) + dfdy(ineut,ineut) = -ratdum(irnep) + dfdy(ineut,iprot) = ratdum(irpen) ! photodisintegration protons jacobian elements - dfdy(iprot,ihe4) = 2.0d0 * y(ife52) * ratdum(ir6f54) & - + 2.0d0 * ratdum(iralf1) & - + 2.0d0 * y(ife54) * ratdum(irfe56_aux4) - - dfdy(iprot,ife52) = 2.0d0 * y(ihe4) * ratdum(ir6f54) - - dfdy(iprot,ife54) = -2.0d0 * y(iprot) * y(iprot) * ratdum(ir3f54) & - - 2.0d0 * y(iprot) * y(iprot) * ratdum(ir5f54) & - + 2.0d0 * y(ihe4) * ratdum(irfe56_aux4) - - dfdy(iprot,ife56) = -2.0d0 * y(iprot) * y(iprot) * ratdum(irfe56_aux3) - - dfdy(iprot,ini56) = 2.0d0 * ratdum(ir4f54) - - dfdy(iprot,ineut) = 2.0d0 * y(ihe4) * dratdumdy1(iralf1) & - - 4.0d0 * y(ineut) * y(iprot)*y(iprot) * ratdum(iralf2) & - - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy1(iralf2) & - + ratdum(irnep) - - dfdy(iprot,iprot) = -4.0d0 * y(ife54) * y(iprot) * ratdum(ir3f54) & - - 2.0d0 * y(ife54) * y(iprot)*y(iprot)*dratdumdy1(ir3f54) & - + 2.0d0 * y(ini56) * dratdumdy1(ir4f54) & - - 4.0d0 * y(ife54) * y(iprot) * ratdum(ir5f54) & - - 2.0d0 * y(ife54) * y(iprot)*y(iprot)*dratdumdy1(ir5f54) & - + 2.0d0 * y(ihe4) * y(ife52) * dratdumdy1(ir6f54) & - + 2.0d0 * y(ihe4) * dratdumdy2(iralf1) & - - 4.0d0 * y(ineut)*y(ineut) * y(iprot) * ratdum(iralf2) & - - 2.0d0 * y(ineut)*y(ineut) * y(iprot)*y(iprot) * dratdumdy2(iralf2) & - - ratdum(irpen) & - - 4.0d0 * y(ife56) * y(iprot) * ratdum(irfe56_aux3) & - - 2.0d0 * y(ife56) * y(iprot) * y(iprot) * dratdumdy1(irfe56_aux3) & - + 2.0d0 * y(ihe4) * y(ife54) * dratdumdy1(irfe56_aux4) + dfdy(iprot,ineut) = ratdum(irnep) + dfdy(iprot,iprot) = -ratdum(irpen) + + dfdy(1:species(plus_co56),1:species(plus_co56)) = & + dfdy(1:species(plus_co56),1:species(plus_co56)) + & + real(reduced_dfdy_qp,kind=dp) end subroutine approx21_dfdy diff --git a/net/private/net_eval.f90 b/net/private/net_eval.f90 index b499cc41e8..57477ba2f7 100644 --- a/net/private/net_eval.f90 +++ b/net/private/net_eval.f90 @@ -607,7 +607,7 @@ subroutine get_rates_with_screening(n, ierr) if (dbg) write(*,*) 'call eval_using_rate_tables' call eval_using_rate_tables( & g% num_reactions, g% reaction_id, g% rate_table, g% rattab_f1, nrattab, & - n% ye, n% logT, n% temp, n% rho, n% rate_factors, g% logttab, & + n% ye, n% logT, n% temp, n% rho, n% rate_factors, g% ttab, g% logttab, & n% rate_raw, n% rate_raw_dT, n% rate_raw_dRho, ierr) if (ierr /= 0) then if (dbg) write(*,*) 'ierr from eval_using_rate_tables' diff --git a/net/test/src/test_net_do_one.f90 b/net/test/src/test_net_do_one.f90 index 213cb76ace..42fb2e32f0 100644 --- a/net/test/src/test_net_do_one.f90 +++ b/net/test/src/test_net_do_one.f90 @@ -136,6 +136,10 @@ subroutine do1_net(handle, symbolic) call mesa_error(__FILE__, __LINE__) end if + if (net_file == 'approx21.net') then + call test_approx21_reduced_flow_derivatives(n) + end if + if (symbolic .and. .not. qt) then write (*, *) 'nonzero d_dxdt_dx entries' k = 0 @@ -180,6 +184,99 @@ subroutine do1_net(handle, symbolic) end subroutine do1_net + subroutine test_approx21_reduced_flow_derivatives(n) + use net_approx21, only: eval_approx21_reduced_flows + type(Net_Info), intent(in) :: n + + integer :: i, j, i_max, j_max + real(qp), parameter :: temp = 2.0e9_qp, rel_step = 1.0e-5_qp, & + tolerance = 1.0e-7_qp + real(qp) :: deriv, error, max_error, scale, step + real(qp), allocatable :: y(:), y_plus(:), y_minus(:) + real(qp), allocatable :: rate(:), rate_plus(:), rate_minus(:), drate(:) + real(qp), allocatable :: dydt(:), dydt_plus(:), dydt_minus(:), drate_dydt(:) + real(qp), allocatable :: dfdy(:,:) + + allocate( & + y(size(n% y)), y_plus(size(n% y)), y_minus(size(n% y)), & + rate(size(n% rate_screened)), rate_plus(size(n% rate_screened)), & + rate_minus(size(n% rate_screened)), drate(size(n% rate_screened)), & + dydt(size(n% y)), dydt_plus(size(n% y)), dydt_minus(size(n% y)), & + drate_dydt(size(n% y)), dfdy(size(n% y),size(n% y))) + + do i = 1, size(y) + y(i) = 0.2_qp + 0.01_qp*real(mod(11*i,17),kind=qp) + end do + do i = 1, size(rate) + rate(i) = 0.75_qp + 0.01_qp*real(mod(19*i,43),kind=qp) + end do + + call eval_approx21_reduced_flows(y, rate, temp, dydt, dfdy) + max_error = 0.0_qp + i_max = 0 + j_max = 0 + do j = 1, size(y) + step = rel_step*max(1.0_qp,abs(y(j))) + y_plus = y + y_minus = y + y_plus(j) = y_plus(j) + step + y_minus(j) = y_minus(j) - step + call eval_approx21_reduced_flows(y_plus, rate, temp, dydt_plus) + call eval_approx21_reduced_flows(y_minus, rate, temp, dydt_minus) + do i = 1, size(y) + deriv = (dydt_plus(i) - dydt_minus(i))/(2.0_qp*step) + scale = max(1.0_qp,abs(deriv),abs(dfdy(i,j))) + error = abs(deriv - dfdy(i,j))/scale + if (error > max_error) then + max_error = error + i_max = i + j_max = j + end if + end do + end do + if (max_error > tolerance) then + write (*, '(a,2i6,1x,es14.6)') & + 'bad approx21 reduced-flow composition derivative', i_max, j_max, max_error + call mesa_error(__FILE__, __LINE__) + end if + + max_error = 0.0_qp + i_max = 0 + j_max = 0 + do j = 1, size(rate) + drate = 0.0_qp + drate(j) = 1.0_qp + call eval_approx21_reduced_flows( & + y, rate, temp, dydt, drate=drate, reduced_drate=drate_dydt) + step = rel_step*max(1.0_qp,abs(rate(j))) + rate_plus = rate + rate_minus = rate + rate_plus(j) = rate_plus(j) + step + rate_minus(j) = rate_minus(j) - step + call eval_approx21_reduced_flows(y, rate_plus, temp, dydt_plus) + call eval_approx21_reduced_flows(y, rate_minus, temp, dydt_minus) + do i = 1, size(y) + deriv = (dydt_plus(i) - dydt_minus(i))/(2.0_qp*step) + scale = max(1.0_qp,abs(deriv),abs(drate_dydt(i))) + error = abs(deriv - drate_dydt(i))/scale + if (error > max_error) then + max_error = error + i_max = i + j_max = j + end if + end do + end do + if (max_error > tolerance) then + write (*, '(a,2i6,1x,es14.6)') & + 'bad approx21 reduced-flow rate derivative', i_max, j_max, max_error + call mesa_error(__FILE__, __LINE__) + end if + + deallocate( & + y, y_plus, y_minus, rate, rate_plus, rate_minus, drate, & + dydt, dydt_plus, dydt_minus, drate_dydt, dfdy) + end subroutine test_approx21_reduced_flow_derivatives + subroutine show_results( & g, n, logT, logRho, species, num_reactions, xin, & eps_nuc, d_eps_nuc_dRho, d_eps_nuc_dT, d_eps_nuc_dx, & diff --git a/net/test/test_output b/net/test/test_output index 6a9d138e8d..f6a752fd81 100644 --- a/net/test/test_output +++ b/net/test/test_output @@ -194,19 +194,19 @@ test_one_zone_burn_small_net number of species 21 large final abundances 1.0000000000000000D-02 - c12 1 9.8714342254068388D-01 - he4 2 1.2802186836379657D-02 + c12 1 9.8714342211786166D-01 + he4 2 1.2802187260317500D-02 - xsum 1.0000000000006330D+00 + xsum 1.0000000000006328D+00 test_one_zone_burn_const_P number of species 21 large final abundances 1.0000000000000000D-02 - o16 1 6.9559810129898159D-01 - si28 2 2.0468825030127175D-01 - s32 3 4.4368437194011109D-02 - mg24 4 3.7858761619376645D-02 + o16 1 6.9559810129898514D-01 + si28 2 2.0468825030127566D-01 + s32 3 4.4368437194014287D-02 + mg24 4 3.7858761619371552D-02 xsum 1.0000000000000000D+00 diff --git a/rates/private/rates_support.f90 b/rates/private/rates_support.f90 index 57703b6812..12250b2e7d 100644 --- a/rates/private/rates_support.f90 +++ b/rates/private/rates_support.f90 @@ -18,7 +18,7 @@ ! *********************************************************************** module rates_support - use const_def, only: dp, use_mesa_temp_cache, missing_value, ln10 + use const_def, only: dp, qp, use_mesa_temp_cache, missing_value, ln10 use math_lib use rates_def use utils_lib, only: mv, switch_str, mesa_error @@ -32,12 +32,12 @@ module rates_support subroutine do_get_raw_rates( & num_reactions, reaction_id, rattab, rattab_f1, nT8s, & - ye, logtemp_in, btemp, bden, raw_rate_factor, logttab, & + ye, logtemp_in, btemp, bden, raw_rate_factor, ttab, logttab, & rate_raw, rate_raw_dT, rate_raw_dRho, ierr) integer, intent(in) :: num_reactions, reaction_id(:), nT8s real(dp), intent(in) :: & ye, logtemp_in, btemp, bden, raw_rate_factor(:), & - rattab(:,:), logttab(:) + rattab(:,:), ttab(:), logttab(:) real(dp), pointer, intent(in) :: rattab_f1(:) real(dp), intent(inout), dimension(:) :: rate_raw, rate_raw_dT, rate_raw_dRho integer, intent(out) :: ierr @@ -155,7 +155,13 @@ subroutine get_rates_from_table(r1, r2) do while (logtemp > logttab(k+1) .and. k+1 < nrattab) k = k+1 end do - dt = logtemp - logttab(k) + if (logtemp_in >= max_safe_logT_for_rates) then + dt = logtemp - logttab(k) + else + ! Preserve local changes lost by subtracting absolute logarithms. + dt = real(log(real(btemp,kind=qp)/real(ttab(k),kind=qp))/ & + real(ln10,kind=qp),kind=dp) + end if do i = r1,r2 diff --git a/rates/public/rates_lib.f90 b/rates/public/rates_lib.f90 index 6d3378ed6d..4de9078233 100644 --- a/rates/public/rates_lib.f90 +++ b/rates/public/rates_lib.f90 @@ -769,18 +769,18 @@ end subroutine eval_weak_reaction_info subroutine eval_using_rate_tables( & num_reactions, reaction_id, rattab, rattab_f1, nT8s, & - ye, logtemp, btemp, bden, raw_rate_factor, logttab, & + ye, logtemp, btemp, bden, raw_rate_factor, ttab, logttab, & rate_raw, rate_raw_dT, rate_raw_dRho, ierr) use rates_support, only : do_get_raw_rates integer, intent(in) :: num_reactions, reaction_id(:), nT8s real(dp), intent(in) :: & ye, logtemp, btemp, bden, raw_rate_factor(:), & - rattab(:,:), logttab(:) + rattab(:,:), ttab(:), logttab(:) real(dp), pointer :: rattab_f1(:) real(dp), intent(out), dimension(:) :: rate_raw, rate_raw_dT, rate_raw_dRho integer, intent(out) :: ierr call do_get_raw_rates(num_reactions, reaction_id, rattab, rattab_f1, nT8s, & - ye, logtemp, btemp, bden, raw_rate_factor, logttab, & + ye, logtemp, btemp, bden, raw_rate_factor, ttab, logttab, & rate_raw, rate_raw_dT, rate_raw_dRho, ierr) end subroutine eval_using_rate_tables diff --git a/star/defaults/controls.defaults b/star/defaults/controls.defaults index 3f15424e37..22514a41d9 100644 --- a/star/defaults/controls.defaults +++ b/star/defaults/controls.defaults @@ -9149,6 +9149,8 @@ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ignore_species_in_max_correction ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ! solver_use_T_reference + ! ~~~~~~~~~~~~~~~~~~~~~~ ! num_times_solver_reuse_mtx ! ~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -9165,6 +9167,13 @@ tiny_corr_factor = 2 ignore_min_corr_coeff_for_scale_max_correction = .false. ignore_species_in_max_correction = .false. + + ! During solver iterations, calculate the trial temperature from the + ! temperature at the start of the solver call and the logarithmic + ! temperature correction. This avoids losing small corrections when + ! adding them to the absolute value of lnT. + + solver_use_T_reference = .true. num_times_solver_reuse_mtx = 0 diff --git a/star/defaults/controls_dev.defaults b/star/defaults/controls_dev.defaults index fd96bba1f5..2168283919 100644 --- a/star/defaults/controls_dev.defaults +++ b/star/defaults/controls_dev.defaults @@ -24,6 +24,16 @@ use_TDC_Y_face_seeded_newton = .false. + ! solver_damping_on_retry + ! ~~~~~~~~~~~~~~~~~~~~~~~ + ! If true, enable the global correction coefficient limit immediately + ! on retries. If false, retries outside relaxation use the same correction + ! coefficients as the original timestep. Relaxation retries remain damped. + ! + ! :: + + solver_damping_on_retry = .true. + ! Hydro matrix solver options ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! hydro_matrix_solver selects the mesa-star hydro matrix solve. diff --git a/star/private/alloc.f90 b/star/private/alloc.f90 index 675345b21a..bdcf779255 100644 --- a/star/private/alloc.f90 +++ b/star/private/alloc.f90 @@ -1188,6 +1188,9 @@ subroutine star_info_arrays(s, c_in, action_in, ierr) call do2(s% xa_sub_xa_start, c% xa_sub_xa_start, species, 'xa_sub_xa_start') if (failed('xa_sub_xa_start')) exit + call do1(s% solver_T_reference, c% solver_T_reference) + if (failed('solver_T_reference')) exit + call do1(s% lnd_start, c% lnd_start) if (failed('lnd_start')) exit call do1(s% lnPgas_start, c% lnPgas_start) diff --git a/star/private/ctrls_io.f90 b/star/private/ctrls_io.f90 index 7b6cbff920..c520f7c141 100644 --- a/star/private/ctrls_io.f90 +++ b/star/private/ctrls_io.f90 @@ -361,7 +361,7 @@ module ctrls_io P_theta_for_velocity_time_centering, L_theta_for_velocity_time_centering, & max_logT_for_include_P_and_L_in_velocity_time_centering, & steps_before_use_TDC, use_P_d_1_div_rho_form_of_work_when_time_centering_velocity, compare_TDC_to_MLT, & - use_TDC_Y_face_seeded_newton, & + use_TDC_Y_face_seeded_newton, solver_damping_on_retry, & hydro_matrix_solver, & remesh_for_TDC_pulsations_log_core_zoning, velocity_logT_lower_bound, & max_dt_yrs_for_velocity_logT_lower_bound, velocity_tau_lower_bound, velocity_q_upper_bound, & @@ -403,6 +403,7 @@ module ctrls_io tiny_corr_coeff_limit, scale_correction_norm, corr_param_factor, num_times_solver_reuse_mtx, & scale_max_correction, ignore_min_corr_coeff_for_scale_max_correction, & ignore_too_large_correction, ignore_species_in_max_correction, & + solver_use_T_reference, & corr_norm_jump_limit, max_corr_jump_limit, resid_norm_jump_limit, max_resid_jump_limit, RSP2_use_mass_interp_face_values, & corr_coeff_limit, tiny_corr_factor, solver_test_partials_call_number, solver_test_partials_iter_number, & max_tries1, solver_max_tries_before_reject, max_tries_for_retry, max_tries_after_5_retries, solver_test_partials_sink_name, & @@ -2017,6 +2018,7 @@ subroutine store_controls(s, ierr) s% ignore_min_corr_coeff_for_scale_max_correction = ignore_min_corr_coeff_for_scale_max_correction s% ignore_too_large_correction = ignore_too_large_correction s% ignore_species_in_max_correction = ignore_species_in_max_correction + s% solver_use_T_reference = solver_use_T_reference s% corr_norm_jump_limit = corr_norm_jump_limit s% max_corr_jump_limit = max_corr_jump_limit @@ -2105,6 +2107,7 @@ subroutine store_controls(s, ierr) s% include_mlt_in_velocity_time_centering = include_mlt_in_velocity_time_centering s% compare_TDC_to_MLT = compare_TDC_to_MLT s% use_TDC_Y_face_seeded_newton = use_TDC_Y_face_seeded_newton + s% solver_damping_on_retry = solver_damping_on_retry s% hydro_matrix_solver = hydro_matrix_solver s% TDC_hydro_use_mass_interp_face_values = TDC_hydro_use_mass_interp_face_values s% TDC_hydro_nz = TDC_hydro_nz @@ -3731,6 +3734,7 @@ subroutine set_controls_for_writing(s, ierr) ignore_min_corr_coeff_for_scale_max_correction = s% ignore_min_corr_coeff_for_scale_max_correction ignore_too_large_correction = s% ignore_too_large_correction ignore_species_in_max_correction = s% ignore_species_in_max_correction + solver_use_T_reference = s% solver_use_T_reference corr_norm_jump_limit = s% corr_norm_jump_limit max_corr_jump_limit = s% max_corr_jump_limit @@ -3819,6 +3823,7 @@ subroutine set_controls_for_writing(s, ierr) include_mlt_in_velocity_time_centering = s% include_mlt_in_velocity_time_centering compare_TDC_to_MLT = s% compare_TDC_to_MLT use_TDC_Y_face_seeded_newton = s% use_TDC_Y_face_seeded_newton + solver_damping_on_retry = s% solver_damping_on_retry hydro_matrix_solver = s% hydro_matrix_solver TDC_hydro_use_mass_interp_face_values = s% TDC_hydro_use_mass_interp_face_values TDC_hydro_nz = s% TDC_hydro_nz diff --git a/star/private/evolve.f90 b/star/private/evolve.f90 index dece07d13d..fa588dd751 100644 --- a/star/private/evolve.f90 +++ b/star/private/evolve.f90 @@ -672,6 +672,8 @@ integer function do_step_part2(id, first_try) end if skip_global_corr_coeff_limit = (first_try .or. & + (.not. s% doing_relax .and. .not. s% solver_damping_on_retry .and. & + s% retry_cnt > 0) .or. & s% model_number_for_last_retry /= s% model_number) ! last alternative is for redo's s% doing_struct_burn_mix = .true. diff --git a/star/private/solver_support.f90 b/star/private/solver_support.f90 index 70d05ca43b..a07f0b1951 100644 --- a/star/private/solver_support.f90 +++ b/star/private/solver_support.f90 @@ -181,14 +181,13 @@ end subroutine dump_eval_equ end subroutine eval_equations - subroutine sizequ(s, nvar, equ_norm, equ_max, k_max, j_max, ierr) ! equ = residuals type (star_info), pointer :: s integer, intent(in) :: nvar real(dp), intent(out) :: equ_norm, equ_max integer, intent(out) :: k_max, j_max, ierr - integer :: j, k, num_terms, n, nz, nvar_hydro, nvar_chem, skip_eqn1, skip_eqn2, skip_eqn3 + integer :: j, k, num_terms, nz, skip_eqn1, skip_eqn2 real(dp) :: sumequ, absq logical :: dbg @@ -204,85 +203,34 @@ subroutine sizequ(s, nvar, equ_norm, equ_max, k_max, j_max, ierr) ! equ = resid dbg = s% solver_check_everything - nvar_hydro = min(nvar, s% nvar_hydro) - nvar_chem = s% nvar_chem - nz = s% nz - n = nz num_terms = 0 sumequ = 0 skip_eqn1 = 0 skip_eqn2 = 0 - skip_eqn3 = 0 if (s% convergence_ignore_equL_residuals) skip_eqn1 = s% i_equL if (s% convergence_ignore_alpha_RTI_residuals) skip_eqn2 = s% i_dalpha_RTI_dt - if (s% do_burn .or. s% do_mix) then - num_terms = num_terms + nvar*nz - if (skip_eqn1 > 0) num_terms = num_terms - nz - if (skip_eqn2 > 0) num_terms = num_terms - nz - if (skip_eqn3 > 0) num_terms = num_terms - nz - do k = 1, nz - do j = 1, nvar - if (j == skip_eqn1 .or. j == skip_eqn2 .or. j == skip_eqn3) cycle - if (is_bad(s% equ(j,k)) .or. is_bad(s% residual_weight(j,k))) then - ierr = 1 - return - end if - absq = abs(s% equ(j,k)*s% residual_weight(j,k)) - sumequ = sumequ + absq - if (absq > equ_max) then - equ_max = absq - j_max = j - k_max = k - end if - end do - end do - else - if (skip_eqn1 == 0 .and. skip_eqn2 == 0) then - num_terms = num_terms + nvar_hydro*nz - else if (skip_eqn1 > 0 .and. skip_eqn2 > 0) then - num_terms = num_terms + (nvar_hydro-2)*nz - else - num_terms = num_terms + (nvar_hydro-1)*nz - end if - do k = 1, nz - do j = 1, nvar_hydro - if (j == skip_eqn1 .or. j == skip_eqn2) cycle - absq = abs(s% equ(j,k)*s% residual_weight(j,k)) - sumequ = sumequ + absq - if (is_bad(sumequ)) then - if (dbg) then - write(*,3) trim(s% nameofequ(j)) // ' sumequ', j, k, sumequ - call mesa_error(__FILE__,__LINE__,'sizeq 1') - end if - ierr = -1 - if (s% report_ierr) & - write(*,3) 'bad equ(j,k)*s% residual_weight(j,k) ' // trim(s% nameofequ(j)), & - j, k, s% equ(j,k)*s% residual_weight(j,k) - if (s% stop_for_bad_nums) call mesa_error(__FILE__,__LINE__,'sizeq 2') - return - end if - if (absq > equ_max) then - equ_max = absq - j_max = j - k_max = k - end if - end do - end do - end if - if (s% do_burn .or. s% do_mix) then - num_terms = num_terms + nvar_chem*nz - do k = 1, nz - do j = nvar_hydro+1, nvar - absq = abs(s% equ(j,k)*s% residual_weight(j,k)) - sumequ = sumequ + absq - if (absq > equ_max) then - equ_max = absq - j_max = j - k_max = k - end if - end do + do k = 1, nz + do j = 1, nvar + if (j == skip_eqn1 .or. j == skip_eqn2) cycle + if (is_bad(s% equ(j,k)) .or. is_bad(s% residual_weight(j,k))) then + ierr = 1 + return + end if + num_terms = num_terms + 1 + absq = abs(s% equ(j,k)*s% residual_weight(j,k)) + sumequ = sumequ + absq + if (absq > equ_max) then + equ_max = absq + j_max = j + k_max = k + end if end do + end do + + if (num_terms == 0) then + ierr = -1 + return end if equ_norm = sumequ/num_terms @@ -1138,7 +1086,6 @@ subroutine set1(k,report,ierr) if (do_lnT) then s% lnT(k) = x(i_lnT) - s% T(k) = exp(s% lnT(k)) s% dxh_lnT(k) = s% solver_dx(i_lnT,k) if (abs(s% lnT(k) - s% lnT_start(k)) > & ln10*s% hydro_mtx_max_allowed_abs_dlogT .and. & @@ -1173,7 +1120,11 @@ subroutine set1(k,report,ierr) ierr = -1 return end if - s% T(k) = exp(s% lnT(k)) + if (s% solver_use_T_reference) then + s% T(k) = s% solver_T_reference(k)*exp(s% solver_dx(i_lnT,k)) + else + s% T(k) = exp(s% lnT(k)) + end if if (is_bad_num(s% T(k))) then s% retry_message = 'bad num for T' if (s% stop_for_bad_nums) then diff --git a/star/private/struct_burn_mix.f90 b/star/private/struct_burn_mix.f90 index 584f758ac0..417fa966f6 100644 --- a/star/private/struct_burn_mix.f90 +++ b/star/private/struct_burn_mix.f90 @@ -335,6 +335,7 @@ subroutine save_start_values(s, ierr) end if do k=1,s% nz + s% solver_T_reference(k) = s% T(k) do j=1,s% nvar_hydro s% xh_start(j,k) = s% xh(j,k) end do diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common index 86bc7b9733..62bf1b0ec2 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common +++ b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_common @@ -21,11 +21,6 @@ change_v_flag = .true. new_v_flag = .true. - ! Can be adopted to avoid issues with surface for high Z or during advanced burning - !relax_to_this_tau_factor = 1.5d6 - !dlogtau_factor = 0.1d0 - !relax_initial_tau_factor = .true. - ! If using a big network, comment out and used a local rate_tables directory num_special_rate_factors = 2 reaction_for_special_factor(1) = 'r_c12_ag_o16' @@ -49,6 +44,7 @@ &eos use_Skye = .true. use_PC = .false. + mass_fraction_limit_for_skye = 1d-16 / ! end of eos namelist &kap @@ -188,9 +184,9 @@ max_timestep_factor = 1.05d0 timestep_factor_for_retries = 0.75 - limit_for_rel_error_in_energy_conservation = 1d-7 - hard_limit_for_rel_error_in_energy_conservation = 1d-6 - + limit_for_rel_error_in_energy_conservation = 1d-2 !1d-7 + hard_limit_for_rel_error_in_energy_conservation = 1d-2 !1d-6 + never_skip_hard_limits = .true. min_xa_hard_limit = -1d-5 min_xa_hard_limit_for_highT = -3d-5 @@ -299,23 +295,17 @@ ! solver ! damped newton and structure only - scale_max_correction = 0.1d0 - ignore_species_in_max_correction = .true. - - ! might help to turn off gold2 in 8-10 Msun degenerate cores - ! for O-Ne flames. - use_gold2_tolerances = .true. - gold2_tol_max_residual2 = 5d-7 - gold2_tol_max_residual3 = 5d-4 + ! scale_max_correction = 0.1d0 + ! ignore_species_in_max_correction = .true. + solver_use_T_reference = .true. + solver_damping_on_retry = .false. + fix_d_eos_dxa_partials = .false. + use_gold2_tolerances = .false. use_gold_tolerances = .true. - gold_tol_max_residual2 = 5d-4 - !gold_tol_max_residual3 = 1d-4 ! Default in controls is 1d-5 tol_correction_high_T_limit = 1d9 ! Switch to lower tol at high temp for large Mass - solver_iters_timestep_limit = 20 gold_solver_iters_timestep_limit = 20 - iter_for_resid_tol2 = 10 max_abs_rel_run_E_err = 1d-2 @@ -325,8 +315,8 @@ convergence_ignore_equL_residuals = .true. make_gradr_sticky_in_solver_iters = .true. xa_scale = 1d-5 - iter_for_resid_tol2 = 10 - min_timestep_limit = 1d-12 ! (seconds) ! 1d-20 if things are sticky + min_timestep_limit = 1d-16 ! (seconds) ! 1d-20 if things are sticky + warn_rates_for_high_temp = .true. max_safe_logT_for_rates = 10.5d0 @@ -342,7 +332,7 @@ ! 4d9 has been well tested with the approx21 network ! lower to 2.5d9 or worst case 1d9 for large networks to help with ! numerical stability and speed, see MESA VI (Jermyn 2023) - op_split_burn = .true. + op_split_burn = .false. op_split_burn_min_T = 4d9 burn_steps_limit = 150 burn_steps_hard_limit = 250 diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_make_late_pre_zams b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_make_late_pre_zams index f8ed6c7c95..b81a1d79b6 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_make_late_pre_zams +++ b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_make_late_pre_zams @@ -24,6 +24,8 @@ &controls + solver_damping_on_retry = .false. + ! limit max_model_number as part of test_suite max_model_number = 10000 use_gold2_tolerances = .true. diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_pgstar b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_pgstar index afb609442f..0a4eaab7a9 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_pgstar +++ b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_pgstar @@ -17,7 +17,6 @@ pgstar_show_age = .false. Grid1_win_flag = .true. Grid1_win_width = 12 -Grid1_win_width = 17 Grid1_win_aspect_ratio = 0.666 Grid1_win_aspect_ratio = 0.53 diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_to_cc b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_to_cc index f4dba2a266..6ffe1b515c 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/inlist_to_cc +++ b/star/test_suite/20M_pre_ms_to_core_collapse/inlist_to_cc @@ -21,9 +21,6 @@ &controls -! to avoid a host of retries we relax our limits on split_burn -! for speed/stability in test_suite, at reduced time resolution. -op_split_burn_min_T = 1d9 ! prevent development of radial pulses during advanced burning !drag_coefficient = 1d0 @@ -56,12 +53,8 @@ ignore_too_large_correction = .true. ! for conv_vel's ! mixing ! timesteps - limit_for_rel_error_in_energy_conservation = 1d-3 - hard_limit_for_rel_error_in_energy_conservation = 1d-2 - delta_XSi_cntr_limit = 0.0025 - ! mesh ! solver diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_c_burn.mod b/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_c_burn.mod index b7eb399f4b..a793ba568c 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_c_burn.mod +++ b/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_c_burn.mod @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:27b3dbcdbeb520e582aa3aa0aa84505c8bfadc5392909415cc86fefe772c0d24 -size 771196 +oid sha256:a3b4791775e03f50defaf85166c79f5cbf5ff29698865517b7e3f3bb6528ea85 +size 772762 diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_he_burn.mod b/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_he_burn.mod index a4573cab3f..cceb3ae590 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_he_burn.mod +++ b/star/test_suite/20M_pre_ms_to_core_collapse/standard_after_core_he_burn.mod @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:85949dee2e0eb16554fa788be125bb45510dcdfbf511dd0fb8a1dd2bc4d44f00 -size 697726 +oid sha256:34acb9c120bea1367afa5cefbc14311f7008da89357e13a0d4d57c70a74ac3cd +size 662162 diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/standard_late_pre_zams.mod b/star/test_suite/20M_pre_ms_to_core_collapse/standard_late_pre_zams.mod index 16ff276464..820eff9c3b 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/standard_late_pre_zams.mod +++ b/star/test_suite/20M_pre_ms_to_core_collapse/standard_late_pre_zams.mod @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:0048c37c1cc1f072de1566c8e512a20754e8e6a336838b4b2be5a537084ce762 -size 201118 +oid sha256:17018c81c254717fcdfb62fd4a09c592086ec7396bc97b6a7a7f81b62a873144 +size 201104 diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/standard_lgTmax.mod b/star/test_suite/20M_pre_ms_to_core_collapse/standard_lgTmax.mod index e728df3c5c..dae8dfcb24 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/standard_lgTmax.mod +++ b/star/test_suite/20M_pre_ms_to_core_collapse/standard_lgTmax.mod @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:93614de8fd66541bec86d84ac636946c2354ae22b29ca40c47b9efaa35e0af33 -size 1494046 +oid sha256:efc4a9e68153a4386d40f0a9c0eff0b7b71093336cc605ba023986dadb2c0831 +size 1535902 diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/standard_removed_envelope.mod b/star/test_suite/20M_pre_ms_to_core_collapse/standard_removed_envelope.mod index 9c2fe229df..b6626959ad 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/standard_removed_envelope.mod +++ b/star/test_suite/20M_pre_ms_to_core_collapse/standard_removed_envelope.mod @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:8d26e81be912518fcb8bd8f7f5e4e27957123e4cec1eee96b06bc027851a6999 -size 697726 +oid sha256:9ae668ef056d9a4ec8f78082e79aa69c932fac39a0475ca0dedc6c51e47a3a4e +size 662162 diff --git a/star/test_suite/20M_pre_ms_to_core_collapse/standard_zams.mod b/star/test_suite/20M_pre_ms_to_core_collapse/standard_zams.mod index c0d7a472fc..5768b6dabe 100644 --- a/star/test_suite/20M_pre_ms_to_core_collapse/standard_zams.mod +++ b/star/test_suite/20M_pre_ms_to_core_collapse/standard_zams.mod @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:d3b979b176822035c37f0b8f93e2eecb8e38eba6310945da57d845ec27cc3149 -size 312206 +oid sha256:8bba1ae15aa87ab3d64d68a4d123a79a8053918a48f140bc924785e7ea0d90bb +size 312982 diff --git a/star_data/private/star_controls.inc b/star_data/private/star_controls.inc index b4e3b4174a..9265b21985 100644 --- a/star_data/private/star_controls.inc +++ b/star_data/private/star_controls.inc @@ -985,6 +985,7 @@ logical :: ignore_min_corr_coeff_for_scale_max_correction logical :: ignore_too_large_correction logical :: ignore_species_in_max_correction + logical :: solver_use_T_reference real(dp) :: corr_norm_jump_limit real(dp) :: max_corr_jump_limit real(dp) :: resid_norm_jump_limit diff --git a/star_data/private/star_controls_dev.inc b/star_data/private/star_controls_dev.inc index 336bad4f9e..f887d8a3bd 100644 --- a/star_data/private/star_controls_dev.inc +++ b/star_data/private/star_controls_dev.inc @@ -1,6 +1,7 @@ logical :: compare_TDC_to_MLT logical :: use_TDC_Y_face_seeded_newton + logical :: solver_damping_on_retry character(32) :: hydro_matrix_solver logical :: TDC_use_density_form_for_eddy_viscosity diff --git a/star_data/public/star_data_step_work.inc b/star_data/public/star_data_step_work.inc index 1991b965df..a378d49447 100644 --- a/star_data/public/star_data_step_work.inc +++ b/star_data/public/star_data_step_work.inc @@ -809,6 +809,7 @@ real(dp), pointer :: xh_start(:,:) ! (nvar_hydro,nz) real(dp), pointer :: xa_start(:,:) ! (species,nz) real(dp), pointer :: xa_sub_xa_start(:,:) ! (species,nz) + real(dp), pointer :: solver_T_reference(:) ! (nz) real(dp), pointer :: lnd_start(:) ! (nz) real(dp), pointer :: Peos_start(:) ! (nz) real(dp), pointer :: Peos_face_start(:) ! (nz)