diff --git a/source/xnet_conditions.f90 b/source/xnet_conditions.f90 index 91f2686d..289ea721 100644 --- a/source/xnet_conditions.f90 +++ b/source/xnet_conditions.f90 @@ -132,35 +132,8 @@ Subroutine t9rhofind_vector(kstep,tf,nf,t9f,rhof,mask_in) Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - - ! For constant conditions (nh = 1), set temperature and density - If ( nh(izb) == 1 ) Then - t9f(izb) = t9h(1,izb) - rhof(izb) = rhoh(1,izb) - nf(izb) = 1 - - ! Otherwise, calculate T9 and rho by interpolation - Else - Do n = 1, nh(izb) - If ( tf(izb) <= th(n,izb) ) Exit - EndDo - nf(izb) = n - If ( n > 1 .and. n <= nh(izb) ) Then - rdt = 1.0 / (th(n,izb)-th(n-1,izb)) - dt = tf(izb) - th(n-1,izb) - dt9 = t9h(n,izb) - t9h(n-1,izb) - drho = rhoh(n,izb) - rhoh(n-1,izb) - t9f(izb) = dt*rdt*dt9 + t9h(n-1,izb) - rhof(izb) = dt*rdt*drho + rhoh(n-1,izb) - ElseIf ( n == 1 ) Then - t9f(izb) = t9h(1,izb) - rhof(izb) = rhoh(1,izb) - Else - t9f(izb) = t9h(nh(izb),izb) - rhof(izb) = rhoh(nh(izb),izb) - Write(lun_stdout,*) 'Time beyond thermodynamic range',tf(izb),' >',th(nh(izb),izb) - EndIf - EndIf + Call t9rhofind_scalar(kstep,tf(izb),nf(izb),t9f(izb),rhof(izb), & + & nh(izb),th(:,izb),t9h(:,izb),rhoh(:,izb)) EndIf EndDo diff --git a/source/xnet_controls.f90 b/source/xnet_controls.f90 index 9a33e44b..cd7ef042 100644 --- a/source/xnet_controls.f90 +++ b/source/xnet_controls.f90 @@ -205,7 +205,7 @@ Subroutine read_controls(data_dir) Allocate (zone_id(3,nzevolve)) Allocate (lzactive(nzevolve)) Allocate (iweak(nzevolve),lun_ev(nzevolve),lun_ts(nzevolve)) - Allocate (kmon(2,nzevolve),ktot(5,nzevolve)) + Allocate (kmon(5,nzevolve),ktot(5,nzevolve)) !$omp parallel default(shared) zb_offset = (tid-1) * nzbatchmx zb_lo = zb_offset + 1 diff --git a/source/xnet_data.f90 b/source/xnet_data.f90 index b51b188b..5dad49ee 100644 --- a/source/xnet_data.f90 +++ b/source/xnet_data.f90 @@ -108,19 +108,19 @@ Subroutine partf(t9,mask_in) !----------------------------------------------------------------------------------------------- ! This routine calculates the nuclear partition functions as a function of temperature. !----------------------------------------------------------------------------------------------- - Use xnet_controls, Only: idiag, iheat, lun_diag, nzevolve, zb_lo, zb_hi, lzactive + Use xnet_controls, Only: idiag, iheat, lun_diag, zb_lo, zb_hi, lzactive Use xnet_types, Only: dp Use xnet_util, Only: safe_exp Implicit None ! Input variables - Real(dp), Intent(in) :: t9(nzevolve) + Real(dp), Intent(in) :: t9(zb_lo:zb_hi) ! Optional variables Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Integer :: i, ii, izb + Integer :: i, ii, k, izb Real(dp) :: rdt9 Logical, Pointer :: mask(:) @@ -139,27 +139,39 @@ Subroutine partf(t9,mask_in) ii = i ! Linear interpolation in log-space + gg(0,izb) = 1.0 ! placeholder for non-nuclei, gamma-rays, etc. Select Case (ii) Case (1) - gg(1:ny,izb) = g(1,1:ny) + Do k = 1, ny + gg(k,izb) = g(1,k) + EndDo Case (ng+1) - gg(1:ny,izb) = g(ng,1:ny) + Do k = 1, ny + gg(k,izb) = g(ng,k) + EndDo Case Default rdt9 = (t9(izb)-t9i(ii-1)) / (t9i(ii)-t9i(ii-1)) - gg(1:ny,izb) = safe_exp( rdt9*log(g(ii,1:ny)) + (1.0-rdt9)*log(g(ii-1,1:ny)) ) + Do k = 1, ny + gg(k,izb) = safe_exp( rdt9*log(g(ii,k)) + (1.0-rdt9)*log(g(ii-1,k)) ) + EndDo End Select - gg(0,izb) = 1.0 ! placeholder for non-nuclei, gamma-rays, etc. If ( iheat > 0 ) Then + dlngdt9(0,izb) = 0.0 Select Case (ii) Case (1) - dlngdt9(1:ny,izb) = log(g(2,1:ny)/g(1,1:ny)) / (t9i(2)-t9i(1)) + Do k = 1, ny + dlngdt9(k,izb) = log(g(2,k)/g(1,k)) / (t9i(2)-t9i(1)) + EndDo Case (ng+1) - dlngdt9(1:ny,izb) = log(g(ng,1:ny)/g(ng-1,1:ny)) / (t9i(ng)-t9i(ng-1)) + Do k = 1, ny + dlngdt9(k,izb) = log(g(ng,k)/g(ng-1,k)) / (t9i(ng)-t9i(ng-1)) + EndDo Case Default - dlngdt9(1:ny,izb) = log(g(ii,1:ny)/g(ii-1,1:ny)) / (t9i(ii)-t9i(ii-1)) + Do k = 1, ny + dlngdt9(k,izb) = log(g(ii,k)/g(ii-1,k)) / (t9i(ii)-t9i(ii-1)) + EndDo End Select - dlngdt9(0,izb) = 0.0 EndIf EndIf EndDo diff --git a/source/xnet_evolve.f90 b/source/xnet_evolve.f90 index 79785b94..0280a2f6 100644 --- a/source/xnet_evolve.f90 +++ b/source/xnet_evolve.f90 @@ -54,12 +54,6 @@ Subroutine full_net(kstep) start_timer = xnet_wtime() timer_xnet = timer_xnet - start_timer - ! Initialize counters - kstep = 0 - Do izb = zb_lo, zb_hi - kmon(:,izb) = 0 - ktot(:,izb) = 0 - EndDo ! Set reaction controls not read in from control idiag0 = idiag @@ -76,10 +70,16 @@ Subroutine full_net(kstep) t9t(izb) = t9(izb) rhoo(izb) = rho(izb) rhot(izb) = rho(izb) - yet(izb) = ye(izb) yeo(izb) = ye(izb) - yo(:,izb) = y(:,izb) - yt(:,izb) = y(:,izb) + yet(izb) = ye(izb) + Do k = 1, ny + yo(k,izb) = y(k,izb) + yt(k,izb) = y(k,izb) + EndDo + Do k = 1, 5 + kmon(k,izb) = 0 + ktot(k,izb) = 0 + EndDo EndDo en0 = 0.0 @@ -100,16 +100,20 @@ Subroutine full_net(kstep) ! Output initial abundances and conditions Call ts_output(0,delta_en,edot) - ! Start evolution + ! Initialize timestep loop flags + kstep = 0 + mykstep = 0 + lzsolve = lzactive(zb_lo:zb_hi) + lzoutput = lzsolve Do izb = zb_lo, zb_hi - If ( lzactive(izb) ) Then + If ( lzsolve(izb) ) Then its(izb) = 0 Else its(izb) = -1 EndIf - lzsolve(izb) = lzactive(izb) - mykstep(izb) = 0 EndDo + + ! Start evolution Do kstep = 1, kstmx ! Determine if this is an output step @@ -129,47 +133,55 @@ Subroutine full_net(kstep) Call solve_be(kstep,its) End Select - Do izb = zb_lo, zb_hi - izone = izb + szbatch - zb_lo + ! If convergence is successful, output timestep results + If ( idiag >= 1 ) Then + Do izb = zb_lo, zb_hi + izone = izb + szbatch - zb_lo + If ( its(izb) == 0 .and. idiag >= 1 ) Then + Write(lun_diag,"(2(a,i5),5es14.7)") & + & 'KStep ',kstep,' Zone ',izone,t(izb),tdel(izb),t9(izb),rho(izb),ye(izb) + If ( idiag >= 3 ) Then + Write(lun_diag,"(a8)") 'delta Y' + Write(lun_diag,"(8x,a5,4es12.4)") & + & (nname(k),y(k,izb),yo(k,izb),(y(k,izb)-yo(k,izb)),(tdel(izb)*ydot(k,izb)),k=1,ny) + If ( iheat > 0 ) Write(lun_diag,"(8x,a5,4es12.4)") & + & 'T9',t9(izb),t9o(izb),t9(izb)-t9o(izb),tdel(izb)*t9dot(izb) + EndIf + ElseIf ( its(izb) == 1 ) Then + Write(lun_diag,"(a,2(i5,a))") 'KStep ',kstep,' Zone ',izone,' Inter!!!' + EndIf + EndDo + EndIf - ! If convergence is successful, output timestep results + Do izb = zb_lo, zb_hi If ( its(izb) == 0 ) Then - If ( idiag >= 1 ) Write(lun_diag,"(2(a,i5),5es14.7)") & - & 'KStep ',kstep,' Zone ',izone,t(izb),tdel(izb),t9(izb),rho(izb),ye(izb) - If ( idiag >= 3 ) Then - Write(lun_diag,"(a8)") 'delta Y' - Write(lun_diag,"(8x,a5,4es12.4)") & - & (nname(k),y(k,izb),yo(k,izb),(y(k,izb)-yo(k,izb)),(tdel(izb)*ydot(k,izb)),k=1,ny) - If ( iheat > 0 ) Write(lun_diag,"(8x,a5,4es12.4)") & - & 'T9',t9(izb),t9o(izb),t9(izb)-t9o(izb),tdel(izb)*t9dot(izb) - EndIf + enold(izb) = enm(izb) Call benuc(yt(:,izb),enb(izb),enm(izb)) delta_en(izb) = enm(izb) - en0(izb) edot(izb) = -(enm(izb)-enold(izb)) / tdel(izb) - ! If reduced timesteps fail to successfully integrate, warn and flag to remove from loop - ElseIf ( its(izb) == 1 ) Then - If ( idiag >= 0 ) Write(lun_diag,"(a,2(i5,a))") 'KStep ',kstep,' Zone ',izone,' Inter!!!' - its(izb) = 2 - EndIf - EndDo - - lzoutput = ( lzsolve .and. its <= 0 ) - Call ts_output(kstep,delta_en,edot,mask_in = lzoutput) - - Do izb = zb_lo, zb_hi - izone = izb + szbatch - zb_lo - ! If this zone reaches the stop time, flag it to remove from loop If ( t(izb) >= tstop(izb) ) Then mykstep(izb) = kstep its(izb) = -1 + lzsolve(izb) = .false. EndIf + + ! If reduced timesteps fail to successfully integrate, flag it to remove from loop + ElseIf ( its(izb) == 1 ) Then + its(izb) = 2 + lzsolve(izb) = .false. + lzoutput(izb) = .false. + ElseIf ( its(izb) < 0 ) Then + lzsolve(izb) = .false. + lzoutput(izb) = .false. + EndIf EndDo + Call ts_output(kstep,delta_en,edot,mask_in = lzoutput) + ! Test if all zones have stopped - lzsolve = ( its == 0 ) If ( .not. any( lzsolve ) ) Exit EndDo diff --git a/source/xnet_ffn.f90 b/source/xnet_ffn.f90 index bae33e2d..722a5b87 100644 --- a/source/xnet_ffn.f90 +++ b/source/xnet_ffn.f90 @@ -20,11 +20,14 @@ Module xnet_ffn Real(dp), Allocatable :: ffn_ec(:,:), ffn_beta(:,:) ! dim(nffn,ngrid) Real(dp), Allocatable :: ffnsum(:,:), ffnenu(:,:) ! dim(nffn,ngrid) Real(dp), Allocatable :: ffn_qval(:) + Real(dp), Allocatable :: phasei(:,:), dphaseidt9(:,:) ! dim(nffn,ngrid) Integer, Allocatable :: has_logft(:) Real(dp), Allocatable :: rffn(:,:) ! FFN reaction rates Real(dp), Allocatable :: dlnrffndt9(:,:) ! log FFN reaction rates + Real(dp), Parameter :: lrfmin = -30.0, rfmin = 1.0e-30 + Contains Subroutine read_ffn_data(nffn,data_dir) @@ -80,7 +83,7 @@ Subroutine read_ffn_data_logft(nffn,data_dir) !----------------------------------------------------------------------------------------------- ! This routine allocates and loads the data structures for FFN reaction rates. !----------------------------------------------------------------------------------------------- - Use xnet_controls, Only: lun_diag + Use xnet_controls, Only: lun_diag, nzevolve Implicit None ! Input variables @@ -96,11 +99,14 @@ Subroutine read_ffn_data_logft(nffn,data_dir) Allocate (ffnsum(nffn,ngrid),ffnenu(nffn,ngrid)) Allocate (ffn_ec(nffn,ngrid),ffn_beta(nffn,ngrid)) Allocate (ffn_qval(nffn)) + Allocate (phasei(nffn,nzevolve),dphaseidt9(nffn,nzevolve)) ffnsum = 0.0 ffnenu = 0.0 ffn_ec = 0.0 ffn_beta = 0.0 ffn_qval = 0.0 + phasei = 0.0 + dphaseidt9 = 0.0 Open(newunit=lun_ffn, file=trim(data_dir)//"/netweak", status='old') Do i = 1, nffn @@ -119,9 +125,11 @@ Subroutine read_ffn_data_logft(nffn,data_dir) EndDo Close(lun_ffn) - Where (ffn_ec < -30d0) - ffn_ec = 99.d0 - End Where + Do i = 1, nffn + Do j = 1, ngrid + If ( ffn_ec(i,j) < lrfmin ) ffn_ec(i,j) = 99.0 + EndDo + EndDo Return End Subroutine read_ffn_data_logft @@ -131,7 +139,7 @@ Subroutine ffn_rate(nffn,t9,ene,rf,dlnrfdt9,mask_in) ! This routine calculates the reaction rates for FFN weak rates !----------------------------------------------------------------------------------------------- Use xnet_constants, Only: ln_2, ln_10, bok, m_e - Use xnet_controls, Only: iheat, nzevolve, zb_lo, zb_hi, lzactive + Use xnet_controls, Only: iheat, zb_lo, zb_hi, lzactive Use xnet_conditions, Only: etae Use xnet_types, Only: dp Use xnet_controls, Only: lun_diag, idiag @@ -140,7 +148,7 @@ Subroutine ffn_rate(nffn,t9,ene,rf,dlnrfdt9,mask_in) ! Input variables Integer, Intent(in) :: nffn ! Number of FFN rates - Real(dp), Intent(in) :: t9(nzevolve) ! Temperature [GK] + Real(dp), Intent(in) :: t9(zb_lo:zb_hi) ! Temperature [GK] Real(dp), Intent(in) :: ene(zb_lo:zb_hi) ! Electron Density [g cm^{-3}] ! Output variables @@ -151,11 +159,10 @@ Subroutine ffn_rate(nffn,t9,ene,rf,dlnrfdt9,mask_in) Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Real(dp), Parameter :: lrfmin = -30.0, rfmin=1d-30 Real(dp) :: r1, r2, dr1, dr2, dr1_ec, dr2_ec - Real(dp) :: rf_beta,rf_ec, phasei, cheme, dphase_dt - Real(dp) :: enel, dt9, dene, rdt9, rdene, drbeta_dt,drec_dt - Integer :: i,izb, k, le1, lt1, i1, i2, i3, i4 + Real(dp) :: rf_beta, rf_ec, cheme + Real(dp) :: enel, dt9, dene, rdt9, rdene, drbeta_dt, drec_dt + Integer :: i, izb, k, le1, lt1, i1, i2, i3, i4 Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -165,6 +172,20 @@ Subroutine ffn_rate(nffn,t9,ene,rf,dlnrfdt9,mask_in) EndIf If ( .not. any(mask) ) Return + ! Pre-calculate phase space integrals and derivatives + Do izb = zb_lo, zb_hi + Do k = 1, nffn + If ( mask(izb) .and. has_logft(k) > 0 ) Then + cheme = etae(izb)*bok*t9(izb) + m_e + If ( has_logft(k) == 1 ) Then + Call effphase(t9(izb),+cheme,ffn_qval(k),phasei(k,izb),dphaseidt9(k,izb)) + ElseIf ( has_logft(k) == 2 ) Then + Call effphase(t9(izb),-cheme,ffn_qval(k),phasei(k,izb),dphaseidt9(k,izb)) + EndIf + EndIf + EndDo + EndDo + Do izb = zb_lo, zb_hi If ( mask(izb) ) Then @@ -188,84 +209,77 @@ Subroutine ffn_rate(nffn,t9,ene,rf,dlnrfdt9,mask_in) i3 = nt9grid*le1 + lt1 i4 = i3 + 1 - cheme = etae(izb)*bok*t9(izb)+m_e Do k = 1, nffn - - If (has_logft(k) > 0) then - - ! Calculate phase space integral and derivative - If ( has_logft(k) == 1 ) Then - Call effphase(t9(izb),cheme,ffn_qval(k),phasei,dphase_dt) - ElseIf ( has_logft(k) == 2 ) Then - Call effphase(t9(izb),-cheme,ffn_qval(k),phasei,dphase_dt) - Else - ! TODO: move this out of the loop - Write(*,*) "Problem in logft option. Should be 1 or 2, but is ",has_logft(k) - Stop - EndIf - - dr1_ec = ffn_ec(k,i2) - ffn_ec(k,i1) - dr2_ec = ffn_ec(k,i4) - ffn_ec(k,i3) - r1 = ffn_ec(k,i1) + rdt9*dr1_ec - r2 = ffn_ec(k,i3) + rdt9*dr2_ec - rf_ec = r1 + rdene*(r2 - r1) ! logft - If ( phasei > rfmin .and. rf_ec > lrfmin ) then - rf_ec = ln_2 * phasei / 10.0**rf_ec ! turn into rate - Else - rf_ec = 0.0 - EndIf - - dr1 = ffn_beta(k,i2) - ffn_beta(k,i1) - dr2 = ffn_beta(k,i4) - ffn_beta(k,i3) - r1 = ffn_beta(k,i1) + rdt9*dr1 - r2 = ffn_beta(k,i3) + rdt9*dr2 - rf_beta = r1 + rdene*(r2 - r1) - If ( rf_beta > lrfmin ) Then - rf_beta = 10.0**rf_beta - Else - rf_beta = 0.0 - EndIf - - rf(k,izb) = rf_beta + rf_ec - - If ( rf_beta < rfmin .and. rf_ec < rfmin ) Then - rf(k,izb) = 0.0 - dlnrfdt9(k,izb) = 0.0 - Else - ! Temperature derivative - If (rf_ec < rfmin) Then - dlnrfdt9(k,izb) = ln_10 * ( rdene*dr2 + (1.0-rdene)*dr1 ) / dt9 - Else - drbeta_dt = rf_beta*ln_10 * ( rdene*dr2 + (1.0-rdene)*dr1 ) / dt9 - drec_dt = rf_ec * ( -ln_10 * ( rdene*dr2_ec + (1.0-rdene)*dr1_ec ) / dt9 & - & + dphase_dt / phasei ) - dlnrfdt9(k,izb) = (drbeta_dt+drec_dt)/rf(k,izb) - EndIf - - If ( idiag >= 5 ) Then - Write(lun_diag,*) "FFN rate, deriv., phase space: ", & - & k,rf(k,izb),dlnrfdt9(k,izb),phasei - EndIf - - EndIf + If ( has_logft(k) > 0 ) Then + + dr1_ec = ffn_ec(k,i2) - ffn_ec(k,i1) + dr2_ec = ffn_ec(k,i4) - ffn_ec(k,i3) + r1 = ffn_ec(k,i1) + rdt9*dr1_ec + r2 = ffn_ec(k,i3) + rdt9*dr2_ec + rf_ec = r1 + rdene*(r2 - r1) ! logft + If ( phasei(k,izb) > rfmin .and. rf_ec > lrfmin ) then + rf_ec = ln_2 * phasei(k,izb) / 10.0**rf_ec ! turn into rate + Else + rf_ec = 0.0 + EndIf + + dr1 = ffn_beta(k,i2) - ffn_beta(k,i1) + dr2 = ffn_beta(k,i4) - ffn_beta(k,i3) + r1 = ffn_beta(k,i1) + rdt9*dr1 + r2 = ffn_beta(k,i3) + rdt9*dr2 + rf_beta = r1 + rdene*(r2 - r1) + If ( rf_beta > lrfmin ) Then + rf_beta = 10.0**rf_beta + Else + rf_beta = 0.0 + EndIf + + rf(k,izb) = rf_beta + rf_ec + + If ( rf_beta < rfmin .and. rf_ec < rfmin ) Then + rf(k,izb) = 0.0 + dlnrfdt9(k,izb) = 0.0 + Else + If (rf_ec < rfmin) Then + dlnrfdt9(k,izb) = ln_10 * ( rdene*dr2 + (1.0-rdene)*dr1 ) / dt9 + Else + drbeta_dt = rf_beta*ln_10 * ( rdene*dr2 + (1.0-rdene)*dr1 ) / dt9 + drec_dt = rf_ec * ( -ln_10 * ( rdene*dr2_ec + (1.0-rdene)*dr1_ec ) / dt9 & + & + dphaseidt9(k,izb) / phasei(k,izb) ) + dlnrfdt9(k,izb) = ( drbeta_dt + drec_dt ) / rf(k,izb) + EndIf + EndIf Else - dr1 = ffnsum(k,i2) - ffnsum(k,i1) - dr2 = ffnsum(k,i4) - ffnsum(k,i3) - r1 = ffnsum(k,i1) + rdt9*dr1 - r2 = ffnsum(k,i3) + rdt9*dr2 - rf(k,izb) = r1 + rdene*(r2 - r1) - If ( rf(k,izb) < lrfmin ) Then - rf(k,izb) = 0.0 - dlnrfdt9(k,izb) = 0.0 - Else - rf(k,izb) = 10.0**rf(k,izb) - dlnrfdt9(k,izb) = ln_10 * ( rdene*dr2 + (1.0-rdene)*dr1 ) / dt9 - EndIf + dr1 = ffnsum(k,i2) - ffnsum(k,i1) + dr2 = ffnsum(k,i4) - ffnsum(k,i3) + r1 = ffnsum(k,i1) + rdt9*dr1 + r2 = ffnsum(k,i3) + rdt9*dr2 + rf(k,izb) = r1 + rdene*(r2 - r1) + If ( rf(k,izb) < lrfmin ) Then + rf(k,izb) = 0.0 + dlnrfdt9(k,izb) = 0.0 + Else + rf(k,izb) = 10.0**rf(k,izb) + dlnrfdt9(k,izb) = ln_10 * ( rdene*dr2 + (1.0-rdene)*dr1 ) / dt9 + EndIf EndIf EndDo EndIf EndDo + If ( idiag >= 5 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Do k = 1, nffn + If ( has_logft(k) > 0 ) Then + Write(lun_diag,"(a,i5,4es23.15)") "FFN rate, deriv., phase space: ", & + & k,rf(k,izb),dlnrfdt9(k,izb),phasei(k,izb),dphaseidt9(k,izb) + EndIf + EndDo + EndIf + EndDo + EndIf + Return End Subroutine ffn_rate diff --git a/source/xnet_integrate.f90 b/source/xnet_integrate.f90 index fc2f25a0..daf79ed4 100644 --- a/source/xnet_integrate.f90 +++ b/source/xnet_integrate.f90 @@ -21,11 +21,12 @@ Subroutine timestep(kstep,mask_in) !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: ny, nname Use xnet_abundances, Only: y, yo, yt, ydot - Use xnet_conditions, Only: t, tt, tdel, tdel_next, tdel_old, t9, t9o, t9t, rho, rhot, & - & t9dot, cv, nt, ntt, ints, intso, tstop, tdelstart, nh, th, t9h, rhoh, t9rhofind - Use xnet_controls, Only: changemx, changemxt, idiag, iheat, iweak, lun_diag, yacc, & - & szbatch, zb_lo, zb_hi, lzactive, iscrn, ymin + Use xnet_conditions, Only: t, tt, tdel, tdel_next, tdel_old, t9, t9o, t9t, rho, & + & rhot, t9dot, cv, nt, ntt, ints, intso, tstop, tdelstart, nh, th, t9h, rhoh, t9rhofind + Use xnet_controls, Only: changemx, changemxt, idiag, iheat, iscrn, iweak, lun_diag, yacc, & + & ymin, szbatch, zb_lo, zb_hi, lzactive Use xnet_types, Only: dp + Use xnet_util, Only: xnet_terminate Implicit None ! Input variables @@ -35,11 +36,14 @@ Subroutine timestep(kstep,mask_in) Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Integer :: i, j, izb, izone - Real(dp) :: ydotoy(ny), t9dotot9 - Real(dp) :: changest, changeth, dtherm(zb_lo:zb_hi) - Real(dp) :: tdel_stop, tdel_deriv, tdel_fine, tdel_t9, tdel_init - Logical :: mask_deriv(zb_lo:zb_hi), mask_profile(zb_lo:zb_hi) + Real(dp), Parameter :: changeth = 0.1 + Real(dp) :: changest, yfloor + Real(dp) :: rtau_y(ny), changey, tdel_dy(zb_lo:zb_hi) + Real(dp) :: rtau_t9, changet9, tdel_dt9(zb_lo:zb_hi) + Real(dp) :: tdel_stop(zb_lo:zb_hi), tdel_init + Real(dp) :: dtherm(zb_lo:zb_hi) + Integer :: i, j, k, izb, izone + Logical :: mask_init(zb_lo:zb_hi) Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -50,188 +54,184 @@ Subroutine timestep(kstep,mask_in) If ( .not. any(mask) ) Return ! Retain old values of timestep and thermo and calculate remaining time - changeth = 0.1 Do izb = zb_lo, zb_hi If ( mask(izb) ) Then tdel_old(izb) = tdel(izb) + mask_init(izb) = ( abs(tdel_old(izb)) < tiny(0.0) ) + Else + mask_init(izb) = .false. EndIf - mask_deriv(izb) = ( mask(izb) .and. abs(tdel_old(izb)) < tiny(0.0) ) EndDo - Call cross_sect(mask_in = mask_deriv) - Call yderiv(mask_in = mask_deriv) + Call cross_sect(mask_in = mask_init) + Call yderiv(mask_in = mask_init) Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - izone = izb + szbatch - zb_lo + tdel_stop(izb) = tstop(izb) - t(izb) + tdel_dy(izb) = 0.0 + tdel_dt9(izb) = 0.0 + If ( mask_init(izb) ) Then + intso(izb) = 0 + If ( tdelstart(izb) < 0.0 ) Then + changest = 0.5 + Else + changest = 0.01 + EndIf - tdel_stop = tstop(izb) - t(izb) - tdel_fine = 0.0 - tdel_t9 = 0.0 - If ( tdelstart(izb) < 0.0 ) Then - changest = 0.5 + ! For unevolved initial abundances, as often found in test problems, derivatives may + ! produce large jumps from zero abundance in the first timestep, though generally + ! inconsequential; an additional factor, changest, is used to increase accuracy. + changey = changest*min(0.1,changemx) + changet9 = changest*min(0.01,changemxt) + tdel_init = min(changest*abs(tdelstart(izb)),tdel_stop(izb)) + If ( nh(izb) > 2 ) Then + tdel_next(izb) = 1.0e-4*tdel_stop(izb) + Else + tdel_next(izb) = tdel_stop(izb) + EndIf Else - changest = 0.01 + changey = changemx + changet9 = changemxt + tdel_init = 0.0 EndIf - tdel_init = min(changest*abs(tdelstart(izb)),tdel_stop) - - ! If this is not the initial timestep, calculate timestep from changes in last timestep. - If ( tdel_old(izb) > 0.0 ) Then - Where ( y(:,izb) > yacc ) - ydotoy = abs((y(:,izb)-yo(:,izb))/y(:,izb)) - ElseWhere ( y(:,izb) > ymin ) - ydotoy = abs((y(:,izb)-yo(:,izb))/yacc) - ElseWhere - ydotoy = 0.0 - EndWhere - ints(izb) = maxloc(ydotoy,dim=1) - If ( abs(ydotoy(ints(izb))) > 0.0 ) Then - tdel_deriv = tdel_old(izb) * changemx/ydotoy(ints(izb)) - Else - tdel_deriv = tdel_next(izb) - EndIf - If ( iheat > 0 ) Then - t9dotot9 = abs((t9(izb)-t9o(izb))/t9(izb)) - If ( t9dotot9 > 0.0 ) Then - tdel_t9 = tdel_old(izb) * changemxt/t9dotot9 + ! Estimate timestep from relevant timescales + If ( tdel_old(izb) >= 0.0 ) Then + ! Calculate timescales for abundance changes, Y/(dY/dt) + Do k = 1, ny + If ( y(k,izb) > ymin ) Then + yfloor = max(y(k,izb),yacc) + If ( mask_init(izb) ) Then + ! Calculate timescale directly from derivatives. + rtau_y(k) = abs(ydot(k,izb)/yfloor) + Else + ! Calculate timescale from changes in last timestep. + rtau_y(k) = abs((y(k,izb)-yo(k,izb))/yfloor) / tdel_old(izb) + EndIf Else - tdel_t9 = tdel_next(izb) + rtau_y(k) = 0.0 EndIf - Else - tdel_t9 = tdel_next(izb) - EndIf - tdel(izb) = min( tdel_deriv, tdel_stop, tdel_next(izb), tdel_t9 ) - - ! If this is an initial timestep, yo does not exist, so calculate timestep directly from derivatives. - ElseIf ( abs(tdel_old(izb)) < tiny(0.0) ) Then - If ( nh(izb) > 2 ) tdel_stop = 1.0e-4*tdel_stop - intso(izb) = 0 - - Where ( y(:,izb) > yacc ) - ydotoy = abs(ydot(:,izb)/y(:,izb)) - ElseWhere ( y(:,izb) > ymin ) - ydotoy = abs(ydot(:,izb)/yacc) - ElseWhere - ydotoy = 0.0 - EndWhere + EndDo ! If derivatives are non-zero, use Y/(dY/dt). - ints(izb) = maxloc(ydotoy,dim=1) - If ( abs(ydotoy(ints(izb))) > 0.0 ) Then - tdel_deriv = changemx/ydotoy(ints(izb)) - - ! For unevolved initial abundances, as often found in test problems, - ! tdel_deriv may produce large jumps from zero abundance in the first - ! timestep. While generally inconsequential, tdel_fine limits these - ! to the accuracy abundance limit. - tdel_fine = changest*min(0.1d0,changemx)/ydotoy(ints(izb)) - - ! If derivatives are zero, take a small step. + ints(izb) = maxloc(rtau_y,dim=1) + If ( abs(rtau_y(ints(izb))) > 0.0 ) Then + tdel_dy(izb) = changey/rtau_y(ints(izb)) Else - tdel_deriv = tdel_stop - tdel_fine = tdel_stop + tdel_dy(izb) = tdel_next(izb) EndIf + + ! Calculate timescale for temperature changes, T/(dT/dt) If ( iheat > 0 ) Then - tdel_t9 = changest*min(0.01d0,changemxt) * abs(t9(izb)/t9dot(izb)) + If ( mask_init(izb) ) Then + rtau_t9 = abs(t9dot(izb)/t9(izb)) + Else + rtau_t9 = abs((t9(izb)-t9o(izb))/t9(izb)) / tdel_old(izb) + EndIf + Else + rtau_t9 = 0.0 + EndIf + If ( rtau_t9 > 0.0 ) Then + tdel_dt9(izb) = changet9/rtau_t9 Else - tdel_t9 = tdel_stop + tdel_dt9(izb) = tdel_next(izb) EndIf - tdel(izb) = min( tdel_stop, tdel_deriv, tdel_fine, tdel_t9 ) + + ! Estimate timestep + tdel(izb) = min( tdel_stop(izb), tdel_next(izb), tdel_dy(izb), tdel_dt9(izb) ) ! Use the user-defined initial timestep if it is larger than the estimated timestep tdel(izb) = max( tdel_init, tdel(izb) ) - ! Keep timestep constant Else + ! Keep timestep constant tdel(izb) = -tdel_old(izb) EndIf - If ( idiag >= 2 ) Write(lun_diag,"(a4,2i5,7es23.8,i5)") & - & 'tdel',kstep,izone,tdel(izb),tdel_deriv,tdel_old(izb),tdel_stop,tdel_next(izb),tdel_fine,tdel_t9,ints(izb) - !If ( idiag >= 2 ) Write(lun_diag,"(a5,i4,2es12.4)") (nname(k),k,y(k,izb),ydotoy(k),k=1,ny) - - ! Retain the index of the species setting the timestep - If ( ints(izb) /= intso(izb) ) Then - If ( idiag >= 2 ) Write(lun_diag,"(a4,a5,3es23.15)") 'ITC ',nname(ints(izb)),y(ints(izb),izb),t(izb),tdel(izb) - intso(izb) = ints(izb) - EndIf + ! Update trial time + tt(izb) = min(t(izb) + tdel(izb), tstop(izb)) EndIf EndDo + If ( idiag >= 2 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + izone = izb + szbatch - zb_lo + Write(lun_diag,"(a4,2i5,6es12.4,i5)") & + & 'tdel',kstep,izone,tdel(izb),tdel_old(izb),tdel_stop(izb),tdel_next(izb), & + & tdel_dy(izb),tdel_dt9(izb),ints(izb) + !Write(lun_diag,"(a5,i4,2es12.4)") (nname(k),k,y(k,izb),rtau_y(k),k=1,ny) + + ! Retain the index of the species setting the timestep + If ( ints(izb) /= intso(izb) ) Then + Write(lun_diag,"(a4,a5,3es23.15)") 'ITC ',nname(ints(izb)),y(ints(izb),izb),t(izb),tdel(izb) + intso(izb) = ints(izb) + EndIf + EndIf + EndDo + EndIf + ! For varying temperature and density, capture thermodynamic features If ( iheat == 0 ) Then ! Make sure to not skip any features in the temperature or density profiles by checking ! for profile monotonicity between t and t+del Do izb = zb_lo, zb_hi - mask_profile(izb) = ( mask(izb) .and. nh(izb) > 1 ) - If ( mask_profile(izb) ) Then - tt(izb) = min(t(izb) + tdel(izb), tstop(izb)) - EndIf - EndDo - Call t9rhofind(kstep,tt(zb_lo:zb_hi),ntt(zb_lo:zb_hi),t9t(zb_lo:zb_hi),rhot(zb_lo:zb_hi), & - & mask_in = mask_profile) - Do izb = zb_lo, zb_hi - If ( mask_profile(izb) .and. ntt(izb)-1 > nt(izb) ) Then - Do j = nt(izb), ntt(izb)-1 - If ( t9h(j,izb) > t9(izb) .and. t9h(j,izb) > t9t(izb) ) Then - tdel(izb) = th(j,izb) - t(izb) - Exit - ElseIf ( t9h(j,izb) < t9(izb) .and. t9h(j,izb) < t9t(izb) ) Then - tdel(izb) = th(j,izb) - t(izb) - Exit - ElseIf ( rhoh(j,izb) > rho(izb) .and. rhoh(j,izb) > rhot(izb) ) Then - tdel(izb) = th(j,izb) - t(izb) - Exit - ElseIf ( rhoh(j,izb) < rho(izb) .and. rhoh(j,izb) < rhot(izb) ) Then - tdel(izb) = th(j,izb) - t(izb) + If ( mask(izb) .and. nh(izb) > 1 ) Then + Call t9rhofind(kstep,tt(izb),ntt(izb),t9t(izb),rhot(izb), & + & nh(izb),th(:,izb),t9h(:,izb),rhoh(:,izb)) + If ( ntt(izb)-1 > nt(izb) ) Then + Do j = nt(izb), ntt(izb)-1 + If ( t9h(j,izb) > t9(izb) .and. t9h(j,izb) > t9t(izb) ) Then + tdel(izb) = th(j,izb) - t(izb) + Exit + ElseIf ( t9h(j,izb) < t9(izb) .and. t9h(j,izb) < t9t(izb) ) Then + tdel(izb) = th(j,izb) - t(izb) + Exit + ElseIf ( rhoh(j,izb) > rho(izb) .and. rhoh(j,izb) > rhot(izb) ) Then + tdel(izb) = th(j,izb) - t(izb) + Exit + ElseIf ( rhoh(j,izb) < rho(izb) .and. rhoh(j,izb) < rhot(izb) ) Then + tdel(izb) = th(j,izb) - t(izb) + Exit + EndIf + EndDo + tt(izb) = min(t(izb) + tdel(izb), tstop(izb)) + EndIf + + ! Limit timestep if fractional density change is larger than changeth (10% by default) + ! or fraction temperature change is larger than 0.1*changeth (1% by default) + dtherm(izb) = 0.0 + Do i = 1, 10 + Call t9rhofind(kstep,tt(izb),ntt(izb),t9t(izb),rhot(izb), & + & nh(izb),th(:,izb),t9h(:,izb),rhoh(:,izb)) + If ( t9(izb) > 0.0 ) Then + dtherm(izb) = 10.0*abs(t9t(izb)-t9(izb))/t9(izb) + abs(rhot(izb)-rho(izb))/rho(izb) + EndIf + If ( dtherm(izb) < changeth ) Then Exit + Else + tdel(izb) = 0.5*tdel(izb) + tt(izb) = min(t(izb) + tdel(izb), tstop(izb)) EndIf EndDo EndIf EndDo - - ! Limit timestep if fractional density change is larger than changeth (10% by default) - ! or fraction temperature change is larger than 0.1*changeth (1% by default) - dtherm = 0.0 - Do i = 1, 10 - Do izb = zb_lo, zb_hi - If ( mask_profile(izb) ) Then - tt(izb) = min(t(izb) + tdel(izb), tstop(izb)) - EndIf - EndDo - Call t9rhofind(kstep,tt(zb_lo:zb_hi),ntt(zb_lo:zb_hi),t9t(zb_lo:zb_hi),rhot(zb_lo:zb_hi), & - & mask_in = mask_profile) - Do izb = zb_lo, zb_hi - If ( mask_profile(izb) .and. t9(izb) > 0.0 ) Then - dtherm(izb) = 10.0*abs(t9t(izb)-t9(izb))/t9(izb) + abs(rhot(izb)-rho(izb))/rho(izb) - EndIf - EndDo - If ( all( dtherm < changeth ) ) Exit - Do izb = zb_lo, zb_hi - If ( dtherm(izb) >= changeth ) Then - tdel(izb) = 0.5*tdel(izb) - EndIf - EndDo - EndDo If ( idiag >= 2 ) Then Do izb = zb_lo, zb_hi - If ( mask_profile(izb) ) Then + If ( mask(izb) .and. nh(izb) > 1 ) Then izone = izb + szbatch - zb_lo - If ( dtherm(izb) >= changeth ) Write(lun_diag,"(a,i2,a,i5,3es23.15)") & - & 'Error in Thermo variations after ',i,' reductions',izone,tdel(izb),t9t(izb),rhot(izb) + If ( dtherm(izb) >= changeth ) Then + Write(lun_diag,"(a,i2,a,i5,3es24.16)") & + & 'Error in Thermo variations after ',i,' reductions',izone,tdel(izb),t9t(izb),rhot(izb) + Call xnet_terminate('Error in Thermo variations') + EndIf Write(lun_diag,"(a5,2i5,2es23.15)") 'T9del',kstep,izone,tdel(izb),dtherm(izb) EndIf EndDo EndIf EndIf - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - tt(izb) = min(t(izb) + tdel(izb), tstop(izb)) - EndIf - EndDo - Return End Subroutine timestep @@ -239,12 +239,12 @@ Subroutine update_iweak(t9,mask_in) !----------------------------------------------------------------------------------------------- ! This routine updates the value of the iweak flag to control the treatment of strong reactions. !----------------------------------------------------------------------------------------------- - Use xnet_controls, Only: iweak0, iweak, lun_stdout, t9min, nzevolve, zb_lo, zb_hi, lzactive + Use xnet_controls, Only: iweak0, iweak, lun_stdout, t9min, zb_lo, zb_hi, lzactive Use xnet_types, Only: dp Implicit None ! Input variables - Real(dp), Intent(in) :: t9(nzevolve) + Real(dp), Intent(in) :: t9(zb_lo:zb_hi) ! Optional variables Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) @@ -305,12 +305,9 @@ Subroutine update_eos(mask_in) start_timer = xnet_wtime() timer_eos = timer_eos - start_timer - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - call eos_interface(t9t(izb),rhot(izb),yt(:,izb),yet(izb),cv(izb),etae(izb),detaedt9(izb), & - & xext(izb),aext(izb),zext(izb)) - EndIf - EndDo + Call eos_interface(t9t(zb_lo:zb_hi),rhot(zb_lo:zb_hi),yt(:,zb_lo:zb_hi), & + & yet(zb_lo:zb_hi),cv(zb_lo:zb_hi),etae(zb_lo:zb_hi),detaedt9(zb_lo:zb_hi), & + & xext(zb_lo:zb_hi),aext(zb_lo:zb_hi),zext(zb_lo:zb_hi),mask_in = mask) stop_timer = xnet_wtime() timer_eos = timer_eos + stop_timer @@ -329,7 +326,8 @@ Subroutine yderiv(mask_in) & n22, n31, n32, n33, n41, n42, n43, n44, csect1, csect2, csect3, csect4, n1i, n2i, n3i, n4i Use xnet_abundances, Only: yt, ydot Use xnet_conditions, Only: cv, t9t, t9dot - Use xnet_controls, Only: idiag, iheat, ktot, lun_diag, nzbatchmx, szbatch, zb_lo, zb_hi, lzactive + Use xnet_controls, Only: idiag, iheat, ktot, lun_diag, nzbatchmx, szbatch, zb_lo, zb_hi, & + & lzactive Use xnet_timers, Only: xnet_wtime, start_timer, stop_timer, timer_deriv Use xnet_types, Only: dp Implicit None @@ -339,7 +337,7 @@ Subroutine yderiv(mask_in) ! Local variables Integer :: i, j, i0, i1, la1, le1, la2, le2, la3, le3, la4, le4, izb, izone - Real(dp) :: s1, s2, s3, s4, s11, s22, s33, s44 + Real(dp) :: s1, s2, s3, s4, s11, s22, s33, s44, sdot Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -352,16 +350,11 @@ Subroutine yderiv(mask_in) start_timer = xnet_wtime() timer_deriv = timer_deriv - start_timer + ! From the cross sections and the counting array, calculate the reaction rates + ! Calculate Ydot and T9dot for each nucleus, summing over the reactions which affect it. Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - ! From the cross sections and the counting array, calculate the reaction rates - b1(:,izb) = a1*csect1(mu1,izb) - b2(:,izb) = a2*csect2(mu2,izb) - b3(:,izb) = a3*csect3(mu3,izb) - b4(:,izb) = a4*csect4(mu4,izb) - - ! Calculate Ydot for each nucleus, summing over the reactions which affect it. Do i0 = 1, ny la1 = la(1,i0) la2 = la(2,i0) @@ -375,6 +368,7 @@ Subroutine yderiv(mask_in) ! Sum over the reactions with 1 reactant s1 = 0.0 Do i1 = la1, le1 + b1(i1,izb) = a1(i1)*csect1(mu1(i1),izb) s11 = b1(i1,izb)*yt(n11(i1),izb) s1 = s1 + s11 EndDo @@ -382,6 +376,7 @@ Subroutine yderiv(mask_in) ! Sum over the reactions with 2 reactants s2 = 0.0 Do i1 = la2, le2 + b2(i1,izb) = a2(i1)*csect2(mu2(i1),izb) s22 = b2(i1,izb)*yt(n21(i1),izb)*yt(n22(i1),izb) s2 = s2 + s22 EndDo @@ -389,6 +384,7 @@ Subroutine yderiv(mask_in) ! Sum over the reactions with 3 reactants s3 = 0.0 Do i1 = la3, le3 + b3(i1,izb) = a3(i1)*csect3(mu3(i1),izb) s33 = b3(i1,izb)*yt(n31(i1),izb)*yt(n32(i1),izb)*yt(n33(i1),izb) s3 = s3 + s33 EndDo @@ -396,6 +392,7 @@ Subroutine yderiv(mask_in) ! Sum over the reactions with 4 reactants s4 = 0.0 Do i1 = la4, le4 + b4(i1,izb) = a4(i1)*csect4(mu4(i1),izb) s44 = b4(i1,izb)*yt(n41(i1),izb)*yt(n42(i1),izb)*yt(n43(i1),izb)*yt(n44(i1),izb) s4 = s4 + s44 EndDo @@ -403,23 +400,27 @@ Subroutine yderiv(mask_in) ! Sum the 4 components of Ydot ydot(i0,izb) = s1 + s2 + s3 + s4 EndDo + EndIf + EndDo - If ( iheat > 0 ) Then - ! Surprisingly, this seems to perform better than the DGEMV below - s1 = 0.0 + If ( iheat > 0 ) Then + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + sdot = 0.0 Do i0 = 1, ny - s1 = s1 + mex(i0)*ydot(i0,izb) + sdot = sdot - mex(i0)*ydot(i0,izb) / cv(izb) EndDo - t9dot(izb) = -s1 / cv(izb) + t9dot(izb) = sdot EndIf + EndDo + EndIf + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + ktot(4,izb) = ktot(4,izb) + 1 EndIf EndDo - !If ( iheat > 0 ) Then - ! Call dgemv('T',ny,nzbatchmx,1.0,ydot(1,zb_lo),ny,mex,1,0.0,t9dot(zb_lo),1) - ! Where ( mask ) - ! t9dot = -t9dot / cv - ! EndWhere - !EndIf ! Separate loop for diagnostics so compiler can properly vectorize If ( idiag >= 5 ) Then @@ -482,12 +483,6 @@ Subroutine yderiv(mask_in) EndDo EndIf - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - ktot(4,izb) = ktot(4,izb) + 1 - EndIf - EndDo - stop_timer = xnet_wtime() timer_deriv = timer_deriv + stop_timer @@ -510,18 +505,23 @@ Subroutine cross_sect(mask_in) Use xnet_screening, Only: h1, h2, h3, h4, dh1dt9, dh2dt9, dh3dt9, dh4dt9, screening Use xnet_timers, Only: xnet_wtime, start_timer, stop_timer, timer_csect Use xnet_types, Only: dp + Use xnet_util, Only: safe_exp Implicit None ! Optional variables Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables + Logical, Parameter :: use_blas = .true. Integer, Parameter :: dgemm_nzbatch = 200 ! Min number of zones to use dgemm, otherwise use dgemv Real(dp) :: t09(7,zb_lo:zb_hi), dt09(7,zb_lo:zb_hi) Real(dp) :: ene(zb_lo:zb_hi), ytot, abar, zbar, z2bar, zibar Real(dp) :: ascrn, rhot2, rhot3 Real(dp) :: rpf1, rpf2, rpf3, rpf4 Real(dp) :: dlnrpf1dt9, dlnrpf2dt9, dlnrpf3dt9, dlnrpf4dt9 + Real(dp) :: p1, p2, p3, p4, s1, s2, s3, s4 + Real(dp) :: lambda1, lambda2, lambda3, lambda4 + Real(dp) :: dlam1dt9, dlam2dt9, dlam3dt9, dlam4dt9 Integer :: j, k, izb, izone, nzmask Integer :: nr1, nr2, nr3, nr4 Logical, Pointer :: mask(:) @@ -532,6 +532,7 @@ Subroutine cross_sect(mask_in) mask(zb_lo:) => lzactive(zb_lo:zb_hi) EndIf If ( .not. any(mask) ) Return + nzmask = count(mask) nr1 = nreac(1) nr2 = nreac(2) @@ -541,9 +542,6 @@ Subroutine cross_sect(mask_in) ! Update thermodynamic state Call update_eos(mask_in = mask) - ! Check for any changes to iweak - Call update_iweak(t9t,mask_in = mask) - ! Calculate the screening terms If ( iscrn >= 1 ) Then ascrn = 1.0 @@ -555,9 +553,13 @@ Subroutine cross_sect(mask_in) start_timer = xnet_wtime() timer_csect = timer_csect - start_timer + ! Check for any changes to iweak + Call update_iweak(t9t(zb_lo:zb_hi),mask_in = mask) + + ! Calculate partition functions for each nucleus at t9t + Call partf(t9t(zb_lo:zb_hi),mask_in = mask) + ! Calculate necessary thermodynamic moments - ene = 0.0 - t09 = 0.0 Do izb = zb_lo, zb_hi If ( mask(izb) ) Then ene(izb) = yet(izb)*rhot(izb) @@ -569,275 +571,251 @@ Subroutine cross_sect(mask_in) t09(5,izb) = t9t(izb) t09(6,izb) = t9t(izb)**(+5.0/3.0) t09(7,izb) = log(t9t(izb)) + + ! Calculate reaction rate coefficient derivatives + dt09(1,izb) = 0.0 + dt09(2,izb) = -t9t(izb)**(-2) + dt09(3,izb) = -t9t(izb)**(-4.0/3.0) / 3.0 + dt09(4,izb) = +t9t(izb)**(-2.0/3.0) / 3.0 + dt09(5,izb) = 1.0 + dt09(6,izb) = +t9t(izb)**(+2.0/3.0) * 5.0/3.0 + dt09(7,izb) = 1.0 / t9t(izb) EndIf EndDo - ! Calculate the REACLIB exponent polynomials, adding screening terms - nzmask = count(mask) - If ( nzmask >= dgemm_nzbatch ) Then - If ( nr1 > 0 ) Call dgemm('T','N',nr1,nzbatchmx,7,1.0,rc1,7,t09,7,ascrn,h1(:,zb_lo:zb_hi),nr1) - If ( nr2 > 0 ) Call dgemm('T','N',nr2,nzbatchmx,7,1.0,rc2,7,t09,7,ascrn,h2(:,zb_lo:zb_hi),nr2) - If ( nr3 > 0 ) Call dgemm('T','N',nr3,nzbatchmx,7,1.0,rc3,7,t09,7,ascrn,h3(:,zb_lo:zb_hi),nr3) - If ( nr4 > 0 ) Call dgemm('T','N',nr4,nzbatchmx,7,1.0,rc4,7,t09,7,ascrn,h4(:,zb_lo:zb_hi),nr4) + ! If there are any FFN reactions, calculate their rates + If ( nffn > 0 ) Call ffn_rate(nffn,t9t(zb_lo:zb_hi),ene, & + & rffn(:,zb_lo:zb_hi),dlnrffndt9(:,zb_lo:zb_hi),mask_in = mask) + + ! If there are any neutrino-nucleus reactions, calculate their rates + If ( nnnu > 0 ) Call nnu_rate(nnnu,tt(zb_lo:zb_hi),rnnu(:,:,zb_lo:zb_hi),mask_in = mask) + + ! Calculate the REACLIB exponent polynomials and derivatives, adding screening terms + If ( use_blas ) Then + If ( nzmask >= dgemm_nzbatch ) Then + If ( nr1 > 0 ) Call dgemm('T','N',nr1,nzbatchmx,7,1.0,rc1,7,t09,7,ascrn,h1(:,zb_lo:zb_hi),nr1) + If ( nr2 > 0 ) Call dgemm('T','N',nr2,nzbatchmx,7,1.0,rc2,7,t09,7,ascrn,h2(:,zb_lo:zb_hi),nr2) + If ( nr3 > 0 ) Call dgemm('T','N',nr3,nzbatchmx,7,1.0,rc3,7,t09,7,ascrn,h3(:,zb_lo:zb_hi),nr3) + If ( nr4 > 0 ) Call dgemm('T','N',nr4,nzbatchmx,7,1.0,rc4,7,t09,7,ascrn,h4(:,zb_lo:zb_hi),nr4) + + If ( iheat > 0 ) Then + If ( nr1 > 0 ) Call dgemm('T','N',nr1,nzbatchmx,7,1.0,rc1,7,dt09,7,ascrn,dh1dt9(:,zb_lo:zb_hi),nr1) + If ( nr2 > 0 ) Call dgemm('T','N',nr2,nzbatchmx,7,1.0,rc2,7,dt09,7,ascrn,dh2dt9(:,zb_lo:zb_hi),nr2) + If ( nr3 > 0 ) Call dgemm('T','N',nr3,nzbatchmx,7,1.0,rc3,7,dt09,7,ascrn,dh3dt9(:,zb_lo:zb_hi),nr3) + If ( nr4 > 0 ) Call dgemm('T','N',nr4,nzbatchmx,7,1.0,rc4,7,dt09,7,ascrn,dh4dt9(:,zb_lo:zb_hi),nr4) + EndIf + Else + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + If ( nr1 > 0 ) Call dgemv('T',7,nr1,1.0,rc1,7,t09(:,izb),1,ascrn,h1(:,izb),1) + If ( nr2 > 0 ) Call dgemv('T',7,nr2,1.0,rc2,7,t09(:,izb),1,ascrn,h2(:,izb),1) + If ( nr3 > 0 ) Call dgemv('T',7,nr3,1.0,rc3,7,t09(:,izb),1,ascrn,h3(:,izb),1) + If ( nr4 > 0 ) Call dgemv('T',7,nr4,1.0,rc4,7,t09(:,izb),1,ascrn,h4(:,izb),1) + If ( iheat > 0 ) Then + If ( nr1 > 0 ) Call dgemv('T',7,nr1,1.0,rc1,7,dt09(:,izb),1,ascrn,dh1dt9(:,izb),1) + If ( nr2 > 0 ) Call dgemv('T',7,nr2,1.0,rc2,7,dt09(:,izb),1,ascrn,dh2dt9(:,izb),1) + If ( nr3 > 0 ) Call dgemv('T',7,nr3,1.0,rc3,7,dt09(:,izb),1,ascrn,dh3dt9(:,izb),1) + If ( nr4 > 0 ) Call dgemv('T',7,nr4,1.0,rc4,7,dt09(:,izb),1,ascrn,dh4dt9(:,izb),1) + EndIf + EndIf + EndDo + EndIf Else Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - If ( nr1 > 0 ) Call dgemv('T',7,nr1,1.0,rc1,7,t09(:,izb),1,ascrn,h1(:,izb),1) - If ( nr2 > 0 ) Call dgemv('T',7,nr2,1.0,rc2,7,t09(:,izb),1,ascrn,h2(:,izb),1) - If ( nr3 > 0 ) Call dgemv('T',7,nr3,1.0,rc3,7,t09(:,izb),1,ascrn,h3(:,izb),1) - If ( nr4 > 0 ) Call dgemv('T',7,nr4,1.0,rc4,7,t09(:,izb),1,ascrn,h4(:,izb),1) + Do k = 1, nr1 + lambda1 = 0.0 + dlam1dt9 = 0.0 + Do j = 1, 7 + lambda1 = lambda1 + rc1(j,k)*t09(j,izb) + dlam1dt9 = dlam1dt9 + rc1(j,k)*dt09(j,izb) + EndDo + h1(k,izb) = ascrn*h1(k,izb) + lambda1 + dh1dt9(k,izb) = ascrn*dh1dt9(k,izb) + dlam1dt9 + EndDo + Do k = 1, nr2 + lambda2 = 0.0 + dlam2dt9 = 0.0 + Do j = 1, 7 + lambda2 = lambda2 + rc2(j,k)*t09(j,izb) + dlam2dt9 = dlam2dt9 + rc2(j,k)*dt09(j,izb) + EndDo + h2(k,izb) = ascrn*h2(k,izb) + lambda2 + dh2dt9(k,izb) = ascrn*dh2dt9(k,izb) + dlam2dt9 + EndDo + Do k = 1, nr3 + lambda3 = 0.0 + dlam3dt9 = 0.0 + Do j = 1, 7 + lambda3 = lambda3 + rc3(j,k)*t09(j,izb) + dlam3dt9 = dlam3dt9 + rc3(j,k)*dt09(j,izb) + EndDo + h3(k,izb) = ascrn*h3(k,izb) + lambda3 + dh3dt9(k,izb) = ascrn*dh3dt9(k,izb) + dlam3dt9 + EndDo + Do k = 1, nr4 + lambda4 = 0.0 + dlam4dt9 = 0.0 + Do j = 1, 7 + lambda4 = lambda4 + rc4(j,k)*t09(j,izb) + dlam4dt9 = dlam4dt9 + rc4(j,k)*dt09(j,izb) + EndDo + h4(k,izb) = ascrn*h4(k,izb) + lambda4 + dh4dt9(k,izb) = ascrn*dh4dt9(k,izb) + dlam4dt9 + EndDo EndIf EndDo EndIf - ! Calculate partition functions for each nucleus at t9t - Call partf(t9t,mask_in = mask) - - ! If there are any FFN reactions, calculate their rates - If ( nffn > 0 ) Call ffn_rate(nffn,t9t,ene,rffn,dlnrffndt9,mask_in = mask) - - ! If there are any neutrino-nucleus reactions, calculate their rates - If ( nnnu > 0 ) Call nnu_rate(nnnu,tt,rnnu,mask_in = mask) - ! Calculate cross sections Do izb = zb_lo, zb_hi If ( mask(izb) ) Then ! Calculate the csect for reactions with 1 reactant Do k = 1, nr1 - If ( irev1(k) == 1 ) Then - rpf1 = ( gg(n1i(2,k),izb) * gg(n1i(3,k),izb) * gg(n1i(4,k),izb) * gg(n1i(5,k),izb) ) & - & / ( gg(n1i(1,k),izb) ) + If ( (iweak(izb) < 0 .and. iwk1(k) == 0) .or. & + & (iweak(izb) == 0 .and. iwk1(k) /= 0) ) Then + csect1(k,izb) = 0.0 + dcsect1dt9(k,izb) = 0.0 Else rpf1 = 1.0 - EndIf - If ( iweak(izb) > 0 ) Then - If ( iwk1(k) == 0 .or. iwk1(k) == 4 ) Then - csect1(k,izb) = rpf1 * exp(h1(k,izb)) - ElseIf ( iwk1(k) == 1 ) Then - csect1(k,izb) = rpf1 * exp(h1(k,izb)) * ene(izb) - ElseIf ( iwk1(k) == 2 .or. iwk1(k) == 3 ) Then ! FFN reaction - csect1(k,izb) = rffn(iffn(k),izb) - ElseIf ( iwk1(k) == 7 ) Then ! Electron neutrino capture - csect1(k,izb) = rpf1 * rnnu(innu(k),1,izb) - ElseIf ( iwk1(k) == 8 ) Then ! Electron anti-neutrino capture - csect1(k,izb) = rpf1 * rnnu(innu(k),2,izb) + dlnrpf1dt9 = 0.0 + If ( irev1(k) == 1 ) Then + p1 = gg(n1i(1,k),izb) + s1 = dlngdt9(n1i(1,k),izb) + Do j = 2, 5 + rpf1 = rpf1 * gg(n1i(j,k),izb) + dlnrpf1dt9 = dlnrpf1dt9 + dlngdt9(n1i(j,k),izb) + EndDo + rpf1 = rpf1 / p1 + dlnrpf1dt9 = dlnrpf1dt9 - s1 EndIf - ElseIf ( iweak(izb) < 0 ) Then - If ( iwk1(k) == 0 ) Then - csect1(k,izb) = 0.0 - ElseIf ( iwk1(k) == 1 ) Then - csect1(k,izb) = rpf1 * exp(h1(k,izb)) * ene(izb) - ElseIf ( iwk1(k) == 4 ) Then - csect1(k,izb) = rpf1 * exp(h1(k,izb)) + If ( iwk1(k) == 1 ) Then + csect1(k,izb) = ene(izb) * rpf1 * safe_exp(h1(k,izb)) ElseIf ( iwk1(k) == 2 .or. iwk1(k) == 3 ) Then ! FFN reaction csect1(k,izb) = rffn(iffn(k),izb) ElseIf ( iwk1(k) == 7 ) Then ! Electron neutrino capture csect1(k,izb) = rpf1 * rnnu(innu(k),1,izb) ElseIf ( iwk1(k) == 8 ) Then ! Electron anti-neutrino capture csect1(k,izb) = rpf1 * rnnu(innu(k),2,izb) - EndIf - Else - If ( iwk1(k) == 0 ) Then - csect1(k,izb) = rpf1 * exp(h1(k,izb)) Else - csect1(k,izb) = 0.0 + csect1(k,izb) = rpf1 * safe_exp(h1(k,izb)) + EndIf + If ( iheat > 0 ) Then + If ( iwk1(k) == 2 .or. iwk1(k) == 3 ) Then ! FFN reaction + dcsect1dt9(k,izb) = csect1(k,izb)*dlnrffndt9(iffn(k),izb) + ElseIf ( iwk1(k) == 7 .or. iwk1(k) == 8 ) Then ! NNU reaction + dcsect1dt9(k,izb) = csect1(k,izb)*dlnrpf1dt9 + Else + dcsect1dt9(k,izb) = csect1(k,izb)*(dh1dt9(k,izb)+dlnrpf1dt9) + EndIf EndIf EndIf EndDo ! Calculate the csect for reactions with 2 reactants Do k = 1, nr2 - If ( irev2(k) == 1 ) Then - rpf2 = ( gg(n2i(3,k),izb) * gg(n2i(4,k),izb) * gg(n2i(5,k),izb) * gg(n2i(6,k),izb) ) & - & / ( gg(n2i(1,k),izb) * gg(n2i(2,k),izb) ) + If ( (iweak(izb) < 0 .and. iwk2(k) == 0) .or. & + & (iweak(izb) == 0 .and. iwk2(k) /= 0) ) Then + csect2(k,izb) = 0.0 + dcsect2dt9(k,izb) = 0.0 Else rpf2 = 1.0 - EndIf - If ( iweak(izb) > 0 ) Then - If ( iwk2(k) == 1 ) Then - csect2(k,izb) = rhot(izb) * rpf2 * exp(h2(k,izb)) * ene(izb) - Else - csect2(k,izb) = rhot(izb) * rpf2 * exp(h2(k,izb)) - EndIf - ElseIf ( iweak(izb) < 0 ) Then - If ( iwk2(k) == 0 ) Then - csect2(k,izb) = 0.0 - ElseIf ( iwk2(k) == 1 ) Then - csect2(k,izb) = rhot(izb) * rpf2 * exp(h2(k,izb)) * ene(izb) - Else - csect2(k,izb) = rhot(izb) * rpf2 * exp(h2(k,izb)) + dlnrpf2dt9 = 0.0 + If ( irev2(k) == 1 ) Then + p2 = 1.0 + s2 = 0.0 + Do j = 1, 2 + p2 = p2 * gg(n2i(j,k),izb) + s2 = s2 + dlngdt9(n2i(j,k),izb) + EndDo + Do j = 3, 6 + rpf2 = rpf2 * gg(n2i(j,k),izb) + dlnrpf2dt9 = dlnrpf2dt9 + dlngdt9(n2i(j,k),izb) + EndDo + rpf2 = rpf2 / p2 + dlnrpf2dt9 = dlnrpf2dt9 - s2 EndIf - Else - If ( iwk2(k) == 0 ) Then - csect2(k,izb) = rhot(izb) * rpf2 * exp(h2(k,izb)) + If ( iwk2(k) == 1 ) Then + csect2(k,izb) = rhot(izb) * ene(izb) * rpf2 * safe_exp(h2(k,izb)) Else - csect2(k,izb) = 0.0 + csect2(k,izb) = rhot(izb) * rpf2 * safe_exp(h2(k,izb)) EndIf + If ( iheat > 0 ) dcsect2dt9(k,izb) = csect2(k,izb)*(dh2dt9(k,izb)+dlnrpf2dt9) EndIf EndDo ! Calculate the csect for reactions with 3 reactants Do k = 1, nr3 - If ( irev3(k) == 1 ) Then - rpf3 = ( gg(n3i(4,k),izb) * gg(n3i(5,k),izb) * gg(n3i(6,k),izb) ) & - & / ( gg(n3i(1,k),izb) * gg(n3i(2,k),izb) * gg(n3i(3,k),izb) ) + If ( (iweak(izb) < 0 .and. iwk3(k) == 0) .or. & + & (iweak(izb) == 0 .and. iwk3(k) /= 0) ) Then + csect3(k,izb) = 0.0 + dcsect3dt9(k,izb) = 0.0 Else rpf3 = 1.0 - EndIf - rhot2 = rhot(izb)**2 - If ( iweak(izb) > 0 ) Then - If ( iwk3(k) == 1 ) Then - csect3(k,izb) = rhot2 * rpf3 * exp(h3(k,izb)) * ene(izb) - Else - csect3(k,izb) = rhot2 * rpf3 * exp(h3(k,izb)) - EndIf - ElseIf ( iweak(izb) < 0 ) Then - If ( iwk3(k) == 0 ) Then - csect3(k,izb) = 0.0 - ElseIf ( iwk3(k) == 1 ) Then - csect3(k,izb) = rhot2 * rpf3 * exp(h3(k,izb)) * ene(izb) - Else - csect3(k,izb) = rhot2 * rpf3 * exp(h3(k,izb)) + dlnrpf3dt9 = 0.0 + If ( irev3(k) == 1 ) Then + p3 = 1.0 + s3 = 0.0 + Do j = 1, 3 + p3 = p3 * gg(n3i(j,k),izb) + s3 = s3 + dlngdt9(n3i(j,k),izb) + EndDo + Do j = 4, 6 + rpf3 = rpf3 * gg(n3i(j,k),izb) + dlnrpf3dt9 = dlnrpf3dt9 + dlngdt9(n3i(j,k),izb) + EndDo + rpf3 = rpf3 / p3 + dlnrpf3dt9 = dlnrpf3dt9 - s3 EndIf - Else - If ( iwk3(k) == 0 ) Then - csect3(k,izb) = rhot2 * rpf3 * exp(h3(k,izb)) + rhot2 = rhot(izb)**2 + If ( iwk3(k) == 1 ) Then + csect3(k,izb) = rhot2 * ene(izb) * rpf3 * safe_exp(h3(k,izb)) Else - csect3(k,izb) = 0.0 + csect3(k,izb) = rhot2 * rpf3 * safe_exp(h3(k,izb)) EndIf + If ( csect3(k,izb) < 1.0e-20 ) csect3(k,izb) = 0.0 + If ( iheat > 0 ) dcsect3dt9(k,izb) = csect3(k,izb)*(dh3dt9(k,izb)+dlnrpf3dt9) EndIf - If ( csect3(k,izb) < 1.0e-20 ) csect3(k,izb) = 0.0 EndDo ! Calculate the csect for reactions with 4 reactants Do k = 1, nr4 - If ( irev4(k) == 1 ) Then - rpf4 = ( gg(n4i(5,k),izb) * gg(n4i(6,k),izb) ) & - & / ( gg(n4i(1,k),izb) * gg(n4i(2,k),izb) * gg(n4i(3,k),izb) * gg(n4i(4,k),izb ) ) + If ( (iweak(izb) < 0 .and. iwk4(k) == 0) .or. & + & (iweak(izb) == 0 .and. iwk4(k) /= 0) ) Then + csect4(k,izb) = 0.0 + dcsect4dt9(k,izb) = 0.0 Else rpf4 = 1.0 - EndIf - rhot3 = rhot(izb)**3 - If ( iweak(izb) > 0 ) Then - If ( iwk4(k) == 1 ) Then - csect4(k,izb) = rhot3 * rpf4 * exp(h4(k,izb)) * ene(izb) - Else - csect4(k,izb) = rhot3 * rpf4 * exp(h4(k,izb)) - EndIf - ElseIf ( iweak(izb) < 0 ) Then - If ( iwk4(k) == 0 ) Then - csect4(k,izb) = 0.0 - ElseIf ( iwk4(k) == 1 ) Then - csect4(k,izb) = rhot3 * rpf4 * exp(h4(k,izb)) * ene(izb) - Else - csect4(k,izb) = rhot3 * rpf4 * exp(h4(k,izb)) + dlnrpf4dt9 = 0.0 + If ( irev4(k) == 1 ) Then + p4 = 1.0 + s4 = 0.0 + Do j = 1, 4 + p4 = p4 * gg(n4i(j,k),izb) + s4 = s4 + dlngdt9(n4i(j,k),izb) + EndDo + Do j = 5, 6 + rpf4 = rpf4 * gg(n4i(j,k),izb) + dlnrpf4dt9 = dlnrpf4dt9 + dlngdt9(n4i(j,k),izb) + EndDo + rpf4 = rpf4 / p4 + dlnrpf4dt9 = dlnrpf4dt9 - s4 EndIf - Else - If ( iwk4(k) == 0 ) Then - csect4(k,izb) = rhot3 * rpf4 * exp(h4(k,izb)) + rhot3 = rhot(izb)**3 + If ( iwk4(k) == 1 ) Then + csect4(k,izb) = rhot3 * ene(izb) * rpf4 * safe_exp(h4(k,izb)) Else - csect4(k,izb) = 0.0 + csect4(k,izb) = rhot3 * rpf4 * safe_exp(h4(k,izb)) EndIf + If ( iheat > 0 ) dcsect4dt9(k,izb) = csect4(k,izb)*(dh4dt9(k,izb)+dlnrpf4dt9) EndIf EndDo - EndIf - EndDo - - If ( iheat > 0 ) Then - - ! Calculate reaction rate coefficient derivatives - dt09 = 0.0 - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - dt09(1,izb) = 0.0 - dt09(2,izb) = -t9t(izb)**(-2) - dt09(3,izb) = -t9t(izb)**(-4.0/3.0) / 3.0 - dt09(4,izb) = +t9t(izb)**(-2.0/3.0) / 3.0 - dt09(5,izb) = 1.0 - dt09(6,izb) = +t9t(izb)**(+2.0/3.0) * 5.0/3.0 - dt09(7,izb) = 1.0 / t9t(izb) - EndIf - EndDo - ! Calculate the derivatives of REACLIB exponent polynomials, adding screening terms - If ( nzmask >= dgemm_nzbatch ) Then - If ( nr1 > 0 ) Call dgemm('T','N',nr1,nzbatchmx,7,1.0,rc1,7,dt09,7,ascrn,dh1dt9(:,zb_lo:zb_hi),nr1) - If ( nr2 > 0 ) Call dgemm('T','N',nr2,nzbatchmx,7,1.0,rc2,7,dt09,7,ascrn,dh2dt9(:,zb_lo:zb_hi),nr2) - If ( nr3 > 0 ) Call dgemm('T','N',nr3,nzbatchmx,7,1.0,rc3,7,dt09,7,ascrn,dh3dt9(:,zb_lo:zb_hi),nr3) - If ( nr4 > 0 ) Call dgemm('T','N',nr4,nzbatchmx,7,1.0,rc4,7,dt09,7,ascrn,dh4dt9(:,zb_lo:zb_hi),nr4) - Else - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - If ( nr1 > 0 ) Call dgemv('T',7,nr1,1.0,rc1,7,dt09(:,izb),1,ascrn,dh1dt9(:,izb),1) - If ( nr2 > 0 ) Call dgemv('T',7,nr2,1.0,rc2,7,dt09(:,izb),1,ascrn,dh2dt9(:,izb),1) - If ( nr3 > 0 ) Call dgemv('T',7,nr3,1.0,rc3,7,dt09(:,izb),1,ascrn,dh3dt9(:,izb),1) - If ( nr4 > 0 ) Call dgemv('T',7,nr4,1.0,rc4,7,dt09(:,izb),1,ascrn,dh4dt9(:,izb),1) - EndIf - EndDo + ! Increment counter + ktot(5,izb) = ktot(5,izb) + 1 EndIf - - ! Calculate cross-section derivatives - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - - ! 1-reactant reactions - Do k = 1, nr1 - If ( irev1(k) == 1 ) Then - dlnrpf1dt9 = ( dlngdt9(n1i(2,k),izb) + dlngdt9(n1i(3,k),izb) & - & + dlngdt9(n1i(4,k),izb) + dlngdt9(n1i(5,k),izb) ) & - & - ( dlngdt9(n1i(1,k),izb) ) - Else - dlnrpf1dt9 = 0.0 - EndIf - If ( iwk1(k) == 2 .or. iwk1(k) == 3 ) Then ! FFN reaction - dcsect1dt9(k,izb) = rffn(iffn(k),izb)*dlnrffndt9(iffn(k),izb) - ElseIf ( iwk1(k) == 7 ) Then ! Electron neutrino capture - dcsect1dt9(k,izb) = rnnu(innu(k),1,izb)*dlnrpf1dt9 - ElseIf ( iwk1(k) == 8 ) Then ! Electron anti-neutrino capture - dcsect1dt9(k,izb) = rnnu(innu(k),2,izb)*dlnrpf1dt9 - Else - dcsect1dt9(k,izb) = csect1(k,izb)*(dh1dt9(k,izb)+dlnrpf1dt9) - EndIf - EndDo - - ! 2-reactant reactions - Do k = 1, nr2 - If ( irev2(k) == 1 ) Then - dlnrpf2dt9 = ( dlngdt9(n2i(3,k),izb) + dlngdt9(n2i(4,k),izb) & - + dlngdt9(n2i(5,k),izb) + dlngdt9(n2i(6,k),izb) ) & - & - ( dlngdt9(n2i(1,k),izb) + dlngdt9(n2i(2,k),izb) ) - Else - dlnrpf2dt9 = 0.0 - EndIf - dcsect2dt9(k,izb) = csect2(k,izb)*(dh2dt9(k,izb)+dlnrpf2dt9) - EndDo - - ! 3-reactant reactions - Do k = 1, nr3 - If ( irev3(k) == 1 ) Then - dlnrpf3dt9 = ( dlngdt9(n3i(4,k),izb) + dlngdt9(n3i(5,k),izb) + dlngdt9(n3i(6,k),izb) ) & - & - ( dlngdt9(n3i(1,k),izb) + dlngdt9(n3i(2,k),izb) + dlngdt9(n3i(3,k),izb) ) - Else - dlnrpf3dt9 = 0.0 - EndIf - dcsect3dt9(k,izb) = csect3(k,izb)*(dh3dt9(k,izb)+dlnrpf3dt9) - EndDo - - ! 4-reactant reactions - Do k = 1, nr4 - If ( irev4(k) == 1 ) Then - dlnrpf4dt9 = ( dlngdt9(n4i(5,k),izb) + dlngdt9(n4i(6,k),izb) ) & - & - ( dlngdt9(n4i(1,k),izb) + dlngdt9(n4i(2,k),izb) & - & + dlngdt9(n4i(3,k),izb) + dlngdt9(n4i(4,k),izb) ) - Else - dlnrpf4dt9 = 0.0 - EndIf - dcsect4dt9(k,izb) = csect4(k,izb)*(dh4dt9(k,izb)+dlnrpf4dt9) - EndDo - EndIf - EndDo - EndIf + EndDo If ( idiag >= 6 ) Then Do izb = zb_lo, zb_hi @@ -882,12 +860,6 @@ Subroutine cross_sect(mask_in) EndDo EndIf - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - ktot(5,izb) = ktot(5,izb) + 1 - EndIf - EndDo - stop_timer = xnet_wtime() timer_csect = timer_csect + stop_timer diff --git a/source/xnet_integrate_bdf.f90 b/source/xnet_integrate_bdf.f90 index bf6dd650..1481079e 100644 --- a/source/xnet_integrate_bdf.f90 +++ b/source/xnet_integrate_bdf.f90 @@ -142,7 +142,7 @@ Subroutine solve_bdf(kstep,ierr) ! On output, = BDF_STATUS_FAIL if zone fails to converge ! Local variables - Integer :: kts, izb, izone + Integer :: kts, k, izb, izone ! Initial setup and allocations If ( kstep == 1 ) Call bdf_reset @@ -150,8 +150,9 @@ Subroutine solve_bdf(kstep,ierr) start_timer = xnet_wtime() timer_tstep = timer_tstep - start_timer - ! Copy status flag to local copy Do izb = zb_lo, zb_hi + + ! Copy status flag to local copy ierr_ts(izb) = ierr(izb) ! If the zone has previously converged or failed, do not iterate @@ -160,13 +161,11 @@ Subroutine solve_bdf(kstep,ierr) Else ierr_nr(izb) = BDF_STATUS_FAIL EndIf - EndDo - - ! Make adjustments to account for any change in stepsize or order from previous step - If ( kstep /= 1 ) Call bdf_adjust(kstep) - Do izb = zb_lo, zb_hi - ydot0(:,izb) = ydot(:,izb) + ! Save initial derivatives so they can be reloaded later if necessary + Do k = 1, ny + ydot0(k,izb) = ydot(k,izb) + EndDo If ( iheat > 0 ) t9dot0(izb) = t9dot(izb) ! Reset maximum stepsize for next step @@ -177,6 +176,10 @@ Subroutine solve_bdf(kstep,ierr) nerr_ts(izb) = 0 ncf_ts(izb) = 0 EndDo + + ! Make adjustments to account for any change in stepsize or order from previous step + If ( kstep /= 1 ) Call bdf_adjust(kstep) + Do kts = 1, max_it_ts Do izb = zb_lo, zb_hi @@ -204,7 +207,12 @@ Subroutine solve_bdf(kstep,ierr) tdel_next(izb) = eta(izb) * tdel(izb) If ( ierr_nr(izb) == BDF_STATUS_SUCCESS ) Then - yt(:,izb) = z(1:ny,0,izb) + ierr_ts(izb) = BDF_STATUS_SUCCESS + Do k = 1, ny + yt(k,izb) = z(k,0,izb) + yo(k,izb) = y(k,izb) + y(k,izb) = yt(k,izb) + EndDo If ( iheat > 0 ) t9t(izb) = z(neq,0,izb) nto(izb) = nt(izb) nt(izb) = ntt(izb) @@ -216,35 +224,31 @@ Subroutine solve_bdf(kstep,ierr) rho(izb) = rhot(izb) yeo(izb) = ye(izb) ye(izb) = yet(izb) - yo(:,izb) = y(:,izb) - y(:,izb) = yt(:,izb) + ElseIf ( ierr_nr(izb) == BDF_STATUS_FAIL .or. ierr_ts(izb) == BDF_STATUS_FAIL ) Then + ierr_ts(izb) = BDF_STATUS_FAIL EndIf ! Record iteration counters kmon(1,izb) = nit_ts(izb) ktot(1,izb) = ktot(1,izb) + nit_ts(izb) kmon(2,izb) = nit_nrslv(izb) + + ! Copy local status flag to output + ierr(izb) = ierr_ts(izb) EndDo ! Log TS success/failure Do izb = zb_lo, zb_hi izone = izb + szbatch - zb_lo If ( ierr_nr(izb) == BDF_STATUS_SUCCESS ) Then - ierr_ts(izb) = BDF_STATUS_SUCCESS If ( idiag >= 2 ) Write(lun_diag,"(a,2i5,2i3)") & 'BDF TS Success',kstep,izone,nit_ts(izb),nit_nr(izb) ElseIf ( ierr_nr(izb) == BDF_STATUS_FAIL .or. ierr_ts(izb) == BDF_STATUS_FAIL ) Then - ierr_ts(izb) = BDF_STATUS_FAIL If ( idiag >= 0 ) Write(lun_diag,"(a,2i5,4es12.4,2i3)") & 'BDF TS Fail',kstep,izone,t(izb),tdel(izb),t9t(izb),rhot(izb),nit_nr(izb),nit_ts(izb) EndIf EndDo - ! Copy local status flag to output - Do izb = zb_lo, zb_hi - ierr(izb) = ierr_ts(izb) - EndDo - stop_timer = xnet_wtime() timer_tstep = timer_tstep + stop_timer @@ -382,7 +386,7 @@ Subroutine bdf_reset Implicit None ! Local variables - Integer :: izb + Integer :: i, j, k, izb Do izb = zb_lo, zb_hi bdf_active(izb) = .false. @@ -403,37 +407,60 @@ Subroutine bdf_reset q_age(izb) = 0 nscon(izb) = 0 - sddat(:,:,izb) = 0.0_dp + Do k = 1, 3 + Do i = 1, 5 + sddat(i,k,izb) = 0.0_dp + EndDo + EndDo + + Do j = 0, max_order + Do i = 1, neq + z(i,j,izb) = 0.0_dp + z0(i,j,izb) = 0.0_dp + zt0(i,j,izb) = 0.0_dp + EndDo + EndDo - z(:,:,izb) = 0.0_dp - z(1:ny,0,izb) = yt(:,izb) - z(1:ny,1,izb) = ydot(:,izb) * tdel(izb) + Do i = 1, ny + z(i,0,izb) = yt(i,izb) + z(i,1,izb) = ydot(i,izb) * tdel(izb) + EndDo If ( iheat > 0 ) Then z(neq,0,izb) = t9t(izb) z(neq,1,izb) = t9dot(izb) * tdel(izb) EndIf - z0(:,:,izb) = 0.0_dp - zt0(:,:,izb) = 0.0_dp - hvec(:,izb) = tdel(izb) - lvec(:,izb) = 0.0_dp + Do i = 0, max_order + hvec(i,izb) = tdel(izb) + lvec(i,izb) = 0.0_dp + EndDo dt_scale(izb) = tdel(izb) - etaq(:,izb) = 0.0_dp + + Do i = -1, 1 + etaq(i,izb) = 0.0_dp + EndDo eta(izb) = 0.0_dp eta_next(izb) = eta_max gam(izb) = tdel(izb) gamhat(izb) = tdel(izb) gamratio(izb) = 1.0_dp - ewt(:,izb) = 0.0_dp - acor(:,izb) = 0.0_dp - acorp(:,izb) = 0.0_dp + Do i = 1, neq + ewt(i,izb) = 0.0_dp + acor(i,izb) = 0.0_dp + acorp(i,izb) = 0.0_dp + EndDo acnrm(izb) = 0.0_dp crate(izb) = 0.0_dp - tq(:,izb) = 0.0_dp + + Do i = -1, 2 + tq(i,izb) = 0.0_dp + EndDo tq2save(izb) = 0.0_dp - ydot0(:,izb) = ydot(:,izb) + Do i = 1, ny + ydot0(i,izb) = ydot(i,izb) + EndDo t9dot0(izb) = t9dot(izb) EndDo @@ -489,15 +516,21 @@ Subroutine bdf_predict(kstep) Integer, Intent(in) :: kstep ! Local variables - Integer :: i, j, izb, izone + Integer :: i, j, k, izb, izone Do izb = zb_lo, zb_hi If ( bdf_active(izb) ) Then - z0(:,:,izb) = z(:,:,izb) - zt0(:,:,izb) = z(:,:,izb) + Do j = 0, max_order + Do i = 1, neq + z0(i,j,izb) = z(i,j,izb) + zt0(i,j,izb) = z(i,j,izb) + EndDo + EndDo Do i = 0, q(izb) Do j = i+1, q(izb) - zt0(:,i,izb) = zt0(:,i,izb) + z(:,j,izb) * A_pascal(j,i) + Do k = 1, neq + zt0(k,i,izb) = zt0(k,i,izb) + z(k,j,izb) * A_pascal(j,i) + EndDo EndDo EndDo EndIf @@ -560,7 +593,9 @@ Subroutine bdf_update(kstep) ! Compute lvec and tq lvec(0,izb) = 1.0_dp lvec(1,izb) = 1.0_dp - lvec(2:max_order,izb) = 0.0_dp + Do i = 2, max_order + lvec(i,izb) = 0.0_dp + EndDo xi_inv = 1.0_dp xistar_inv = 1.0_dp alpha0 = -1.0_dp @@ -620,8 +655,12 @@ Subroutine bdf_update(kstep) gamhat(izb) = gam(izb) gamratio(izb) = 1.0_dp EndIf + EndIf + EndDo - If ( idiag >= 3 ) Then + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( bdf_active(izb) ) Then izone = izb + szbatch - zb_lo Write(lun_diag,"(a,2i5,i3,7es12.4)") 'BDF Update',kstep,izone,q(izb),tdel(izb),tt(izb) Write(lun_diag,"(2x,a5,6es12.4)") ' H',(hvec(i,izb),i=0,max_order) @@ -637,8 +676,8 @@ Subroutine bdf_update(kstep) EndDo Write(lun_diag,"(2x,a5,2(a5,es14.7))") 'EWT',(ewtname(i),ewt(iewt(i),izb),i=1,2) EndIf - EndIf - EndDo + EndDo + EndIf Return End Subroutine bdf_update @@ -672,19 +711,20 @@ Subroutine bdf_step(kstep) start_timer = xnet_wtime() timer_nraph = timer_nraph - start_timer - ! Load prediction into abundance vector + ! Load prediction into abundance vector and initialize masks Do izb = zb_lo, zb_hi If ( bdf_active(izb) ) Then - acor(:,izb) = 0.0_dp - yt(:,izb) = zt0(1:ny,0,izb) - If ( iheat > 0 ) t9t(izb) = zt0(neq,0,izb) + Do i = 1, ny + acor(i,izb) = 0.0_dp + yt(i,izb) = zt0(i,0,izb) + EndDo + If ( iheat > 0 ) Then + acor(neq,izb) = 0.0_dp + t9t(izb) = zt0(neq,0,izb) + EndIf + nit_nr(izb) = 0 + nit_nrslv(izb) = 0 EndIf - EndDo - If ( iheat == 0 ) Call t9rhofind(kstep,tt(zb_lo:zb_hi),ntt(zb_lo:zb_hi),t9t(zb_lo:zb_hi), & - rhot(zb_lo:zb_hi),mask_in = retry_ts(zb_lo:zb_hi)) - - ! Initialize masks - Do izb = zb_lo, zb_hi ! Determine which zones to rebuild the Newton matrix and factor refactor(izb) = ( bdf_active(izb) .and. & @@ -693,16 +733,14 @@ Subroutine bdf_step(kstep) ! Determine which zones' Jacobians need to be rebuilt rebuild(izb) = ( refactor(izb) .and. j_age(izb) > max_j_age ) - If ( bdf_active(izb) ) Then - nit_nr(izb) = 0 - nit_nrslv(izb) = 0 - EndIf iterate(izb) = bdf_active(izb) converged(izb) = .false. crate(izb) = 1.0_dp delp(izb) = 0.0_dp diag(izb) = 1.0_dp EndDo + If ( iheat == 0 ) Call t9rhofind(kstep,tt(zb_lo:zb_hi),ntt(zb_lo:zb_hi),t9t(zb_lo:zb_hi), & + rhot(zb_lo:zb_hi),mask_in = retry_ts(zb_lo:zb_hi)) ! Do the NR iteration Do kit = 1, max_it_nrslv @@ -731,6 +769,7 @@ Subroutine bdf_step(kstep) Call jacobian_build(mask_in = rebuild) Call jacobian_scale(diag,-gam(zb_lo:zb_hi),mask_in = refactor) Call jacobian_decomp(kstep,mask_in = refactor) + If ( idiag >= 4 ) Then Do izb = zb_lo, zb_hi izone = izb + szbatch - zb_lo @@ -740,7 +779,10 @@ Subroutine bdf_step(kstep) 'BDF Rebuild',kstep,izone,nit_nr(izb),gam(izb),abs(gamratio(izb)-1.0_dp),j_age(izb) EndDo EndIf + Do izb = zb_lo, zb_hi + + ! Reset Jacobian age If ( refactor(izb) ) Then gamhat(izb) = gam(izb) gamratio(izb) = 1.0_dp @@ -752,16 +794,17 @@ Subroutine bdf_step(kstep) j_age(izb) = 0 rebuild(izb) = .false. EndIf - EndDo - ! Calculate RHS - Do izb = zb_lo, zb_hi + ! Calculate RHS If ( iterate(izb) ) Then cstar = 2.0_dp / ( 1.0_dp + gamratio(izb) ) - yrhs(:,izb) = cstar * ( gam(izb)*ydot(:,izb) - zt0(1:ny,1,izb)/lvec(1,izb) - acor(1:ny,izb) ) + Do k = 1, ny + yrhs(k,izb) = cstar * ( gam(izb)*ydot(k,izb) - zt0(k,1,izb)/lvec(1,izb) - acor(k,izb) ) + EndDo If ( iheat > 0 ) t9rhs(izb) = cstar * ( gam(izb)*t9dot(izb) - zt0(neq,1,izb)/lvec(1,izb) - acor(neq,izb) ) EndIf EndDo + If ( idiag >= 4 ) Then Do izb = zb_lo, zb_hi If ( iterate(izb) ) Then @@ -784,16 +827,18 @@ Subroutine bdf_step(kstep) If ( iterate(izb) ) Then ! Prevent abundance from becoming too small - yt(:,izb) = zt0(1:ny,0,izb) + acor(1:ny,izb) + dy(:,izb) - Where ( yt(1:ny,izb) < ymin ) - yt(1:ny,izb) = 0.0_dp - dy(1:ny,izb) = -zt0(1:ny,0,izb) - acor(1:ny,izb) - ElseWhere ( yt(1:ny,izb) * aa(1:ny) > 1.0_dp ) - yt(1:ny,izb) = 1.0_dp / aa(1:ny) - dy(1:ny,izb) = 1.0_dp / aa(1:ny) - zt0(1:ny,0,izb) - acor(1:ny,izb) - EndWhere - acor(1:ny,izb) = acor(1:ny,izb) + dy(:,izb) - dvec(1:ny) = dy(:,izb) + Do k = 1, ny + yt(k,izb) = zt0(k,0,izb) + acor(k,izb) + dy(k,izb) + If ( yt(k,izb) < ymin ) Then + yt(k,izb) = 0.0_dp + dy(k,izb) = -zt0(k,0,izb) - acor(k,izb) + ElseIf ( yt(k,izb) * aa(k) > 1.0_dp ) Then + yt(k,izb) = 1.0_dp / aa(k) + dy(k,izb) = 1.0_dp / aa(k) - zt0(k,0,izb) - acor(k,izb) + EndIf + acor(k,izb) = acor(k,izb) + dy(k,izb) + dvec(k) = dy(k,izb) + EndDo If ( iheat > 0 ) Then dvec(neq) = dt9(izb) acor(neq,izb) = acor(neq,izb) + dt9(izb) @@ -804,6 +849,7 @@ Subroutine bdf_step(kstep) dcon(izb) = del(izb) * min(1.0_dp,crate(izb)) * tq(0,izb) EndIf EndDo + If ( idiag >= 3 ) Then Do izb = zb_lo, zb_hi If ( iterate(izb) ) Then @@ -847,12 +893,17 @@ Subroutine bdf_step(kstep) ! Reset the iteration and try again with a new Jacobian If ( j_age(izb) > 0 ) Then - yt(:,izb) = zt0(1:ny,0,izb) - If ( iheat > 0 ) t9t(izb) = zt0(neq,0,izb) + Do k = 1, ny + yt(k,izb) = zt0(k,0,izb) + acor(k,izb) = 0.0_dp + EndDo + If ( iheat > 0 ) Then + t9t(izb) = zt0(neq,0,izb) + acor(neq,izb) = 0.0_dp + EndIf refactor(izb) = .true. rebuild(izb) = ( abs(gamratio(izb) - 1.0_dp) < dgmaxj ) nit_nr(izb) = 0 - acor(:,izb) = 0.0_dp Else iterate(izb) = .false. nit_nr(izb) = max_it_nr + 1 @@ -927,19 +978,16 @@ Subroutine bdf_check(kstep) If ( nit_ts(izb) >= max_it_ts ) Then ierr_nr(izb) = BDF_STATUS_SKIP ierr_ts(izb) = BDF_STATUS_FAIL - Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Too many TS failures)' ! Solver failed to converge with minimum timestep, give up ElseIf ( tdel(izb) <= spacing(t(izb)) ) Then ierr_nr(izb) = BDF_STATUS_SKIP ierr_ts(izb) = BDF_STATUS_FAIL - Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Timestep too small)' ! Solver failed to converge too many times, give up ElseIf ( ncf_ts(izb) == max_ncf_ts ) Then ierr_nr(izb) = BDF_STATUS_SKIP ierr_ts(izb) = BDF_STATUS_FAIL - Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Too many NR failures)' ! Try again with smaller timestep Else @@ -978,20 +1026,15 @@ Subroutine bdf_check(kstep) ! Prevent stepsize from growing after this step eta_next(izb) = 1.0_dp - If ( idiag >= 2 ) Write(lun_diag,"(a,2i5,2i3,es12.4)") & - 'BDF TS Error',kstep,izone,nit_ts(izb),nerr_ts(izb),error - ! Solver failed error test too many times, give up If ( nerr_ts(izb) == max_nerr_ts ) Then ierr_nr(izb) = BDF_STATUS_SKIP ierr_ts(izb) = BDF_STATUS_FAIL - Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Too many failed error tests)' ! Solver failed error test with minimum timestep, give up ElseIf ( tdel(izb) <= spacing(t(izb)) ) Then ierr_nr(izb) = BDF_STATUS_SKIP ierr_ts(izb) = BDF_STATUS_FAIL - Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Timestep too small)' ! Solver failed error test after multiple attempts, reduce timestep, reduce order, and retry ElseIf ( nerr_ts(izb) > nerr_ts_qreduce ) Then @@ -1039,23 +1082,40 @@ Subroutine bdf_check(kstep) EndIf EndIf EndIf - EndIf - EndDo - ! Record iteration counters - Do izb = zb_lo, zb_hi - If ( bdf_active(izb) ) Then + ! Record iteration counters ktot(2,izb) = ktot(2,izb) + nit_nrslv(izb) nit_ts(izb) = nit_ts(izb) + 1 EndIf EndDo ! Log the failed integration attempts - If ( idiag >= 2 ) Then + If ( idiag >= 1 ) Then Do izb = zb_lo, zb_hi - If ( retry_ts(izb) ) Then + If ( bdf_active(izb) ) Then izone = izb + szbatch - zb_lo - Write(lun_diag,"(a,i5,3i3,3es12.4,i3)") & + If ( nit_nr(izb) > max_it_nr ) Then + If ( nit_ts(izb) >= max_it_ts ) Then + Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Too many TS failures)' + ElseIf ( tdel(izb) <= spacing(t(izb)) ) Then + Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Timestep too small)' + ElseIf ( ncf_ts(izb) == max_ncf_ts ) Then + Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Too many NR failures)' + EndIf + Else ! nit_nr(izb) <= max_it_nr + error = tq(0,izb) * acnrm(izb) + If ( error > 1.0_dp ) Then + If ( idiag >= 2 ) Write(lun_diag,"(a,2i5,2i3,es12.4)") & + 'BDF TS Error',kstep,izone,nit_ts(izb)-1,nerr_ts(izb),error + If ( nerr_ts(izb) == max_nerr_ts ) Then + Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Too many failed error tests)' + ElseIf ( tdel(izb) <= spacing(t(izb)) ) Then + Write(lun_stdout,"(a,2i5,1x,a)") 'BDF TS Fail',kstep,izone,'(Timestep too small)' + EndIf + EndIf + EndIf + + If ( idiag >= 2 .and. retry_ts(izb) ) Write(lun_diag,"(a,i5,3i3,3es12.4,i3)") & 'BDF TS Reduce',izone,nit_ts(izb),ncf_ts(izb),nerr_ts(izb),t(izb),tdel(izb),eta(izb),q(izb) EndIf EndDo @@ -1090,20 +1150,21 @@ Subroutine bdf_correct(kstep) Integer, Intent(in) :: kstep ! Local variables - Integer :: i, j, izb, izone + Integer :: i, j, k, izb, izone Do izb = zb_lo, zb_hi If ( ierr_nr(izb) == BDF_STATUS_SUCCESS ) Then ! Apply corrections to Nordsieck vector - z(1:ny,0,izb) = yt(:,izb) + Do k = 1, ny + z(k,0,izb) = yt(k,izb) + EndDo If ( iheat > 0 ) z(neq,0,izb) = t9t(izb) Do j = 1, q(izb) - z(:,j,izb) = zt0(:,j,izb) + acor(:,izb) * lvec(j,izb) + Do i = 1, neq + z(i,j,izb) = zt0(i,j,izb) + acor(i,izb) * lvec(j,izb) + EndDo EndDo - !Call dgemm('N','T',neq,q(izb)+1,1,1.0_dp,acor(:,izb),neq,lvec(:,izb),max_order+1,zt0(:,:,izb),neq) - !Call dger(neq,q(izb)+1,1.0_dp,acor(:,izb),1,lvec(:,izb),1,zt0(:,:,izb),neq) - !z(:,:,izb) = zt0(:,:,izb) ! Shift timestep history Do j = max_order, 1, -1 @@ -1153,38 +1214,38 @@ Subroutine bdf_eta(kstep) Real(dp), Parameter :: biasq = 6.0_dp Real(dp), Parameter :: biasqp1 = 10.0_dp Real(dp), Parameter :: addon = 1.0e-6_dp - Real(dp) :: c, error, acorhat(neq) + Real(dp) :: c, error(-1:1,zb_lo:zb_hi), acorhat(neq) Real(dp) :: alpha0, alpha1, prod, xi, xiold, hsum, a1 Integer :: i, j, izb, izone Do izb = zb_lo, zb_hi If ( ierr_nr(izb) == BDF_STATUS_SUCCESS ) Then izone = izb + szbatch - zb_lo + Do i = -1, 1 + error(i,izb) = 0.0_dp + EndDo ! Only allow change in order of there were no problems in the solve If ( eta_next(izb) > 1.0_dp ) Then ! Compute eta for each potential change in q etaq(:,izb) = 1.0_dp - error = tq(0,izb) * acnrm(izb) - etaq(0,izb) = 1.0_dp / ( (biasq * error) ** (1.0_dp/(q(izb)+1)) + addon ) - If ( idiag >= 3 ) Write(lun_diag,"(a,2es23.15)") 'BDF Eta(0)',error,etaq(0,izb) + error(0,izb) = tq(0,izb) * acnrm(izb) + etaq(0,izb) = 1.0_dp / ( (biasq * error(0,izb)) ** (1.0_dp/(q(izb)+1)) + addon ) If ( q_age(izb) > q(izb) ) Then ! eta for q-1 If ( q(izb) > 1 ) Then - error = tq(-1,izb) * normw( z(:,q(izb),izb), ewt(:,izb) ) - etaq(-1,izb) = 1.0_dp / ( (biasqm1 * error) ** (1.0_dp/q(izb)) + addon ) - If ( idiag >= 3 ) Write(lun_diag,"(a,2es23.15)") 'BDF Eta(-1)',error,etaq(-1,izb) + error(-1,izb) = tq(-1,izb) * normw( z(:,q(izb),izb), ewt(:,izb) ) + etaq(-1,izb) = 1.0_dp / ( (biasqm1 * error(-1,izb)) ** (1.0_dp/q(izb)) + addon ) EndIf ! eta for q+1 If ( q(izb) < max_order ) Then c = (tq(2,izb) / tq2save(izb)) * (hvec(1,izb) / hvec(2,izb)) ** (q(izb)+1) acorhat(:) = acor(:,izb) - c * acorp(:,izb) - error = tq(1,izb) * normw( acorhat(:), ewt(:,izb) ) - etaq(1,izb) = 1.0_dp / ( (biasqp1 * error) ** (1.0_dp/(q(izb)+2)) + addon ) - If ( idiag >= 3 ) Write(lun_diag,"(a,2es23.15)") 'BDF Eta(+1)',error,etaq(+1,izb) + error(1,izb) = tq(1,izb) * normw( acorhat(:), ewt(:,izb) ) + etaq(1,izb) = 1.0_dp / ( (biasqp1 * error(1,izb)) ** (1.0_dp/(q(izb)+2)) + addon ) EndIf ! Wait a bit before considering another order change @@ -1215,9 +1276,6 @@ Subroutine bdf_eta(kstep) ! Limit stepsize to maximum eta(izb) = min( eta(izb), eta_max ) - - If ( idiag >= 2 ) Write(lun_diag,"(a,2i5,2i3,sp,i3,ss,4es12.4)") & - 'BDF Eta',kstep,izone,q_age(izb),q(izb),deltaq(izb),eta(izb),(etaq(i,izb),i=-1,1) Else ! Defer order change if stepsize is stuck @@ -1230,6 +1288,24 @@ Subroutine bdf_eta(kstep) EndIf EndDo + If ( idiag >= 2 ) Then + Do izb = zb_lo, zb_hi + If ( ierr_nr(izb) == BDF_STATUS_SUCCESS ) Then + If ( eta_next(izb) > 1.0_dp ) Then + izone = izb + szbatch - zb_lo + If ( idiag >= 3 ) Then + Do i = -1, 1 + If ( abs(error(i,izb)) > 0.0_dp ) Write(lun_diag,"(a,sp,i2,s,a,2es23.15)") & + 'BDF Eta(',i,')',error(i,izb),etaq(i,izb) + EndDo + EndIf + If ( idiag >= 2 ) Write(lun_diag,"(a,2i5,2i3,sp,i3,s,4es12.4)") & + 'BDF Eta',kstep,izone,q_age(izb),q(izb),deltaq(izb),eta(izb),(etaq(i,izb),i=-1,1) + EndIf + EndIf + EndDo + EndIf + Return End Subroutine bdf_eta @@ -1284,7 +1360,7 @@ Subroutine bdf_stab(kstep) Do izb = zb_lo, zb_hi If ( ldflag(izb) > 3 ) Then izone = izb + szbatch - zb_lo - Write(lun_diag,"(a,2i5,i3,es12.4,sp,i3)") & + Write(lun_diag,"(a,2i5,i3,es12.4,sp,i3,s)") & 'BDF Stab',kstep,izone,q(izb),eta(izb),ldflag(izb) EndIf EndDo @@ -1324,11 +1400,19 @@ Subroutine bdf_stab_detect(mask,kflag) kflag = 0 loop0: Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - ssdat(:,:) = sddat(:,:,izb) + Do k = 1, 3 + Do i = 1, 5 + ssdat(i,k) = sddat(i,k,izb) + EndDo + EndDo rr = 0.0_dp Do k = 1, 3 - smink = minval(ssdat(:,k)) - smaxk = maxval(ssdat(:,k)) + smink = +huge(0.0_dp) + smaxk = -huge(0.0_dp) + Do i = 1, 5 + smink = min( smink, ssdat(i,k) ) + smaxk = max( smaxk, ssdat(i,k) ) + EndDo If ( smink < small*smaxk ) Then kflag(izb) = -1 Exit @@ -1350,17 +1434,30 @@ Subroutine bdf_stab_detect(mask,kflag) qc(3,k) = 0.0_dp qc(2,k) = ssdat(2,k)*ssdat(5,k) - ssdat(3,k)*ssdat(4,k) qc(1,k) = ssdat(4,k)*ssdat(4,k) - ssdat(3,k)*ssdat(5,k) - qco(:,k) = qc(:,k) + Do i = 1, 5 + qco(i,k) = qc(i,k) + EndDo + EndDo + vmin = +huge(0.0_dp) + vmax = -huge(0.0_dp) + Do i = 1, 4 + vmin = min( vmin, vrat(i) ) + vmax = max( vmax, vrat(i) ) EndDo - vmin = minval(vrat(:)) - vmax = maxval(vrat(:)) If ( vmin < vrrtol*vrrtol ) Then If ( vmax > vrrt2*vrrt2 ) Then kflag(izb) = -2 Cycle loop0 Else - rr = sum(rav(:)) / 3.0_dp - drrmax = maxval( abs(rav(:) - rr) ) + rr = 0.0_dp + Do i = 1, 3 + rr = rr + rav(i) + EndDo + rr = rr / 3.0_dp + drrmax = -huge(0.0_dp) + Do i = 1, 3 + drrmax = max( drrmax, abs(rav(i) - rr) ) + EndDo If ( drrmax > vrrt2 ) Then kflag(izb) = -3 Cycle loop0 @@ -1373,15 +1470,20 @@ Subroutine bdf_stab_detect(mask,kflag) kflag(izb) = -4 Cycle loop0 EndIf - qco(2:5,2) = qco(2:5,2) - (qco(1,2)/qco(1,1))*qco(2:5,1) - qco(1,2) = 0.0_dp - qco(2:5,3) = qco(2:5,3) - (qco(1,3)/qco(1,1))*qco(2:5,1) - qco(1,3) = 0.0_dp + Do k = 2, 3 + qco(1,k) = 0.0_dp + Do i = 2, 5 + qco(i,k) = qco(i,k) - (qco(1,k)/qco(1,1))*qco(i,1) + qco(i,k) = qco(i,k) - (qco(1,k)/qco(1,1))*qco(i,1) + EndDo + EndDo If ( abs(qco(2,2)) < small*ssmax(2) ) Then kflag(izb) = -4 Cycle loop0 EndIf - qco(3:5,3) = qco(3:5,3) - (qco(2,3)/qco(2,2))*qco(3:5,2) + Do i = 3, 5 + qco(i,3) = qco(i,3) - (qco(2,3)/qco(2,2))*qco(i,2) + EndDo If ( abs(qco(4,3)) < small*ssmax(3) ) Then kflag(izb) = -4 Cycle loop0 @@ -1391,31 +1493,44 @@ Subroutine bdf_stab_detect(mask,kflag) kflag(izb) = -5 Cycle loop0 EndIf - qkr(:) = qc(5,:) + rr*(qc(4,:) + rr*rr*(qc(2,:) + rr*qc(1,:))) - sqmax = maxval( abs(qkr(:))/ssmax(:) ) + sqmax = -huge(0.0_dp) + Do k = 1, 3 + qkr(k) = qc(5,k) + rr*(qc(4,k) + rr*rr*(qc(2,k) + rr*qc(1,k))) + sqmax = max( sqmax, abs(qkr(k))/ssmax(k) ) + EndDo If ( sqmax < sqtol ) Then kflag(izb) = 2 Else Do it = 1, 3 - qp(:) = qc(4,:) + rr*rr*(3.0_dp*qc(2,:) + 4.0_dp*rr*qc(1,:)) - Where ( abs(qp(:)) > small*ssmax(:) ) - rrc(:) = rr - qkr(:) / qp(:) - ElseWhere - rrc(:) = rr - EndWhere + sqmin = +huge(0.0_dp) + kmin = 3 Do k = 1, 3 + qp(k) = qc(4,k) + rr*rr*(3.0_dp*qc(2,k) + 4.0_dp*rr*qc(1,k)) + If ( abs(qp(k)) > small*ssmax(k) ) Then + rrc(k) = rr - qkr(k) / qp(k) + Else + rrc(k) = rr + EndIf s = rrc(k) - qjk(:,k) = qc(5,:) + s*(qc(4,:) + s*s*(qc(2,:) + s*qc(1,:))) - sqmx(k) = maxval( abs(qjk(:,k))/ssmax(:) ) + sqmax = -huge(0.0_dp) + Do j = 1, 3 + qjk(j,k) = qc(5,j) + s*(qc(4,j) + s*s*(qc(2,j) + s*qc(1,j))) + sqmax = max( sqmax, abs(qjk(j,k))/ssmax(j) ) + EndDo + sqmx(k) = sqmax + If ( sqmx(k) < sqmin ) Then + sqmin = sqmx(k) + kmin = k + EndIf EndDo - kmin = minloc(sqmx(:), dim=1) - sqmin = sqmx(kmin) rr = rrc(kmin) If ( sqmin < sqtol ) Then kflag(izb) = 3 Exit Else - qkr(:) = qjk(:,kmin) + Do j = 1, 3 + qkr(j) = qjk(j,kmin) + EndDo EndIf EndDo If ( sqmin > sqtol ) Then @@ -1491,12 +1606,21 @@ Subroutine bdf_adjust_order(kstep) Real(dp) :: alpha0, alpha1, prod, xi, xiold, hsum, a1 Integer :: i, j, izb, izone - Do izb = zb_lo, zb_hi + If ( idiag >= 2 ) Then + Do izb = zb_lo, zb_hi + izone = izb + szbatch - zb_lo + If ( deltaq(izb) /= 0 ) Write(lun_diag,"(a,2i5,i3,sp,i3,s)") & + 'BDF Order Change',kstep,izone,q(izb)+deltaq(izb),deltaq(izb) + EndDo + EndIf + Do izb = zb_lo, zb_hi If ( deltaq(izb) == -1 ) Then ! Decrease order - lvec(:,izb) = 0.0_dp + Do i = 0, max_order + lvec(i,izb) = 0.0_dp + EndDo lvec(2,izb) = 1.0_dp hsum = 0.0_dp Do j = 1, q(izb)-2 @@ -1507,19 +1631,24 @@ Subroutine bdf_adjust_order(kstep) EndDo EndDo Do j = 2, q(izb)-1 - z(:,j,izb) = z(:,j,izb) - z(:,q(izb),izb) * lvec(j,izb) + Do i = 1, neq + z(i,j,izb) = z(i,j,izb) - z(i,q(izb),izb) * lvec(j,izb) + EndDo + EndDo + Do i = 1, neq + z(i,q(izb),izb) = 0.0_dp EndDo - !Call dgemm('N','T',neq,q(izb)-2,1,-1.0_dp,z(:,q(izb),izb),neq,lvec(2,izb),max_order+1,z(1,2,izb),neq) - !Call dger(neq,q(izb)-2,-1.0_dp,z(:,q(izb),izb),1,lvec(2,izb),1,z(1,2,izb),neq) - z(:,q(izb),izb) = 0.0_dp q(izb) = q(izb) - 1 q_age(izb) = 0 ElseIf ( deltaq(izb) == +1 ) Then ! Increase order - lvec(:,izb) = 0.0_dp + Do i = 0, max_order + lvec(i,izb) = 0.0_dp + EndDo lvec(2,izb) = 1.0_dp + hsum = 0.0_dp alpha0 = -1.0_dp alpha1 = 1.0_dp prod = 1.0_dp @@ -1537,24 +1666,17 @@ Subroutine bdf_adjust_order(kstep) xiold = xi EndDo a1 = -(alpha0 + alpha1) / prod - z(:,q(izb)+1,izb) = a1 * acor(:,izb) Do j = 2, q(izb) - z(:,j,izb) = z(:,j,izb) + a1 * acor(:,izb) * lvec(j,izb) + Do i = 1, neq + z(i,j,izb) = z(i,j,izb) + a1 * acor(i,izb) * lvec(j,izb) + EndDo + EndDo + Do i = 1, neq + z(i,q(izb)+1,izb) = a1 * acor(i,izb) EndDo - !Call dgemm('N','T',neq,q(izb),1,a1,acor(:,izb),neq,lvec(2,izb),max_order+1,z(1,2,izb),neq) - !Call dger(neq,q(izb),a1,acor(:,izb),1,lvec(2,izb),1,z(1,2,izb),neq) q(izb) = q(izb) + 1 q_age(izb) = 0 EndIf - EndDo - If ( idiag >= 2 ) Then - Do izb = zb_lo, zb_hi - izone = izb + szbatch - zb_lo - If ( deltaq(izb) /= 0 ) Write(lun_diag,"(a,2i5,i3,sp,i3)") & - 'BDF Order Change',kstep,izone,q(izb),deltaq(izb) - EndDo - EndIf - Do izb = zb_lo, zb_hi deltaq(izb) = 0 EndDo @@ -1581,7 +1703,9 @@ Subroutine bdf_rescale(mask) tt(izb) = t(izb) + tdel(izb) ! Scale Nordsieck vector to account for change in stepsize Do j = 1, q(izb) - z(:,j,izb) = z(:,j,izb) * eta(izb)**j + Do i = 1, neq + z(i,j,izb) = z(i,j,izb) * eta(izb)**j + EndDo EndDo dt_scale(izb) = tdel(izb) nscon(izb) = 0 @@ -1598,14 +1722,18 @@ Subroutine bdf_restore(mask) Implicit None ! Input variables - Logical, Intent(in) :: mask(:) + Logical, Intent(in) :: mask(zb_lo:zb_hi) ! Local variables Integer :: i, j, izb Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - z(:,:,izb) = z0(:,:,izb) + Do j = 0, max_order + Do i = 1, neq + z(i,j,izb) = z0(i,j,izb) + EndDo + EndDo EndIf EndDo @@ -1627,7 +1755,9 @@ Subroutine pascal_build Ainv_pascal = 0 ! Set the first column - A_pascal(0:max_order,0) = 1 + Do i = 0, max_order + A_pascal(i,0) = 1 + EndDo ! Build the matrix Do i = 1, max_order @@ -1662,15 +1792,19 @@ Subroutine error_weights( y, wt ) n = size(y) If ( iconvc == 0 .or. iconvc == 1 ) Then - Where ( y(:) < ymin ) - wt(:) = 0.0_dp - ElseWhere ( y(:) < atol(:) ) - wt(:) = 1.0_dp / (rtol(:) * atol(:)) - ElseWhere - wt(:) = 1.0_dp / (rtol(:) * abs(y(:))) - EndWhere + Do i = 1, n + If ( y(i) < ymin ) Then + wt(i) = 0.0_dp + ElseIf ( y(i) < atol(i) ) Then + wt(i) = 1.0_dp / (rtol(i) * atol(i)) + Else + wt(i) = 1.0_dp / (rtol(i) * abs(y(i))) + EndIf + EndDo Else - wt(:) = 1.0_dp / ( rtol(:) * max(abs(y(:)),ymin) + atol(:) ) + Do i = 1, n + wt(i) = 1.0_dp / ( rtol(i) * max(abs(y(i)),ymin) + atol(i) ) + EndDo EndIf Return @@ -1689,27 +1823,34 @@ Function normw( x, wt ) Result( xnrm ) Implicit None ! Input variables - Real(dp), Intent(in) :: x(1:), wt(1:) + Real(dp), Intent(in) :: x(:), wt(:) ! Function variable Real(dp) :: xnrm ! Local variables - Real(dp) :: xw(size(x)) - Integer :: n, incx + Integer :: i, n, incx n = size(x) - xw(:) = x(:) * wt(:) + xnrm = 0.0_dp If ( iconvc == 0 ) Then - xnrm = maxval( abs(xw(:)) ) + Do i = 1, n + xnrm = max( xnrm, abs( x(i) * wt(i) ) ) + EndDo ElseIf ( iconvc == 1 ) Then - xnrm = sum( abs(xw(:)) ) + Do i = 1, n + xnrm = xnrm + abs( x(i) * wt(i) ) + EndDo ElseIf ( iconvc == 2 ) Then - xnrm = sqrt( sum(xw(:)**2) ) + Do i = 1, n + xnrm = xnrm + ( x(i) * wt(i) )**2 + EndDo + xnrm = sqrt( xnrm ) ElseIf ( iconvc == 3 ) Then - xnrm = sqrt( sum(xw(:)**2) / real(n,dp) ) - Else - xnrm = 0.0_dp + Do i = 1, n + xnrm = xnrm + ( x(i) * wt(i) )**2 + EndDo + xnrm = sqrt( xnrm / real(n,dp) ) EndIf Return diff --git a/source/xnet_integrate_be.f90 b/source/xnet_integrate_be.f90 index 6780bc37..7ff2bb98 100644 --- a/source/xnet_integrate_be.f90 +++ b/source/xnet_integrate_be.f90 @@ -20,9 +20,10 @@ Subroutine solve_be(kstep,its) ! the iteration, yt(i) is the trial abundance at the end of the timestep, ydot(i) is the time ! derivative of the trial abundance calculated from reaction rates and tdel is the timestep. !----------------------------------------------------------------------------------------------- - Use xnet_abundances, Only: y, yo, yt, xext - Use xnet_conditions, Only: t, to, tt, tdel, tdel_next, t9, t9o, t9t, rho, rhoo, rhot, & - yeo, ye, yet, nt, nto, ntt, t9rhofind, tdelstart + Use nuclear_data, Only: ny + Use xnet_abundances, Only: y, yo, yt + Use xnet_conditions, Only: t, to, tt, tdel, tdel_next, tdelstart, t9, t9o, t9t, rho, rhoo, & + & rhot, yeo, ye, yet, nt, nto, ntt, t9rhofind Use xnet_controls, Only: idiag, iheat, kitmx, kmon, ktot, lun_diag, lun_stdout, tdel_maxmult, & & szbatch, zb_lo, zb_hi Use xnet_integrate, Only: timestep @@ -38,7 +39,7 @@ Subroutine solve_be(kstep,its) ! Local variables Integer, Parameter :: ktsmx = 10 - Integer :: kts, izb, izone + Integer :: kts, k, izb, izone Integer :: inr(zb_lo:zb_hi) Integer :: mykts(zb_lo:zb_hi) Logical :: lzstep(zb_lo:zb_hi) @@ -55,6 +56,7 @@ Subroutine solve_be(kstep,its) inr(izb) = 0 lzstep(izb) = .true. EndIf + mykts(izb) = 0 EndDo !----------------------------------------------------------------------------------------------- @@ -62,7 +64,6 @@ Subroutine solve_be(kstep,its) ! The possible results for a timestep are a converged set of yt or a failure to converge. ! If convergence fails, retry with the trial timestep reduced by tdel_maxmult. !----------------------------------------------------------------------------------------------- - mykts = 0 Do kts = 1, ktsmx ! Attempt Backward Euler integration over desired timestep @@ -75,9 +76,12 @@ Subroutine solve_be(kstep,its) tdel(izb) = tdel(izb) / tdel_maxmult tt(izb) = t(izb) + tdel(izb) yet(izb) = ye(izb) - yt(:,izb) = y(:,izb) mykts(izb) = kts+1 + Do k = 1, ny + yt(k,izb) = y(k,izb) + EndDo + ! Record number of NR iterations If ( its(izb) == 0 ) Then kmon(2,izb) = kitmx + 1 @@ -94,15 +98,15 @@ Subroutine solve_be(kstep,its) ktot(2,izb) = ktot(2,izb) + inr(izb) EndIf EndIf + lzstep(izb) = ( inr(izb) == 0 ) EndDo ! Reset temperature and density for failed integrations - lzstep = ( inr == 0 ) Call t9rhofind(kstep,tt(zb_lo:zb_hi),ntt(zb_lo:zb_hi), & & t9t(zb_lo:zb_hi),rhot(zb_lo:zb_hi),mask_in = lzstep) If ( iheat > 0 ) Then Do izb = zb_lo, zb_hi - If ( inr(izb) == 0 ) Then + If ( lzstep(izb) ) Then t9t(izb) = t9(izb) EndIf EndDo @@ -112,7 +116,7 @@ Subroutine solve_be(kstep,its) ! as is done for the first timestep. If ( kts == ktsmx-1 ) Then Do izb = zb_lo, zb_hi - If ( inr(izb) == 0 ) Then + If ( lzstep(izb) ) Then tdel(izb) = 0.0 tdelstart(izb) = 0.0 EndIf @@ -135,46 +139,47 @@ Subroutine solve_be(kstep,its) EndDo ! Mark TS convergence only for zones which haven't previously failed or converged - Do izb = zb_lo, zb_hi - If ( inr(izb) > 0 ) Then - nto(izb) = nt(izb) - nt(izb) = ntt(izb) - to(izb) = t(izb) - t(izb) = tt(izb) - t9o(izb) = t9(izb) - t9(izb) = t9t(izb) - rhoo(izb) = rho(izb) - rho(izb) = rhot(izb) - yeo(izb) = ye(izb) - ye(izb) = yet(izb) - yo(:,izb) = y(:,izb) - y(:,izb) = yt(:,izb) - tdel_next(izb) = tdel(izb) * tdel_maxmult - ElseIf ( inr(izb) == 0 ) Then - its(izb) = 1 - EndIf - EndDo - - ! Log TS success/failure - Do izb = zb_lo, zb_hi - izone = izb + szbatch - zb_lo - If ( inr(izb) > 0 ) Then - If ( idiag >= 2 ) Write(lun_diag,"(a,2i5,2i3)") & - & 'BE TS Success',kstep,izone,mykts(izb),inr(izb) - ElseIf ( inr(izb) == 0 ) Then - If ( idiag >= 0 ) Write(lun_diag,"(a,2i5,4es12.4,2i3)") & - & 'BE TS Fail',kstep,izone,t(izb),tdel(izb),t9t(izb),rhot(izb),inr(izb),mykts(izb) - Write(lun_stdout,*) 'Timestep retrys fail after ',mykts(izb),' attempts' - EndIf - EndDo - Do izb = zb_lo, zb_hi If ( inr(izb) >= 0 ) Then kmon(1,izb) = mykts(izb) ktot(1,izb) = ktot(1,izb) + mykts(izb) + If ( inr(izb) > 0 ) Then + tdel_next(izb) = tdel(izb) * tdel_maxmult + nto(izb) = nt(izb) + nt(izb) = ntt(izb) + to(izb) = t(izb) + t(izb) = tt(izb) + t9o(izb) = t9(izb) + t9(izb) = t9t(izb) + rhoo(izb) = rho(izb) + rho(izb) = rhot(izb) + yeo(izb) = ye(izb) + ye(izb) = yet(izb) + Do k = 1, ny + yo(k,izb) = y(k,izb) + y(k,izb) = yt(k,izb) + EndDo + Else + its(izb) = 1 + EndIf EndIf EndDo + ! Log TS success/failure + If ( idiag >= 0 ) Then + Do izb = zb_lo, zb_hi + izone = izb + szbatch - zb_lo + If ( inr(izb) > 0 .and. idiag >= 2 ) Then + Write(lun_diag,"(a,2i5,2i3)") & + & 'BE TS Success',kstep,izone,mykts(izb),inr(izb) + ElseIf ( inr(izb) == 0 ) Then + Write(lun_diag,"(a,2i5,4es12.4,2i3)") & + & 'BE TS Fail',kstep,izone,t(izb),tdel(izb),t9t(izb),rhot(izb),inr(izb),mykts(izb) + Write(lun_stdout,*) 'Timestep retrys fail after ',mykts(izb),' attempts' + EndIf + EndDo + EndIf + stop_timer = xnet_wtime() timer_tstep = timer_tstep + stop_timer @@ -190,7 +195,7 @@ Subroutine step_be(kstep,inr) Use xnet_abundances, Only: y, ydot, yt, xext Use xnet_conditions, Only: cv, rhot, t9, t9dot, t9t, tdel, nh Use xnet_controls, Only: iconvc, idiag, iheat, ijac, kitmx, lun_diag, tolc, tolm, tolt9, ymin, & - & szbatch, zb_lo, zb_hi, iscrn + & szbatch, zb_lo, zb_hi Use xnet_integrate, Only: cross_sect, yderiv Use xnet_jacobian, Only: jacobian_bksub, jacobian_decomp, jacobian_build, jacobian_solve Use xnet_timers, Only: xnet_wtime, start_timer, stop_timer, timer_nraph @@ -208,45 +213,63 @@ Subroutine step_be(kstep,inr) ! Local variables Integer :: irdymx, idymx Integer :: i, k, kit, izb, izone - Real(dp) :: testc, testc2, testm, testn(zb_lo:zb_hi), toln(zb_lo:zb_hi) + Real(dp) :: s1, s2, s3, ytmp + Real(dp) :: testc(zb_lo:zb_hi), testc2(zb_lo:zb_hi), testm(zb_lo:zb_hi) + Real(dp) :: testn(zb_lo:zb_hi), toln(zb_lo:zb_hi) Real(dp) :: xtot(zb_lo:zb_hi), xtot_init(zb_lo:zb_hi), rdt(zb_lo:zb_hi), mult(zb_lo:zb_hi) - Real(dp) :: yrhs(ny,zb_lo:zb_hi), dy(ny,zb_lo:zb_hi), reldy(ny) - Real(dp) :: t9rhs(zb_lo:zb_hi), dt9(zb_lo:zb_hi), relt9 + Real(dp) :: yrhs(ny,zb_lo:zb_hi), dy(ny,zb_lo:zb_hi), reldy(ny,zb_lo:zb_hi) + Real(dp) :: t9rhs(zb_lo:zb_hi), dt9(zb_lo:zb_hi), relt9(zb_lo:zb_hi) Logical :: iterate(zb_lo:zb_hi), eval_rates(zb_lo:zb_hi), rebuild(zb_lo:zb_hi) start_timer = xnet_wtime() timer_nraph = timer_nraph - start_timer - ! Create mask for active zone integrations - iterate = ( inr == 0 ) - ! Calculate initial total mass fraction Do izb = zb_lo, zb_hi - If ( iterate(izb) ) Then - xtot_init(izb) = sum(aa*y(:,izb)) + xext(izb) - 1.0 + If ( inr(izb) == 0 ) Then + + ! Create mask for active zone integrations + iterate(izb) = .true. + + ! Calculate initial total mass fraction + s1 = 0.0 + Do k = 1, ny + s1 = s1 + aa(k)*y(k,izb) + EndDo + xtot_init(izb) = s1 + xext(izb) - 1.0 + + ! Jacobian scaling factors rdt(izb) = 1.0 / tdel(izb) mult(izb) = -1.0 Else + iterate(izb) = .false. xtot_init(izb) = 0.0 rdt(izb) = 0.0 mult(izb) = 0.0 EndIf + If ( iconvc /= 0 ) Then + toln(izb) = tolc + Else + toln(izb) = tolm + EndIf EndDo ! The Newton-Raphson iteration occurs for at most kitmx iterations. Do kit = 1, kitmx + Do izb = zb_lo, zb_hi - ! Rebuild and update LU factorization of the jacobian every ijac iterations - rebuild = ( iterate .and. mod(kit-1,ijac) == 0 ) + ! Rebuild and update LU factorization of the jacobian every ijac iterations + rebuild(izb) = ( iterate(izb) .and. mod(kit-1,ijac) == 0 ) - ! If thermodynamic conditions are changing, rates need to be udpated - If ( iheat > 0 ) Then - eval_rates = iterate - ElseIf ( kit == 1 ) Then - eval_rates = ( iterate .and. nh(zb_lo:zb_hi) > 1 ) - Else - eval_rates = .false. - EndIf + ! If thermodynamic conditions are changing, rates need to be udpated + If ( iheat > 0 ) Then + eval_rates(izb) = iterate(izb) + ElseIf ( kit == 1 ) Then + eval_rates(izb) = ( iterate(izb) .and. nh(izb) > 1 ) + Else + eval_rates(izb) = .false. + EndIf + EndDo ! Calculate the reaction rates and abundance time derivatives Call cross_sect(mask_in = eval_rates) @@ -257,8 +280,22 @@ Subroutine step_be(kstep,inr) ! Calculate equation to zero Do izb = zb_lo, zb_hi If ( iterate(izb) ) Then - yrhs(:,izb) = (y(:,izb)-yt(:,izb))*rdt(izb) + ydot(:,izb) - If ( iheat > 0 ) t9rhs(izb) = (t9(izb)*rdt(izb)-t9t(izb)*rdt(izb)) + t9dot(izb) + If ( kit > 1 ) Then + Do k = 1, ny + yrhs(k,izb) = (y(k,izb)-yt(k,izb))*rdt(izb) + ydot(k,izb) + EndDo + If ( iheat > 0 ) t9rhs(izb) = (t9(izb)-t9t(izb))*rdt(izb) + t9dot(izb) + Else + Do k = 1, ny + yrhs(k,izb) = ydot(k,izb) + EndDo + If ( iheat > 0 ) t9rhs(izb) = t9dot(izb) + EndIf + Else + Do k = 1, ny + yrhs(k,izb) = 0.0 + EndDo + If ( iheat > 0 ) t9rhs(izb) = 0.0 EndIf EndDo If ( idiag >= 4 ) Then @@ -274,95 +311,113 @@ Subroutine step_be(kstep,inr) ! Solve the jacobian and calculate the changes in abundances, dy !Call jacobian_solve(kstep,yrhs,dy,t9rhs,dt9,mask_in = iterate) - !Call jacobian_decomp(kstep,mask_in = iterate) Call jacobian_bksub(kstep,yrhs,dy,t9rhs,dt9,mask_in = iterate) ! Evolve the abundances and calculate convergence tests + !----------------------------------------------------------------------------------------- + ! There are 3 included convergence tests: + ! testc which measures relative changes, + ! testc2 which measures total abundance changes, and + ! testm which tests mass conservation. + ! + ! testc is the most stringent test, and requires the most iterations. + ! testm is the most lax, and therefore the fastest, often requiring only one iteration. + ! Considering the uncertainties of the reaction rates, it is doubtful that the + ! increased precision of testc is truly increased accuracy. + !----------------------------------------------------------------------------------------- Do izb = zb_lo, zb_hi If ( iterate(izb) ) Then + s1 = 0.0 + s2 = 0.0 + s3 = 0.0 Do k = 1, ny - yt(k,izb) = yt(k,izb) + dy(k,izb) - If ( yt(k,izb) < ymin ) Then + ytmp = yt(k,izb) + dy(k,izb) + If ( ytmp < ymin ) Then + dy(k,izb) = -yt(k,izb) yt(k,izb) = 0.0 - reldy(k) = 0.0 + reldy(k,izb) = 0.0 Else - reldy(k) = abs(dy(k,izb) / yt(k,izb)) + yt(k,izb) = ytmp + reldy(k,izb) = abs(dy(k,izb) / yt(k,izb)) EndIf + s1 = s1 + reldy(k,izb) + s2 = s2 + aa(k)*dy(k,izb) + s3 = s3 + aa(k)*yt(k,izb) EndDo + testc(izb) = s1 + testc2(izb) = s2 + xtot(izb) = s3 + xext(izb) - 1.0 + testm(izb) = xtot(izb) - xtot_init(izb) + If ( iheat > 0 ) Then t9t(izb) = t9t(izb) + dt9(izb) If ( abs(t9t(izb)) > tiny(0.0) ) Then - relt9 = abs(dt9(izb) / t9t(izb)) + relt9(izb) = abs(dt9(izb) / t9t(izb)) Else - relt9 = 0.0 + relt9(izb) = 0.0 EndIf Else - relt9 = 0.0 - EndIf - If ( idiag >= 4 ) Then - izone = izb + szbatch - zb_lo - irdymx = maxloc(reldy,dim=1) - idymx = maxloc(dy(:,izb),dim=1) - Write(lun_diag,"(a3,2i5,i3,2(a5,2es23.15))") & - & ' dY',kstep,izone,kit,nname(idymx),dy(idymx,izb),y(idymx,izb),nname(irdymx),reldy(irdymx),y(irdymx,izb) - If ( idiag >= 5 ) Write(lun_diag,"(a5,5es12.4)") & - & (nname(k),yt(k,izb),dy(k,izb),reldy(k),(aa(k)*dy(k,izb)),(aa(k)*yt(k,izb)),k=1,ny) - If ( iheat > 0 ) Write(lun_diag,"(a3,2i5,i3,5x,2es23.15,5x,es12.4)") & - & 'dT9',kstep,izone,kit,dt9(izb),t9t(izb),relt9 + relt9(izb) = 0.0 EndIf - !----------------------------------------------------------------------------------------- - ! There are 3 included convergence tests: - ! testc which measures relative changes, - ! testc2 which measures total abundance changes, and - ! testm which tests mass conservation. - !----------------------------------------------------------------------------------------- - testc = sum(reldy) - testc2 = sum(aa*dy(:,izb)) - xtot(izb) = sum(aa*yt(:,izb)) + xext(izb) - 1.0 - testm = xtot(izb) - xtot_init(izb) - If ( idiag >= 3 ) Write(lun_diag,"(a,3i5,3es14.6)") 'NR',kstep,izone,kit,testm,testc,testc2 - - !----------------------------------------------------------------------------------------- - ! testc is the most stringent test, and requires the most iterations. - ! testm is the most lax, and therefore the fastest, often requiring only one iteration. - ! Considering the uncertainties of the reaction rates, it is doubtful that the - ! increased precision of testc is truly increased accuracy. - !----------------------------------------------------------------------------------------- ! Ordinarily, test for true convergence If ( iconvc /= 0 ) Then - testn(izb) = testc - toln(izb) = tolc + testn(izb) = testc(izb) ! Otherwise, use mass conservation for convergence condition Else - testn(izb) = testm - toln(izb) = tolm + testn(izb) = testm(izb) EndIf ! If converged, exit NR loop - If ( abs(testn(izb)) <= toln(izb) .and. relt9 < tolt9 ) Then + If ( abs(testn(izb)) <= toln(izb) .and. relt9(izb) < tolt9 ) Then inr(izb) = kit iterate(izb) = .false. + Else + iterate(izb) = .true. EndIf EndIf EndDo + If ( idiag >= 3 ) Then + Do izb = zb_lo, zb_hi + If ( inr(izb) >= 0 ) Then + izone = izb + szbatch - zb_lo + If ( idiag >= 4 ) Then + irdymx = maxloc(reldy(:,izb),dim=1) + idymx = maxloc(dy(:,izb),dim=1) + Write(lun_diag,"(a3,2i5,i3,2(a5,2es23.15))") & + & ' dY',kstep,izone,kit,nname(idymx),dy(idymx,izb),y(idymx,izb), & + & nname(irdymx),reldy(irdymx,izb),y(irdymx,izb) + If ( idiag >= 5 ) Write(lun_diag,"(a5,2es23.15,es12.4,2es23.15)") & + & (nname(k),yt(k,izb),dy(k,izb),reldy(k,izb), & + & (aa(k)*dy(k,izb)),(aa(k)*yt(k,izb)),k=1,ny) + If ( iheat > 0 ) Write(lun_diag,"(a3,2i5,i3,5x,2es23.15,5x,es12.4)") & + & 'dT9',kstep,izone,kit,dt9(izb),t9t(izb),relt9(izb) + EndIf + Write(lun_diag,"(a,3i5,3es14.6)") 'NR',kstep,izone,kit,testm(izb),testc(izb),testc2(izb) + EndIf + EndDo + EndIf + ! Check that all zones are converged - If ( .not. any(iterate(:)) ) Exit + If ( .not. any(iterate) ) Exit EndDo If ( idiag >= 2 ) Then Do izb = zb_lo, zb_hi - izone = izb + szbatch - zb_lo - If ( inr(izb) > 0 ) Then - Write(lun_diag,"(a,3i5,3es12.4)") 'Conv',kstep,izone,inr(izb),xtot(izb),testn(izb),toln(izb) - ElseIf ( inr(izb) == 0 ) Then - Write(lun_diag,"(a,3i5,2es10.2)") 'BE Failure',izone,inr(izb),kitmx,xtot(izb),testn(izb) + If ( inr(izb) >= 0 ) Then + izone = izb + szbatch - zb_lo + If ( inr(izb) > 0 ) Then + Write(lun_diag,"(a,3i5,3es12.4)") 'Conv',kstep,izone,inr(izb),xtot(izb),testn(izb),toln(izb) + ElseIf ( inr(izb) == 0 ) Then + Write(lun_diag,"(a,3i5,2es10.2)") 'BE Failure',izone,inr(izb),kitmx,xtot(izb),testn(izb) + EndIf EndIf EndDo EndIf + stop_timer = xnet_wtime() timer_nraph = timer_nraph + stop_timer diff --git a/source/xnet_jacobian_dense.f90 b/source/xnet_jacobian_dense.f90 index eece4a6f..57f3db25 100644 --- a/source/xnet_jacobian_dense.f90 +++ b/source/xnet_jacobian_dense.f90 @@ -17,7 +17,9 @@ Module xnet_jacobian Implicit None Real(dp), Allocatable :: dydotdy(:,:,:) ! dYdot/dY part of jac Real(dp), Allocatable :: jac(:,:,:) ! the Jacobian matrix + Real(dp), Allocatable :: rhs(:,:) ! the right-hand-side of linear system Integer, Allocatable :: indx(:,:) ! Pivots in the LU decomposition + Integer, Allocatable :: info(:) Integer :: msize ! Size of linear system to be solved Real(dp), Allocatable, Target :: diag0(:) @@ -30,24 +32,32 @@ Subroutine read_jacobian_data(data_dir) ! Initializes the Jacobian data. !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: ny - Use xnet_controls, Only: iheat, nzevolve, nzbatchmx + Use xnet_controls, Only: iheat, nzevolve Implicit None ! Input variables Character(*), Intent(in) :: data_dir + ! Calculate array sizes If ( iheat > 0 ) Then msize = ny + 1 Else msize = ny EndIf - Allocate (dydotdy(msize,msize,nzevolve),jac(msize,msize,nzevolve),indx(msize,nzevolve)) + Allocate (dydotdy(msize,msize,nzevolve)) + Allocate (jac(msize,msize,nzevolve)) + Allocate (rhs(msize,nzevolve)) + Allocate (indx(msize,nzevolve)) + Allocate (info(nzevolve)) dydotdy = 0.0 jac = 0.0 + rhs = 0.0 indx = 0 + info = 0 - Allocate (diag0(nzbatchmx),mult1(nzbatchmx)) + Allocate (diag0(nzevolve)) + Allocate (mult1(nzevolve)) diag0 = 0.0 mult1 = 1.0 @@ -56,12 +66,12 @@ End Subroutine read_jacobian_data Subroutine jacobian_finalize Implicit None - Deallocate (dydotdy,jac,indx) + Deallocate (dydotdy,jac,rhs,indx,info) Deallocate (diag0,mult1) Return End Subroutine jacobian_finalize - Subroutine jacobian_scale(diag,mult,mask_in) + Subroutine jacobian_scale(diag_in,mult_in,mask_in) !----------------------------------------------------------------------------------------------- ! This augments a previously calculation Jacobian matrix by multiplying all elements by mult and ! adding diag to the diagonal elements. @@ -72,13 +82,15 @@ Subroutine jacobian_scale(diag,mult,mask_in) Implicit None ! Input variables - Real(dp), Intent(in) :: diag(zb_lo:zb_hi), mult(zb_lo:zb_hi) + Real(dp), Optional, Target, Intent(in) :: diag_in(zb_lo:zb_hi), mult_in(zb_lo:zb_hi) ! Optional variables Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Integer :: i, j, i0, izb, izone + Integer :: i, j, i0, j1, izb, izone + Real(dp) :: rsum + Real(dp), Pointer :: diag(:), mult(:) Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -88,13 +100,29 @@ Subroutine jacobian_scale(diag,mult,mask_in) EndIf If ( .not. any(mask) ) Return + If ( present(diag_in) ) Then + diag(zb_lo:) => diag_in + Else + diag(zb_lo:) => diag0(zb_lo:zb_hi) + EndIf + If ( present(mult_in) ) Then + mult(zb_lo:) => mult_in + Else + mult(zb_lo:) => mult1(zb_lo:zb_hi) + EndIf + Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - jac(:,:,izb) = mult(izb) * dydotdy(:,:,izb) - Do i0 = 1, msize - jac(i0,i0,izb) = jac(i0,i0,izb) + diag(izb) - EndDo - EndIf + Do j1 = 1, msize + If ( mask(izb) ) Then + Do i0 = 1, msize + jac(i0,j1,izb) = mult(izb) * dydotdy(j1,i0,izb) + If ( i0 == j1 ) Then + jac(i0,j1,izb) = jac(i0,j1,izb) + diag(izb) + Else + EndIf + EndDo + EndIf + EndDo EndDo If ( idiag >= 5 ) Then @@ -128,7 +156,8 @@ Subroutine jacobian_build(diag_in,mult_in,mask_in) !----------------------------------------------------------------------------------------------- Use nuclear_data, Only: ny, mex, nname Use reaction_data, Only: a1, a2, a3, a4, b1, b2, b3, b4, la, le, mu1, mu2, mu3, mu4, n11, n21, & - & n22, n31, n32, n33, n41, n42, n43, n44, dcsect1dt9, dcsect2dt9, dcsect3dt9, dcsect4dt9 + & n22, n31, n32, n33, n41, n42, n43, n44, dcsect1dt9, dcsect2dt9, dcsect3dt9, dcsect4dt9, nan, & + & n10, n20, n30, n40 Use xnet_abundances, Only: yt Use xnet_conditions, Only: cv Use xnet_controls, Only: iheat, idiag, ktot, lun_diag, nzbatchmx, szbatch, zb_lo, zb_hi, lzactive @@ -142,12 +171,7 @@ Subroutine jacobian_build(diag_in,mult_in,mask_in) ! Local variables Integer :: i, j, i0, j1, izb, izone - Integer :: la1, le1, la2, le2, la3, le3, la4, le4 - Integer :: i11, i21, i22, i31, i32, i33, i41, i42, i43, i44 - Real(dp) :: s1, s2, s3, s4, r1, r2, r3, r4 - Real(dp) :: y11, y21, y22, y31, y32, y33, y41, y42, y43, y44 - Real(dp) :: dt9dotdy(msize), dr1dt9, dr2dt9, dr3dt9, dr4dt9 - Real(dp), Pointer :: diag(:), mult(:) + Real(dp) :: s1, s2, s3, s4, sdot Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -160,155 +184,77 @@ Subroutine jacobian_build(diag_in,mult_in,mask_in) start_timer = xnet_wtime() timer_jacob = timer_jacob - start_timer - If ( present(diag_in) ) Then - diag(zb_lo:) => diag_in - Else - diag(zb_lo:) => diag0 - EndIf - If ( present(mult_in) ) Then - mult(zb_lo:) => mult_in - Else - mult(zb_lo:) => mult1 - EndIf - ! Build the Jacobian Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - dydotdy(:,:,izb) = 0.0 - Do i0 = 1, ny - la1 = la(1,i0) - la2 = la(2,i0) - la3 = la(3,i0) - la4 = la(4,i0) - le1 = le(1,i0) - le2 = le(2,i0) - le3 = le(3,i0) - le4 = le(4,i0) - Do j1 = la1, le1 - r1 = b1(j1,izb) - i11 = n11(j1) - dydotdy(i0,i11,izb) = dydotdy(i0,i11,izb) + r1 + Do i0 = 1, ny + If ( mask(izb) ) Then + Do j1 = 1, msize + dydotdy(j1,i0,izb) = 0.0 EndDo - Do j1 = la2, le2 - r2 = b2(j1,izb) - i21 = n21(j1) - i22 = n22(j1) - y21 = yt(i21,izb) - y22 = yt(i22,izb) - dydotdy(i0,i21,izb) = dydotdy(i0,i21,izb) + r2 * y22 - dydotdy(i0,i22,izb) = dydotdy(i0,i22,izb) + r2 * y21 + Do j1 = la(1,i0), le(1,i0) + dydotdy(n11(j1),i0,izb) = dydotdy(n11(j1),i0,izb) + b1(j1,izb) EndDo - Do j1 = la3, le3 - r3 = b3(j1,izb) - i31 = n31(j1) - i32 = n32(j1) - i33 = n33(j1) - y31 = yt(i31,izb) - y32 = yt(i32,izb) - y33 = yt(i33,izb) - dydotdy(i0,i31,izb) = dydotdy(i0,i31,izb) + r3 * y32 * y33 - dydotdy(i0,i32,izb) = dydotdy(i0,i32,izb) + r3 * y33 * y31 - dydotdy(i0,i33,izb) = dydotdy(i0,i33,izb) + r3 * y31 * y32 + Do j1 = la(2,i0), le(2,i0) + dydotdy(n21(j1),i0,izb) = dydotdy(n21(j1),i0,izb) + b2(j1,izb) * yt(n22(j1),izb) + dydotdy(n22(j1),i0,izb) = dydotdy(n22(j1),i0,izb) + b2(j1,izb) * yt(n21(j1),izb) EndDo - Do j1 = la4, le4 - r4 = b4(j1,izb) - i41 = n41(j1) - i42 = n42(j1) - i43 = n43(j1) - i44 = n44(j1) - y41 = yt(i41,izb) - y42 = yt(i42,izb) - y43 = yt(i43,izb) - y44 = yt(i44,izb) - dydotdy(i0,i41,izb) = dydotdy(i0,i41,izb) + r4 * y42 * y43 * y44 - dydotdy(i0,i42,izb) = dydotdy(i0,i42,izb) + r4 * y43 * y44 * y41 - dydotdy(i0,i43,izb) = dydotdy(i0,i43,izb) + r4 * y44 * y41 * y42 - dydotdy(i0,i44,izb) = dydotdy(i0,i44,izb) + r4 * y41 * y42 * y43 + Do j1 = la(3,i0), le(3,i0) + dydotdy(n31(j1),i0,izb) = dydotdy(n31(j1),i0,izb) + b3(j1,izb) * yt(n32(j1),izb) * yt(n33(j1),izb) + dydotdy(n32(j1),i0,izb) = dydotdy(n32(j1),i0,izb) + b3(j1,izb) * yt(n33(j1),izb) * yt(n31(j1),izb) + dydotdy(n33(j1),i0,izb) = dydotdy(n33(j1),i0,izb) + b3(j1,izb) * yt(n31(j1),izb) * yt(n32(j1),izb) + EndDo + Do j1 = la(4,i0), le(4,i0) + dydotdy(n41(j1),i0,izb) = dydotdy(n41(j1),i0,izb) + b4(j1,izb) * yt(n42(j1),izb) * yt(n43(j1),izb) * yt(n44(j1),izb) + dydotdy(n42(j1),i0,izb) = dydotdy(n42(j1),i0,izb) + b4(j1,izb) * yt(n43(j1),izb) * yt(n44(j1),izb) * yt(n41(j1),izb) + dydotdy(n43(j1),i0,izb) = dydotdy(n43(j1),i0,izb) + b4(j1,izb) * yt(n44(j1),izb) * yt(n41(j1),izb) * yt(n42(j1),izb) + dydotdy(n44(j1),i0,izb) = dydotdy(n44(j1),i0,izb) + b4(j1,izb) * yt(n41(j1),izb) * yt(n42(j1),izb) * yt(n43(j1),izb) EndDo - EndDo - EndIf - EndDo - If ( iheat > 0 ) Then - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - Do i0 = 1, ny - la1 = la(1,i0) - la2 = la(2,i0) - la3 = la(3,i0) - la4 = la(4,i0) - le1 = le(1,i0) - le2 = le(2,i0) - le3 = le(3,i0) - le4 = le(4,i0) + If ( iheat > 0 ) Then s1 = 0.0 - Do j1 = la1, le1 - dr1dt9 = a1(j1)*dcsect1dt9(mu1(j1),izb) - i11 = n11(j1) - y11 = yt(i11,izb) - s1 = s1 + dr1dt9 * y11 + Do j1 = la(1,i0), le(1,i0) + s1 = s1 + a1(j1) * dcsect1dt9(mu1(j1),izb) * yt(n11(j1),izb) EndDo s2 = 0.0 - Do j1 = la2, le2 - dr2dt9 = a2(j1)*dcsect2dt9(mu2(j1),izb) - i21 = n21(j1) - i22 = n22(j1) - y21 = yt(i21,izb) - y22 = yt(i22,izb) - s2 = s2 + dr2dt9 * y21 * y22 + Do j1 = la(2,i0), le(2,i0) + s2 = s2 + a2(j1) * dcsect2dt9(mu2(j1),izb) * yt(n21(j1),izb) * yt(n22(j1),izb) EndDo s3 = 0.0 - Do j1 = la3, le3 - dr3dt9 = a3(j1)*dcsect3dt9(mu3(j1),izb) - i31 = n31(j1) - i32 = n32(j1) - i33 = n33(j1) - y31 = yt(i31,izb) - y32 = yt(i32,izb) - y33 = yt(i33,izb) - s3 = s3 + dr3dt9 * y31 * y32 * y33 + Do j1 = la(3,i0), le(3,i0) + s3 = s3 + a3(j1) * dcsect3dt9(mu3(j1),izb) * yt(n31(j1),izb) * yt(n32(j1),izb) * yt(n33(j1),izb) EndDo s4 = 0.0 - Do j1 = la4, le4 - dr4dt9 = a4(j1)*dcsect4dt9(mu4(j1),izb) - i41 = n41(j1) - i42 = n42(j1) - i43 = n43(j1) - i44 = n44(j1) - y41 = yt(i41,izb) - y42 = yt(i42,izb) - y43 = yt(i43,izb) - y44 = yt(i44,izb) - s4 = s4 + dr4dt9 * y41 * y42 * y43 * y44 + Do j1 = la(4,i0), le(4,i0) + s4 = s4 + a4(j1) * dcsect4dt9(mu4(j1),izb) * yt(n41(j1),izb) * yt(n42(j1),izb) * yt(n43(j1),izb) * yt(n44(j1),izb) EndDo - dydotdy(i0,ny+1,izb)= s1 + s2 + s3 + s4 - EndDo - - ! The BLAS version of dydotdy(ny+1,i,izb) = -sum(mex*dydotdy(1:ny,i,izb))/cv(izb) - Call dgemv('T',ny,msize,1.0/cv(izb),dydotdy(1,1,izb),msize,mex,1,0.0,dt9dotdy,1) - dydotdy(ny+1,:,izb) = -dt9dotdy + dydotdy(ny+1,i0,izb) = s1 + s2 + s3 + s4 + EndIf EndIf EndDo + EndDo - ! This also works, but could be inefficient if there are few active zones in batch - !Call dgemv('T',ny,msize*nzbatchmx,1.0,dydotdy(1,1,1),msize,mex,1,0.0,dt9dotdy,1) - !ForAll ( izb = 1:nzbatchmx, j1 = 1:msize, mask(izb) ) - ! dydotdy(ny+1,j1,izb) = -dt9dotdy(j1,izb) / cv(izb) - !EndForAll + If ( iheat > 0 ) Then + Do izb = zb_lo, zb_hi + Do j1 = 1, msize + If ( mask(izb) ) Then + sdot = 0.0 + Do i0 = 1, ny + sdot = sdot - mex(i0)*dydotdy(j1,i0,izb) / cv(izb) + EndDo + dydotdy(j1,ny+1,izb) = sdot + EndIf + EndDo + EndDo EndIf + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + ktot(3,izb) = ktot(3,izb) + 1 + EndIf + EndDo + ! Apply the externally provided factors - If ( present(diag_in) ) Then - diag = diag_in - Else - diag = 0.0 - EndIf - If ( present(mult_in) ) Then - mult = mult_in - Else - mult = 1.0 - EndIf - Call jacobian_scale(diag,mult,mask_in = mask) + Call jacobian_scale(diag_in,mult_in,mask_in = mask) If ( idiag >= 5 ) Then Do izb = zb_lo, zb_hi @@ -331,22 +277,6 @@ Subroutine jacobian_build(diag_in,mult_in,mask_in) EndDo EndIf - If ( idiag >= 6 ) Then - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - izone = izb + szbatch - zb_lo - Write(lun_diag,"(a9,i5,2es14.7)") 'JAC_BUILD',izone,diag(izb),mult(izb) - Write(lun_diag,"(14es9.1)") ((jac(i,j,izb),j=1,msize),i=1,msize) - EndIf - EndDo - EndIf - - Do izb = zb_lo, zb_hi - If ( mask(izb) ) Then - ktot(3,izb) = ktot(3,izb) + 1 - EndIf - EndDo - stop_timer = xnet_wtime() timer_jacob = timer_jacob + stop_timer @@ -376,8 +306,7 @@ Subroutine jacobian_solve(kstep,yrhs,dy,t9rhs,dt9,mask_in) Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Real(dp) :: rhs(msize) - Integer :: i, izb, izone, info(zb_lo:zb_hi) + Integer :: i, izb, izone Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -390,14 +319,32 @@ Subroutine jacobian_solve(kstep,yrhs,dy,t9rhs,dt9,mask_in) start_timer = xnet_wtime() timer_solve = timer_solve - start_timer + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Do i = 1, ny + rhs(i,izb) = yrhs(i,izb) + EndDo + If ( iheat > 0 ) rhs(ny+1,izb) = t9rhs(izb) + Else + Do i = 1, msize + rhs(i,izb) = 0.0 + EndDo + EndIf + EndDo + ! Solve the linear system Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - rhs(1:ny) = yrhs(:,izb) - If ( iheat > 0 ) rhs(ny+1) = t9rhs(izb) - Call dgesv(msize,1,jac(1,1,izb),msize,indx(1,izb),rhs,msize,info(izb)) - dy(:,izb) = rhs(1:ny) - If ( iheat > 0 ) dt9(izb) = rhs(ny+1) + Call dgesv(msize,1,jac(1,1,izb),msize,indx(1,izb),rhs(1,izb),msize,info(izb)) + EndIf + EndDo + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Do i = 1, ny + dy(i,izb) = rhs(i,izb) + EndDo + If ( iheat > 0 ) dt9(izb) = rhs(ny+1,izb) EndIf EndDo @@ -433,7 +380,7 @@ Subroutine jacobian_decomp(kstep,mask_in) Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Integer :: i, j, izb, izone, info(zb_lo:zb_hi) + Integer :: i, j, izb, izone Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -494,8 +441,7 @@ Subroutine jacobian_bksub(kstep,yrhs,dy,t9rhs,dt9,mask_in) Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Real(dp) :: rhs(msize) - Integer :: i, izb, izone, info(zb_lo:zb_hi) + Integer :: i, izb, izone Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -509,14 +455,32 @@ Subroutine jacobian_bksub(kstep,yrhs,dy,t9rhs,dt9,mask_in) timer_solve = timer_solve - start_timer timer_bksub = timer_bksub - start_timer + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Do i = 1, ny + rhs(i,izb) = yrhs(i,izb) + EndDo + If ( iheat > 0 ) rhs(ny+1,izb) = t9rhs(izb) + Else + Do i = 1, msize + rhs(i,izb) = 0.0 + EndDo + EndIf + EndDo + ! Solve the LU-decomposed triangular system via back-substitution Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - rhs(1:ny) = yrhs(:,izb) - If ( iheat > 0 ) rhs(ny+1) = t9rhs(izb) - Call dgetrs('N',msize,1,jac(1,1,izb),msize,indx(1,izb),rhs,msize,info(izb)) - dy(:,izb) = rhs(1:ny) - If ( iheat > 0 ) dt9(izb) = rhs(ny+1) + Call dgetrs('N',msize,1,jac(1,1,izb),msize,indx(1,izb),rhs(1,izb),msize,info(izb)) + EndIf + EndDo + + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then + Do i = 1, ny + dy(i,izb) = rhs(i,izb) + EndDo + If ( iheat > 0 ) dt9(izb) = rhs(ny+1,izb) EndIf EndDo diff --git a/source/xnet_nnu.f90 b/source/xnet_nnu.f90 index 6bfb44fd..7f8610d6 100644 --- a/source/xnet_nnu.f90 +++ b/source/xnet_nnu.f90 @@ -26,6 +26,8 @@ Module xnet_nnu Real(dp), Allocatable :: tmevnu(:,:,:) ! Neutrino temperature [MeV] Real(dp), Allocatable :: fluxcms(:,:,:) ! Neutrino fluxes [cm^-2 s^-1] + Real(dp), Allocatable :: ltnu(:,:) ! Interpolated neutrino temperature [MeV] + Real(dp), Allocatable :: fluxnu(:,:) ! Interpolated neutrino fluxes [cm^-2 s^-1] Integer, Allocatable :: irl(:) ! Reaclib Chapter 1 index of reaction Integer, Allocatable :: nuspec(:) ! The neutrino species involved in the reaction @@ -39,6 +41,7 @@ Subroutine read_nnu_data(nnnu,data_dir) !----------------------------------------------------------------------------------------------- ! This routine reads the neutrino cross sections [in units of 10^-42 cm^2] !----------------------------------------------------------------------------------------------- + Use xnet_controls, Only: nzevolve Implicit None ! Input variables @@ -49,7 +52,11 @@ Subroutine read_nnu_data(nnnu,data_dir) Integer :: i, j, lun_nnu Allocate (sigmanu(nnnu,ntnu)) + Allocate (ltnu(nnuspec,nzevolve)) + Allocate (fluxnu(nnuspec,nzevolve)) sigmanu = 0.0 + ltnu = 0.0 + fluxnu = 0.0 Open(newunit=lun_nnu, file=trim(data_dir)//"/netneutr", status='old', action='read') Do i = 1, nnnu @@ -67,45 +74,100 @@ Subroutine nnu_match(nnnu,nr,iwk,innu) ! This routine finds placement in the REACLIB list of each neutrino reaction. ! From this, the neutrino species involved is determined. !----------------------------------------------------------------------------------------------- + Use, Intrinsic :: iso_fortran_env, Only: lun_stdout=>output_unit Use xnet_controls, Only: idiag, lun_diag Implicit None ! Input variables - Integer, Intent(in) :: nnnu, nr, iwk(nr),innu(nr) + Integer, Intent(in) :: nnnu, nr, iwk(nr), innu(nr) ! Local variables Integer :: i, j - Allocate(irl(nnnu),nuspec(nnnu)) + Allocate (irl(nnnu),nuspec(nnnu)) + irl = 0 + nuspec = 0 + Do j = 1, nr - If(innu(j)>0) Then - If (idiag>5) Then - write(lun_diag,*) "[NU]",j,iwk(j) - Endif + If ( innu(j) > 0 ) Then + If ( idiag >= 6 ) Then + Write(lun_diag,"(a,2i5)") "[NU]",j,iwk(j) + EndIf irl(innu(j)) = j - If(iwk(j).eq.7) Then - nuspec(innu(j))=1 - ElseIf(iwk(j).eq.8) Then - nuspec(innu(j))=2 - ElseIf(iwk(j).eq.9 ) Then - nuspec(innu(j))=3 - ElseIf(iwk(j).eq.10) Then - nuspec(innu(j))=4 + If ( iwk(j) == 7 ) Then + nuspec(innu(j)) = 1 + ElseIf (iwk(j) == 8 ) Then + nuspec(innu(j)) = 2 + ElseIf (iwk(j) == 9 ) Then + nuspec(innu(j)) = 3 + ElseIf (iwk(j) == 10 ) Then + nuspec(innu(j)) = 4 Else - Write(6,*) 'nnu match error',j,innu(j),iwk(j) + Write(lun_stdout,"(a,3i5)") 'NNU match error',j,innu(j),iwk(j) EndIf EndIf EndDo - If(idiag==6) Then + If ( idiag>= 6 ) Then Write(lun_diag,"(a)") 'i,IRL,Nuspec' Write(lun_diag,"(3i6)") (i,irl(i),nuspec(i),i=1,nnnu) EndIf Return - End Subroutine nnu_match + Subroutine nnu_flux(tf,nf,ltnuf,fluxf,ts,ns,tnus,fluxs) + Use xnet_types, Only: dp + Use xnet_util, Only: safe_exp + Implicit None + + ! Input variables + Integer, Intent(in) :: ns + Real(dp), Intent(in) :: tf, ts(:), tnus(:,:), fluxs(:,:) + + ! Output variables + Integer, Intent(out) :: nf + Real(dp), Intent(out) :: ltnuf(:), fluxf(:) + + ! Local variables + Real(dp) :: dt, rdt + Integer :: j, n + + ! For constant conditions (ns = 1), do not interpolate in time + If ( ns == 1 ) Then + nf = 1 + + ! Otherwise, intepolate the neutrino fluxes and temperature from time history + Else + Do n = 1, ns + If ( tf <= ts(n) ) Exit + EndDo + nf = n + rdt = ( tf - ts(nf-1)) / ( ts(nf) - ts(nf-1) ) + EndIf + + Do j = 1, nnuspec + If ( nf == 1 ) Then + ltnuf(j) = log(tnus(1,j)) + fluxf(j) = fluxs(1,j) + ElseIf ( nf > 1 .and. nf <= ns ) Then + ltnuf(j) = rdt*log(tnus(nf,j)) + (1.0-rdt)*log(tnus(nf-1,j)) + If ( fluxs(nf,j) == 0.0 .and. fluxs(nf-1,j) == 0.0 ) Then + fluxf(j) = 0.0 + ElseIf ( fluxs(nf,j) == 0.0 .or. fluxs(nf-1,j) == 0.0 ) Then + fluxf(j) = rdt*fluxs(nf,j) + (1.0-rdt)*fluxs(nf-1,j) + Else + fluxf(j) = safe_exp( rdt*log(fluxs(nf,j)) + (1.0-rdt)*log(fluxs(nf-1,j)) ) + EndIf + Else + ltnuf(j) = log(tnus(ns,j)) + fluxf(j) = fluxs(ns,j) + EndIf + EndDo + + Return + End Subroutine nnu_flux + Subroutine nnu_rate(nnnu,time,rate,mask_in) !----------------------------------------------------------------------------------------------- ! This routine calculates the neutrino nucleus reaction rates [in s^-1] as @@ -116,14 +178,14 @@ Subroutine nnu_rate(nnnu,time,rate,mask_in) ! from neutrino luminosities. !----------------------------------------------------------------------------------------------- Use xnet_conditions, Only: nh, th - Use xnet_controls, Only: nzevolve, zb_lo, zb_hi, lzactive, ineutrino, idiag, lun_diag, szbatch + Use xnet_controls, Only: zb_lo, zb_hi, lzactive, ineutrino, idiag, lun_diag, szbatch Use xnet_types, Only: dp Use xnet_util, Only: safe_exp Implicit None ! Input variables Integer, Intent(in) :: nnnu - Real(dp), Intent(in) :: time(nzevolve) + Real(dp), Intent(in) :: time(zb_lo:zb_hi) ! Output variables Real(dp), Intent(out) :: rate(nnnu,nnuspec,zb_lo:zb_hi) @@ -132,10 +194,9 @@ Subroutine nnu_rate(nnnu,time,rate,mask_in) Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables - Real(dp) :: ltnu, flux(nnuspec) Real(dp) :: lsigmanu1, lsigmanu2 - Real(dp) :: rcsnu(nnnu) - Real(dp) :: rdt, rdltnu, ltnu1, ltnu2, lfl1,lfl2 + Real(dp) :: rcsnu + Real(dp) :: rdt, rdltnu, ltnu1, ltnu2, lfl1, lfl2 Real(dp) :: xrate Integer :: izb, it, i, j, k, n, izone Logical, Pointer :: mask(:) @@ -159,45 +220,22 @@ Subroutine nnu_rate(nnnu,time,rate,mask_in) EndIf EndDo Else + + ! Interpolate flux and neutrino temperature from time history Do izb = zb_lo, zb_hi If ( mask(izb) ) Then + Call nnu_flux(time(izb),n,ltnu(:,izb),fluxnu(:,izb), & + th(:,izb),nh(izb),tmevnu(:,:,izb),fluxcms(:,:,izb)) + EndIf + EndDo - ! For constant conditions (nh = 1), do not interpolate in time - If ( nh(izb) == 1 ) Then - n = 1 - - ! Otherwise, intepolate the neutrino fluxes and temperature from time history - Else - Do i = 1, nh(izb) - If ( time(izb) <= th(i,izb) ) Exit - EndDo - n = i - EndIf - - ! Compute neutrino cross sections - Do j = 1, nnuspec + ! Compute neutrino cross sections + Do izb = zb_lo, zb_hi + Do j = 1, nnuspec + If ( mask(izb) ) Then - ! Linear interpolation in log-space - If ( n == 1 ) Then - ltnu = log(tmevnu(1,j,izb)) - flux(j) = fluxcms(1,j,izb) - ElseIf ( n > nh(izb) ) Then - ltnu = log(tmevnu(nh(izb),j,izb)) - flux(j) = fluxcms(nh(izb),j,izb) - Else - rdt = (time(izb)-th(n-1,izb)) / (th(n,izb)-th(n-1,izb)) - ltnu = rdt*log(tmevnu(n,j,izb)) + (1.0-rdt)*log(tmevnu(n-1,j,izb)) - - If ( fluxcms(n,j,izb) ==0. .and. fluxcms(n-1,j,izb)==0. ) then - flux(j) = 0. - Else If ( fluxcms(n,j,izb) ==0. .or. fluxcms(n-1,j,izb)==0. ) then - flux(j) = rdt*fluxcms(n,j,izb) + (1.0-rdt)*fluxcms(n-1,j,izb) - Else - flux (j) = safe_exp( rdt*log(fluxcms(n,j,izb)) + (1.0-rdt)*log(fluxcms(n-1,j,izb)) ) - Endif - EndIf Do i = 1, ntnu - If ( ltnu <= ltnugrid(i) ) Exit + If ( ltnu(j,izb) <= ltnugrid(i) ) Exit EndDo it = i @@ -205,35 +243,34 @@ Subroutine nnu_rate(nnnu,time,rate,mask_in) ! Log interpolation If ( it == 1 ) Then - rcsnu(k) = sigmanu(k,1) + rcsnu = sigmanu(k,1) ElseIf ( it > ntnu ) Then - rcsnu(k) = sigmanu(k,ntnu) + rcsnu = sigmanu(k,ntnu) Else ltnu1 = ltnugrid(it-1) ltnu2 = ltnugrid(it) lsigmanu1 = log(sigmanu(k,it-1)) lsigmanu2 = log(sigmanu(k,it)) - rdltnu = (ltnu-ltnu1) / (ltnu2-ltnu1) - rcsnu(k) = safe_exp( rdltnu*lsigmanu2 + (1.0-rdltnu)*lsigmanu1 ) + rdltnu = (ltnu(j,izb)-ltnu1) / (ltnu2-ltnu1) + rcsnu = safe_exp( rdltnu*lsigmanu2 + (1.0-rdltnu)*lsigmanu1 ) EndIf - ! Calculate rate only for neutrino specie involved in reaction - If (nuspec(k) == j .or. (nuspec(k)==3 .and. j==1) .or. (nuspec(k)==4 .and. j==2) ) Then - xrate = flux(j)*rcsnu(k)*1e-42 - If (xrate > 1.d-80) Then - rate(k,j,izb) = xrate + ! Calculate rate only for neutrino species involved in reaction + If ( nuspec(k) == j .or. & + & ( nuspec(k) == 3 .and. j == 1 ) .or. & + & ( nuspec(k) == 4 .and. j == 2 ) ) Then + xrate = fluxnu(j,izb)*rcsnu*1e-42 ! Not a fan of 1e-42 being hard-coded here + If ( xrate > 1.0e-80 ) Then + rate(k,j,izb) = xrate Else - rate(k,j,izb) = 0.0d0 + rate(k,j,izb) = 0.0 Endif - -! write(*,*) "[NU rate]",k,j,nuspec(k),flux(j),rcsnu(k) -! write(lun_diag,*) "[NU rate]",k,j,flux(j),rcsnu(k) Else rate(k,j,izb) = 0.0 EndIf EndDo - EndDo - EndIf + EndIf + EndDo EndDo EndIf @@ -243,8 +280,8 @@ Subroutine nnu_rate(nnnu,time,rate,mask_in) izone = izb + szbatch - zb_lo Write(lun_diag,"(a,i5)") 'NNU',izone Write(lun_diag,"(a,2i5)") 'Nuspec,Nnnu',nnuspec,nnnu - Write(lun_diag,"(a,4es23.15)") 'Nu Flux',flux - Write(lun_diag,"(i5,5es13.5)") (k, rcsnu(k), rate(k,:,izb), k=1,nnnu) + Write(lun_diag,"(a,4es23.15)") 'Nu Flux',(fluxnu(j,izb), j=1,nnuspec) + Write(lun_diag,"(i5,5es13.5)") (k, rcsnu, (rate(k,j,izb),j=1,nnuspec), k=1,nnnu) EndIf EndDo EndIf diff --git a/source/xnet_screening.f90 b/source/xnet_screening.f90 index 89678438..d5a42717 100644 --- a/source/xnet_screening.f90 +++ b/source/xnet_screening.f90 @@ -26,12 +26,17 @@ Module xnet_screening ! Screening factors from Table 4 of Graboske+ (1973) Real(dp), Allocatable :: zeta2w(:), zeta2i(:), zeta3w(:), zeta3i(:), zeta4w(:), zeta4i(:) ! Reaction charge parameter + Real(dp), Allocatable :: ztilde(:), zinter(:), lambda0(:), gammae(:), dztildedt9(:) + Real(dp), Parameter :: lam_1 = 0.1_dp Real(dp), Parameter :: lam_2 = 0.125_dp Real(dp), Parameter :: lam_3 = 2.0_dp Real(dp), Parameter :: lam_4 = 2.15_dp Real(dp), Parameter :: lam_5 = 4.85_dp Real(dp), Parameter :: lam_6 = 5.0_dp + Real(dp), Parameter :: rdlam12 = 1.0_dp / (lam_1 - lam_2) + Real(dp), Parameter :: rdlam34 = 1.0_dp / (lam_3 - lam_4) + Real(dp), Parameter :: rdlam56 = 1.0_dp / (lam_6 - lam_5) Contains @@ -100,26 +105,35 @@ Subroutine screening_init z4c = real(iz4c,dp) zeta4w = z41*(z42 + z43 + z44) + z42*(z43 + z44) + z43 * z44 zeta4i = z4c**bip1 - z41**bip1 - z42**bip1 - z43**bip1 - z44**bip1 + + Allocate (ztilde(nzevolve)) + Allocate (zinter(nzevolve)) + Allocate (lambda0(nzevolve)) + Allocate (gammae(nzevolve)) + Allocate (dztildedt9(nzevolve)) + ztilde = 0.0 + zinter = 0.0 + lambda0 = 0.0 + gammae = 0.0 + dztildedt9 = 0.0 EndIf Allocate (h1(nr1,nzevolve)) Allocate (h2(nr2,nzevolve)) Allocate (h3(nr3,nzevolve)) Allocate (h4(nr4,nzevolve)) - h1 = 0.0d0 - h2 = 0.0d0 - h3 = 0.0d0 - h4 = 0.0d0 - If ( iheat > 0 ) Then - Allocate (dh1dt9(nr1,nzevolve)) - Allocate (dh2dt9(nr2,nzevolve)) - Allocate (dh3dt9(nr3,nzevolve)) - Allocate (dh4dt9(nr4,nzevolve)) - dh1dt9 = 0.0d0 - dh2dt9 = 0.0d0 - dh3dt9 = 0.0d0 - dh4dt9 = 0.0d0 - EndIf + Allocate (dh1dt9(nr1,nzevolve)) + Allocate (dh2dt9(nr2,nzevolve)) + Allocate (dh3dt9(nr3,nzevolve)) + Allocate (dh4dt9(nr4,nzevolve)) + h1 = 0.0 + h2 = 0.0 + h3 = 0.0 + h4 = 0.0 + dh1dt9 = 0.0 + dh2dt9 = 0.0 + dh3dt9 = 0.0 + dh4dt9 = 0.0 Return End Subroutine screening_init @@ -154,18 +168,12 @@ Subroutine screening(mask_in) Logical, Optional, Target, Intent(in) :: mask_in(zb_lo:zb_hi) ! Local variables + Real(dp) :: fhs(0:izmax+2), dfhsdt9(0:izmax+2) + Real(dp) :: gammaz, gammaz5, lgammaz, lambda + Real(dp) :: hw0, hi0, dlnhw0dt9, dlnhi0dt9 + Real(dp) :: h, hw, hi, hs + Real(dp) :: dhdt9, dhwdt9, dhidt9, dhsdt9 Integer :: j, mu, izb, izone - Real(dp), Dimension(nreac(2)) :: h2w, h2i, h2s, lambda12 - Real(dp), Dimension(nreac(2)) :: dh2wdt9, dh2idt9, dh2sdt9 - Real(dp), Dimension(nreac(2)) :: theta12iw, theta12is, theta12si - Real(dp), Dimension(nreac(3)) :: h3w, h3i, h3s, lambda123 - Real(dp), Dimension(nreac(3)) :: dh3wdt9, dh3idt9, dh3sdt9 - Real(dp), Dimension(nreac(3)) :: theta123iw, theta123is, theta123si - Real(dp), Dimension(nreac(4)) :: h4w, h4i, h4s, lambda1234 - Real(dp), Dimension(nreac(4)) :: dh4wdt9, dh4idt9, dh4sdt9 - Real(dp), Dimension(nreac(4)) :: theta1234iw, theta1234is, theta1234si - Real(dp), Dimension(0:izmax+2) :: gammaz, fhs, fhi, dfhsdt9 - Real(dp) :: ztilde, zinter, lambda0, gammae, dztildedt9 Logical, Pointer :: mask(:) If ( present(mask_in) ) Then @@ -175,199 +183,197 @@ Subroutine screening(mask_in) EndIf If ( .not. any(mask) ) Return + start_timer = xnet_wtime() + timer_prescrn = timer_prescrn - start_timer + + ! Call EOS to get plasma quantities + Call eos_screen(t9t(zb_lo:zb_hi),rhot(zb_lo:zb_hi),yt(:,zb_lo:zb_hi),etae(zb_lo:zb_hi), & + & detaedt9(zb_lo:zb_hi),ztilde(zb_lo:zb_hi),zinter(zb_lo:zb_hi),lambda0(zb_lo:zb_hi), & + & gammae(zb_lo:zb_hi),dztildedt9(zb_lo:zb_hi),xext(zb_lo:zb_hi),aext(zb_lo:zb_hi), & + & zext(zb_lo:zb_hi),mask_in = mask) + + stop_timer = xnet_wtime() + timer_prescrn = timer_prescrn + stop_timer + start_timer = xnet_wtime() timer_scrn = timer_scrn - start_timer Do izb = zb_lo, zb_hi If ( mask(izb) ) Then - ! Call EOS to get plasma quantities - call eos_screen(t9t(izb),rhot(izb),yt(:,izb),etae(izb),detaedt9(izb), & - & ztilde,zinter,lambda0,gammae,dztildedt9,xext(izb),aext(izb),zext(izb)) - ! Calculate screening energies as a function of Z, for prescriptions that follow this approach - gammaz(0) = 0.0d0 - gammaz(1:izmax+2) = gammae * zseq53(1:izmax+2) - fhi(0) = 0.0d0 - fhi(1:izmax+2) = kbi * zinter * lambda0**bi * zseqi(1:izmax+2) - fhs(0) = 0.0d0 - fhs(1:izmax+2) = + cds(1) * gammaz(1:izmax+2) & - & + cds(2) * gammaz(1:izmax+2)**cds(5) / cds(5) & - & + cds(3) * log(gammaz(1:izmax+2)) & - & + cds(4) - dfhsdt9(0) = 0.0d0 - dfhsdt9(1:izmax+2) = + cds(1) * gammaz(1:izmax+2) & - & + cds(2) * gammaz(1:izmax+2)**cds(5) & - & + cds(3) - dfhsdt9(1:izmax+2) = -dfhsdt9(1:izmax+2) / t9t(izb) + fhs(0) = 0.0 + dfhsdt9(0) = 0.0 + Do j = 1, izmax+2 + gammaz = gammae(izb) * zseq53(j) + gammaz5 = gammaz**cds(5) + lgammaz = log(gammaz) + fhs(j) = cds(1)*gammaz + cds(2)/cds(5)*gammaz5 + cds(3)*lgammaz + cds(4) + dfhsdt9(j) = -(cds(1)*gammaz + cds(2)*gammaz5 + cds(3))/t9t(izb) + EndDo ! No screening term for 1-reactant reactions - h1(:,izb) = 0.0d0 - If ( iheat > 0 ) dh1dt9(:,izb) = 0.0d0 - - ! Weak and intermediate screening factors, Table 4 of Graboske et al. (1973) - lambda12 = zeta2w * ztilde * lambda0 - h2w = lambda12 - h2i = kbi * zinter * lambda0**bi * zeta2i - !h2i = fhi(iz2c) - fhi(iz21) - fhi(iz22) - - ! Strong screening from Dewitt & Slattery (2003) using linear mixing. - h2s = fhs(iz21) + fhs(iz22) - fhs(iz2c) - - ! Blending factors - theta12iw = max( 0.0d0, min( 1.0d0, (lambda12 - lam_1) / (lam_2 - lam_1) ) ) - theta12is = max( 0.0d0, min( 1.0d0, (lambda12 - lam_5) / (lam_6 - lam_5) ) ) - theta12si = max( 0.0d0, min( 1.0d0, (lambda12 - lam_3) / (lam_4 - lam_3) ) ) - - ! Select Screening factor for 2 reactant reactions - Where ( iz2c == 0 ) - h2(:,izb) = 0.0d0 - ElseWhere ( lambda12 < lam_1 ) - h2(:,izb) = h2w - ElseWhere ( lambda12 < lam_3 ) - h2(:,izb) = theta12iw*h2i + (1.0d0 - theta12iw)*h2w - ElseWhere ( lambda12 > lam_6 ) - h2(:,izb) = h2s - ElseWhere ( h2i < h2s ) - h2(:,izb) = theta12is*h2s + (1.0d0 - theta12is)*h2i - ElseWhere - h2(:,izb) = theta12si*h2s + (1.0d0 - theta12si)*h2i - EndWhere - If ( iheat > 0 ) Then - dh2wdt9 = +h2w * (dztildedt9/ztilde - 1.5/t9t(izb)) - dh2idt9 = -h2i * (thbim2*dztildedt9/ztilde + bi*1.5/t9t(izb)) - dh2sdt9 = dfhsdt9(iz21) + dfhsdt9(iz22) - dfhsdt9(iz2c) - Where ( iz2c == 0 ) - dh2dt9(:,izb) = 0.0d0 - ElseWhere ( lambda12 < lam_1 ) - dh2dt9(:,izb) = dh2wdt9 - ElseWhere ( lambda12 < lam_3 ) - dh2dt9(:,izb) = theta12iw*dh2idt9 + (1.0d0 - theta12iw)*dh2wdt9 - ElseWhere ( lambda12 > lam_6 ) - dh2dt9(:,izb) = dh2sdt9 - ElseWhere ( h2i < h2s ) - dh2dt9(:,izb) = theta12is*dh2sdt9 + (1.0d0 - theta12is)*dh2idt9 - ElseWhere - dh2dt9(:,izb) = theta12si*dh2sdt9 + (1.0d0 - theta12si)*dh2idt9 - EndWhere - EndIf - - ! Weak and intermediate screening factors, Table 4 of Graboske+ (1973) - lambda123 = zeta3w * ztilde * lambda0 - h3w = lambda123 - !h3i = kbi * zinter * lambda0**bi * zeta3i - h3i = fhi(iz3c) - fhi(iz31) - fhi(iz32) - fhi(iz33) - - ! Strong screening from Dewitt & Slattery (2003) using linear mixing. - h3s = fhs(iz31) + fhs(iz32) + fhs(iz33) - fhs(iz3c) - - ! Blending factors - theta123iw = max( 0.0d0, min( 1.0d0, (lambda123 - lam_1) / (lam_2 - lam_1) ) ) - theta123is = max( 0.0d0, min( 1.0d0, (lambda123 - lam_5) / (lam_6 - lam_5) ) ) - theta123si = max( 0.0d0, min( 1.0d0, (lambda123 - lam_3) / (lam_4 - lam_3) ) ) - - ! Select screening factor for 3 reactant reactions - Where ( iz3c == 0 ) - h3(:,izb) = 0.0d0 - ElseWhere ( lambda123 < lam_1 ) - h3(:,izb) = h3w - ElseWhere ( lambda123 < lam_3 ) - h3(:,izb) = theta123iw*h3i + (1.0d0 - theta123iw)*h3w - ElseWhere ( lambda123 > lam_6 ) - h3(:,izb) = h3s - ElseWhere ( h3i < h3s ) - h3(:,izb) = theta123is*h3s + (1.0d0 - theta123is)*h3i - ElseWhere - h3(:,izb) = theta123si*h3s + (1.0d0 - theta123si)*h3i - EndWhere - - If ( iheat > 0 ) Then - dh3wdt9 = +h3w * (dztildedt9/ztilde - 1.5d0/t9t(izb)) - dh3idt9 = -h3i * (thbim2*dztildedt9/ztilde + bi*1.5d0/t9t(izb)) - dh3sdt9 = dfhsdt9(iz31) + dfhsdt9(iz32) + dfhsdt9(iz33) - dfhsdt9(iz3c) - Where ( iz3c == 0 ) - dh3dt9(:,izb) = 0.0d0 - ElseWhere ( lambda123 < lam_1 ) - dh3dt9(:,izb) = dh3wdt9 - ElseWhere ( lambda123 < lam_3 ) - dh3dt9(:,izb) = theta123iw*dh3idt9 + (1.0d0 - theta123iw)*dh3wdt9 - ElseWhere ( lambda123 > lam_6 ) - dh3dt9(:,izb) = dh3sdt9 - ElseWhere ( h3i < h3s ) - dh3dt9(:,izb) = theta123is*dh3sdt9 + (1.0d0 - theta123is)*dh3idt9 - ElseWhere - dh3dt9(:,izb) = theta123si*dh3sdt9 + (1.0d0 - theta123si)*dh3idt9 - EndWhere - EndIf - - ! Weak and intermediate screening factors, Table 4 of Graboske+ (1973) - lambda1234 = zeta4w * ztilde * lambda0 - h4w = lambda1234 - !h4i = kbi * zinter * lambda0**bi * zeta4i - h4i = fhi(iz4c) - fhi(iz41) - fhi(iz42) - fhi(iz43) - fhi(iz44) - - ! Strong screening from Dewitt & Slattery (2003) using linear mixing. - h4s = fhs(iz41) + fhs(iz42) + fhs(iz43) + fhs(iz44) - fhs(iz4c) - - ! Blending factors - theta1234iw = max( 0.0d0, min( 1.0d0, (lambda1234 - lam_1) / (lam_2 - lam_1) ) ) - theta1234is = max( 0.0d0, min( 1.0d0, (lambda1234 - lam_5) / (lam_6 - lam_5) ) ) - theta1234si = max( 0.0d0, min( 1.0d0, (lambda1234 - lam_3) / (lam_4 - lam_3) ) ) - - ! Select screening factor for 4 reactant reactions - Where ( iz4c == 0 ) - h4(:,izb) = 0.0d0 - ElseWhere ( lambda1234 < lam_1 ) - h4(:,izb) = h4w - ElseWhere ( lambda1234 < lam_3 ) - h4(:,izb) = theta1234iw*h4i + (1.0d0 - theta1234iw)*h4w - ElseWhere ( lambda1234 > lam_6 ) - h4(:,izb) = h4s - ElseWhere ( h4i < h4s ) - h4(:,izb) = theta1234is*h4s + (1.0d0 - theta1234is)*h4i - ElseWhere - h4(:,izb) = theta1234si*h4s + (1.0d0 - theta1234si)*h4i - EndWhere - - If ( iheat > 0 ) Then - dh4wdt9 = +h4w * (dztildedt9/ztilde - 1.5d0/t9t(izb)) - dh4idt9 = -h4i * (thbim2*dztildedt9/ztilde + bi*1.5d0/t9t(izb)) - dh4sdt9 = dfhsdt9(iz41) + dfhsdt9(iz42) + dfhsdt9(iz43) + dfhsdt9(iz44) - dfhsdt9(iz4c) - Where ( iz4c == 0 ) - dh4dt9(:,izb) = 0.0d0 - ElseWhere ( lambda1234 < lam_1 ) - dh4dt9(:,izb) = dh4wdt9 - ElseWhere ( lambda1234 < lam_3 ) - dh4dt9(:,izb) = theta1234iw*dh4idt9 + (1.0d0 - theta1234iw)*dh4wdt9 - ElseWhere ( lambda1234 > lam_6 ) - dh4dt9(:,izb) = dh4sdt9 - ElseWhere ( h4i < h4s ) - dh4dt9(:,izb) = theta1234is*dh4sdt9 + (1.0d0 - theta1234is)*dh4idt9 - ElseWhere - dh4dt9(:,izb) = theta1234si*dh4sdt9 + (1.0d0 - theta1234si)*dh4idt9 - EndWhere - EndIf + Do mu = 1, nreac(1) + h1(mu,izb) = 0.0 + dh1dt9(mu,izb) = 0.0 + EndDo + + hw0 = ztilde(izb) * lambda0(izb) + hi0 = kbi * zinter(izb) * lambda0(izb)**bi + dlnhw0dt9 = dztildedt9(izb)/ztilde(izb) - 1.5/t9t(izb) + dlnhi0dt9 = - thbim2*dztildedt9(izb)/ztilde(izb) - 1.5*bi/t9t(izb) + + ! 2-reactant screening + Do mu = 1, nreac(2) + If ( zeta2w(mu) > 0.0 ) Then + + lambda = hw0 * zeta2w(mu) + hw = lambda + hi = hi0 * zeta2i(mu) + hs = fhs(iz21(mu)) + fhs(iz22(mu)) - fhs(iz2c(mu)) + + dhwdt9 = hw * dlnhw0dt9 + dhidt9 = hi * dlnhi0dt9 + dhsdt9 = dfhsdt9(iz21(mu)) + dfhsdt9(iz22(mu)) - dfhsdt9(iz2c(mu)) + + ! Select Screening factor for 2 reactant reactions + Call screen_blend(lambda,hw,hi,hs,dhwdt9,dhidt9,dhsdt9,h,dhdt9) + h2(mu,izb) = h + dh2dt9(mu,izb) = dhdt9 + Else + h2(mu,izb) = 0.0 + dh2dt9(mu,izb) = 0.0 + EndIf + EndDo + + ! 3-reactant screening + Do mu = 1, nreac(3) + If ( zeta3w(mu) > 0.0 ) Then + + lambda = hw0 * zeta3w(mu) + hw = lambda + hi = hi0 * zeta3i(mu) + hs = fhs(iz31(mu)) + fhs(iz33(mu)) - fhs(iz3c(mu)) + + dhwdt9 = hw * dlnhw0dt9 + dhidt9 = hi * dlnhi0dt9 + dhsdt9 = dfhsdt9(iz31(mu)) + dfhsdt9(iz32(mu)) + dfhsdt9(iz33(mu)) - dfhsdt9(iz3c(mu)) + + ! Select Screening factor for 3 reactant reactions + Call screen_blend(lambda,hw,hi,hs,dhwdt9,dhidt9,dhsdt9,h,dhdt9) + h3(mu,izb) = h + dh3dt9(mu,izb) = dhdt9 + Else + h3(mu,izb) = 0.0 + dh3dt9(mu,izb) = 0.0 + EndIf + EndDo + + ! 4-reactant screening + Do mu = 1, nreac(4) + If ( zeta4w(mu) > 0.0 ) Then + + lambda = hw0 * zeta4w(mu) + hw = lambda + hi = hi0 * zeta4i(mu) + hs = fhs(iz41(mu)) + fhs(iz42(mu)) + fhs(iz43(mu)) + fhs(iz44(mu)) - fhs(iz4c(mu)) + + dhwdt9 = hw * dlnhw0dt9 + dhidt9 = hi * dlnhi0dt9 + dhsdt9 = dfhsdt9(iz41(mu)) + dfhsdt9(iz42(mu)) + dfhsdt9(iz43(mu)) + dfhsdt9(iz44(mu)) - dfhsdt9(iz4c(mu)) + + ! Select Screening factor for 4 reactant reactions + Call screen_blend(lambda,hw,hi,hs,dhwdt9,dhidt9,dhsdt9,h,dhdt9) + h4(mu,izb) = h + dh4dt9(mu,izb) = dhdt9 + Else + h4(mu,izb) = 0.0 + dh4dt9(mu,izb) = 0.0 + EndIf + EndDo + EndIf + EndDo - If ( idiag >= 5 ) Then + If ( idiag >= 5 ) Then + Do izb = zb_lo, zb_hi + If ( mask(izb) ) Then izone = izb + szbatch - zb_lo Write(lun_diag,"(a,i5)") 'SCREEN',izone - Write(lun_diag,"(3a5,i6,5es23.15)") & - & ('H2',(nname(n2i(j,mu)),j=1,2),mu,lambda12(mu), h2(mu,izb),h2w(mu),h2i(mu),h2s(mu),mu=1,nreac(2)) - Write(lun_diag,"(4a5,i6,5es23.15)") & - & ('H3',(nname(n3i(j,mu)),j=1,3),mu,lambda123(mu), h3(mu,izb),h3w(mu),h3i(mu),h3s(mu),mu=1,nreac(3)) - Write(lun_diag,"(5a5,i6,5es23.15)") & - & ('H4',(nname(n4i(j,mu)),j=1,4),mu,lambda1234(mu),h4(mu,izb),h4w(mu),h4i(mu),h4s(mu),mu=1,nreac(4)) + hw0 = ztilde(izb) * lambda0(izb) + hi0 = kbi * zinter(izb) * lambda0(izb)**bi + dlnhw0dt9 = dztildedt9(izb)/ztilde(izb) - 1.5/t9t(izb) + dlnhi0dt9 = - thbim2*dztildedt9(izb)/ztilde(izb) - 1.5*bi/t9t(izb) + fhs(0) = 0.0 + dfhsdt9(0) = 0.0 + Do j = 1, izmax+2 + gammaz = gammae(izb) * zseq53(j) + gammaz5 = gammaz**cds(5) + lgammaz = log(gammaz) + fhs(j) = cds(1)*gammaz + cds(2)/cds(5)*gammaz5 + cds(3)*lgammaz + cds(4) + dfhsdt9(j) = -(cds(1)*gammaz + cds(2)*gammaz5 + cds(3))/t9t(izb) + EndDo + Do mu = 1, nreac(2) + lambda = hw0 * zeta2w(mu) + hw = lambda + hi = hi0 * zeta2i(mu) + hs = fhs(iz21(mu)) + fhs(iz22(mu)) - fhs(iz2c(mu)) + Write(lun_diag,"(3a5,i6,5es23.15)") & + & 'H2',(nname(n2i(j,mu)),j=1,2),mu,lambda,h2(mu,izb),hw,hi,hs + EndDo + Do mu = 1, nreac(3) + lambda = hw0 * zeta3w(mu) + hw = lambda + hi = hi0 * zeta3i(mu) + hs = fhs(iz31(mu)) + fhs(iz32(mu)) + fhs(iz33(mu)) - fhs(iz3c(mu)) + Write(lun_diag,"(4a5,i6,5es23.15)") & + & 'H3',(nname(n3i(j,mu)),j=1,3),mu,lambda,h3(mu,izb),hw,hi,hs + EndDo + Do mu = 1, nreac(4) + lambda = hw0 * zeta4w(mu) + hw = lambda + hi = hi0 * zeta4i(mu) + hs = fhs(iz41(mu)) + fhs(iz42(mu)) + fhs(iz43(mu)) + fhs(iz44(mu)) - fhs(iz4c(mu)) + Write(lun_diag,"(5a5,i6,5es23.15)") & + & 'H4',(nname(n4i(j,mu)),j=1,4),mu,lambda,h4(mu,izb),hw,hi,hs + EndDo If ( iheat > 0 ) Then - Write(lun_diag,"(a7,2a5,i6,4es23.15)") & - & ('dH2/dT9',(nname(n2i(j,mu)),j=1,2),mu,dh2dt9(mu,izb),dh2wdt9(mu),dh2idt9(mu),dh2sdt9(mu),mu=1,nreac(2)) - Write(lun_diag,"(a7,3a5,i6,4es23.15)") & - & ('dH3/dT9',(nname(n3i(j,mu)),j=1,3),mu,dh3dt9(mu,izb),dh3wdt9(mu),dh3idt9(mu),dh3sdt9(mu),mu=1,nreac(3)) - Write(lun_diag,"(a7,4a5,i6,4es23.15)") & - & ('dH4/dT9',(nname(n4i(j,mu)),j=1,4),mu,dh4dt9(mu,izb),dh4wdt9(mu),dh4idt9(mu),dh4sdt9(mu),mu=1,nreac(4)) + Do mu = 1, nreac(2) + lambda = hw0 * zeta2w(mu) + hw = lambda + hi = hi0 * zeta2i(mu) + hs = fhs(iz21(mu)) + fhs(iz22(mu)) - fhs(iz2c(mu)) + dhwdt9 = hw * dlnhw0dt9 + dhidt9 = hi * dlnhi0dt9 + dhsdt9 = dfhsdt9(iz21(mu)) + dfhsdt9(iz22(mu)) - dfhsdt9(iz2c(mu)) + Write(lun_diag,"(a7,2a5,i6,4es23.15)") & + & 'dH2/dT9',(nname(n2i(j,mu)),j=1,2),mu,dh2dt9(mu,izb),dhwdt9,dhidt9,dhsdt9 + EndDo + Do mu = 1, nreac(3) + lambda = hw0 * zeta3w(mu) + hw = lambda + hi = hi0 * zeta3i(mu) + hs = fhs(iz31(mu)) + fhs(iz32(mu)) + fhs(iz33(mu)) - fhs(iz3c(mu)) + dhwdt9 = hw * dlnhw0dt9 + dhidt9 = hi * dlnhi0dt9 + dhsdt9 = dfhsdt9(iz31(mu)) + dfhsdt9(iz32(mu)) + dfhsdt9(iz33(mu)) - dfhsdt9(iz3c(mu)) + Write(lun_diag,"(a7,3a5,i6,4es23.15)") & + & 'dH3/dT9',(nname(n3i(j,mu)),j=1,3),mu,dh3dt9(mu,izb),dhwdt9,dhidt9,dhsdt9 + EndDo + Do mu = 1, nreac(4) + lambda = hw0 * zeta4w(mu) + hw = lambda + hi = hi0 * zeta4i(mu) + hs = fhs(iz41(mu)) + fhs(iz42(mu)) + fhs(iz43(mu)) + fhs(iz44(mu)) - fhs(iz4c(mu)) + dhwdt9 = hw * dlnhw0dt9 + dhidt9 = hi * dlnhi0dt9 + dhsdt9 = dfhsdt9(iz41(mu)) + dfhsdt9(iz42(mu)) + dfhsdt9(iz43(mu)) + dfhsdt9(iz44(mu)) - dfhsdt9(iz4c(mu)) + Write(lun_diag,"(a7,4a5,i6,4es23.15)") & + & 'dH4/dT9',(nname(n4i(j,mu)),j=1,4),mu,dh4dt9(mu,izb),dhwdt9,dhidt9,dhsdt9 + EndDo EndIf EndIf - EndIf - EndDo + EndDo + EndIf stop_timer = xnet_wtime() timer_scrn = timer_scrn + stop_timer @@ -375,4 +381,56 @@ Subroutine screening(mask_in) Return End Subroutine screening + Subroutine screen_blend(lambda,hw,hi,hs,dhwdt9,dhidt9,dhsdt9,h,dhdt9) + !----------------------------------------------------------------------------------------------- + ! This function linearly blends screening prescriptions + !----------------------------------------------------------------------------------------------- + Use xnet_types, Only: dp + Implicit None + + ! Input variables + Real(dp), Intent(in) :: lambda, hw, hi, hs, dhwdt9, dhidt9, dhsdt9 + + ! Output variables + Real(dp), Intent(out) :: h, dhdt9 + + ! Local variables + Real(dp) :: theta + + If ( lambda <= lam_1 ) Then + h = hw + dhdt9 = dhwdt9 + ElseIf ( lambda <= lam_2 ) Then + theta = (lambda - lam_2) * rdlam12 + h = theta*hw + (1.0 - theta)*hi + dhdt9 = theta*dhwdt9 + (1.0 - theta)*dhidt9 + ElseIf ( lambda <= lam_3 ) Then + h = hi + dhdt9 = dhidt9 + ElseIf ( lambda > lam_6 ) Then + h = hs + dhdt9 = dhsdt9 + ElseIf ( hi >= hs ) Then + If ( lambda <= lam_4 ) Then + theta = (lambda - lam_4) * rdlam34 + h = theta*hi + (1.0 - theta)*hs + dhdt9 = theta*dhidt9 + (1.0 - theta)*dhsdt9 + Else + h = hs + dhdt9 = dhsdt9 + EndIf + Else + If ( lambda <= lam_5 ) Then + h = hi + dhdt9 = dhidt9 + Else + theta = (lambda - lam_5) * rdlam56 + h = theta*hs + (1.0 - theta)*hi + dhdt9 = theta*dhsdt9 + (1.0 - theta)*dhidt9 + EndIf + EndIf + + Return + End Subroutine screen_blend + End Module xnet_screening