Fix some incorrect derivative terms in sneut5 (#1633)

AMReX-Astro · Jul 29, 2024 · 83ce77c · 83ce77c
1 parent b42f181
commit 83ce77c
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 6 deletions.
diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H
@@ -686,20 +686,23 @@ void nu_photo(const sneutf_t& sf,
     // and table 2 written out for speed
 
     amrex::Real tau;
-    amrex::Real cc;
+    amrex::Real cc, ccdt;
     amrex::Real c00, c01, c02, c03, c04, c05, c06;
     amrex::Real c10, c11, c12, c13, c14, c15, c16;
     amrex::Real c20, c21, c22, c23, c24, c25, c26;
     amrex::Real dd01, dd02, dd03, dd04, dd05;
     amrex::Real dd11, dd12, dd13, dd14, dd15;
     amrex::Real dd21, dd22, dd23, dd24, dd25;
 
+    amrex::Real taudt = nu_constants::iln10 * sf.tempi;
+
     if (sf.temp < 1.0e8_rt) {
 
         // note: we already bailed above for T < 1.e7, so this is really 1.e7 <= T < 1.e8
 
         tau  =  std::log10(sf.temp * 1.0e-7_rt);
         cc   =  0.5654e0_rt + tau;
+        ccdt =  taudt;
         c00  =  1.008e11_rt;
         c01  =  0.0e0_rt;
         c02  =  0.0e0_rt;
@@ -741,6 +744,7 @@ void nu_photo(const sneutf_t& sf,
 
         tau  =  std::log10(sf.temp * 1.0e-8_rt);
         cc   =  1.5654e0_rt;
+        ccdt =  0.0e0_rt;
         c00  =  9.889e10_rt;
         c01  = -4.524e8_rt;
         c02  = -6.088e6_rt;
@@ -784,6 +788,7 @@ void nu_photo(const sneutf_t& sf,
 
         tau  =  std::log10(sf.t9);
         cc   =  1.5654e0_rt;
+        ccdt =  0.0e0_rt;
         c00  =  9.581e10_rt;
         c01  =  4.107e8_rt;
         c02  =  2.305e8_rt;
@@ -823,8 +828,6 @@ void nu_photo(const sneutf_t& sf,
 
     }
 
-    amrex::Real taudt = nu_constants::iln10 * sf.tempi;
-
     // equation 3.7
 
     const auto [sin1, cos1] = amrex::Math::sincos(nu_constants::fac1 * tau);
@@ -902,7 +905,7 @@ void nu_photo(const sneutf_t& sf,
     amrex::Real xnum = dum * z;
     amrex::Real xnumdt, xnumda, xnumdz;
     if constexpr (do_derivatives) {
-        xnumdt = dumdt * z - dum * z * cc * sf.zetadt;
+        xnumdt = dumdt * z - dum * z * (cc * sf.zetadt + sf.zeta * ccdt);
         xnumda = dumda * z - dum * z * cc * sf.zetada;
         xnumdz = dumdz * z - dum * z * cc * sf.zetadz;
     }
@@ -1383,7 +1386,7 @@ void nu_recomb(const sneutf_t& sf,
         }
 
         amrex::Real c00 = 1.0e0_rt / (1.0e0_rt + f1 * nu + f2 * nu2 + f3 * nu3);
-        amrex::Real c01 = f1 + f2 * 2.0e0_rt * nu + f3 * 3.0e0_rt * nu2;
+        amrex::Real c01 = -(f1 + f2 * 2.0e0_rt * nu + f3 * 3.0e0_rt * nu2) * c00 * c00;
         amrex::Real dum = zeta * c00;
         amrex::Real dumdt, dumda, dumdz;
         if constexpr (do_derivatives) {

diff --git a/unit_test/test_rhs/ci-benchmarks/ecsn.out b/unit_test/test_rhs/ci-benchmarks/ecsn.out
@@ -157,7 +157,7 @@
  J_silicon-28_sulfur-32                      0                      0
  J_phosphorus-31_sulfur-32                      0                      0
  J_sulfur-32_sulfur-32                       0                      0
- J_E_sulfur-32               -1.2066886642e+15       7.8006383744e+14
+ J_E_sulfur-32               -1.2066886642e+15       7.8006373177e+14
  J_hydrogen-1_E                 -0.12580793235       1.2089803642e+12
  J_helium-4_E                 -0.0078065301852       1.4905283338e+12
  J_oxygen-16_E               -5.3990173959e+12       1.2739583515e-06