few more details chainged to address suggestions

lanl · Mar 20, 2024 · 7f3bc17 · 7f3bc17
1 parent 994e2cc
commit 7f3bc17
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 12 deletions.
diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp
@@ -309,6 +309,7 @@ Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain) {
   const bool do_gain = pdo_gain_reducer->val;
   auto progenitor = pmb->packages.Get("progenitor").get();
   const Real inside_pns_threshold = progenitor->Param<Real>("inside_pns_threshold");
+  const Real net_heat_threshold = progenitor->Param<Real>("net_heat_threshold");
 
   parthenon::par_reduce(
       parthenon::LoopPatternMDRange(), "Calculates mass accretion rate (SN)",
@@ -322,9 +323,9 @@ Real CalculateMdot(MeshData<Real> *md, Real rc, bool gain) {
         Real r = std::sqrt(C[1] * C[1] + C[2] * C[2] + C[3] * C[3]);
         if (r <= rc) {
           if (gain) {
-            bool is_netheat = (v(b, internal_variables::GcovHeat(), k, j, i) -
-                                   v(b, internal_variables::GcovCool(), k, j, i) >
-                               1.e-8); // checks that in the gain region
+            bool is_netheat = ((v(b, internal_variables::GcovHeat(), k, j, i) -
+                                v(b, internal_variables::GcovCool(), k, j, i)) >
+                               net_heat_threshold); // checks that in the gain region
             bool is_inside_pns = (r < inside_pns_threshold); // checks that inside PNS
             if (do_gain && (is_inside_pns || is_netheat)) {
 

diff --git a/src/analysis/history.hpp b/src/analysis/history.hpp
@@ -110,15 +110,17 @@ Real ReduceInGain(MeshData<Real> *md, bool is_conserved, int idx = 0) {
   auto geom = Geometry::GetCoordinateSystem(md);
   auto progenitor = pmb->packages.Get("progenitor").get();
   const Real outside_pns_threshold = progenitor->Param<Real>("outside_pns_threshold");
+  const Real net_heat_threshold = progenitor->Param<Real>("net_heat_threshold");
 
   parthenon::par_reduce(
       parthenon::LoopPatternMDRange(),
       "Phoebus History for integrals in gain region (SN)", DevExecSpace(), 0, nblocks - 1,
       kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lresult) {
         Real gdet = geom.DetGamma(CellLocation::Cent, 0, k, j, i);
-        bool is_netheat = (v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i) >
-                           1.e-8); // checks that in the gain region
+        bool is_netheat =
+            ((v(b, iv::GcovHeat(), k, j, i) - v(b, iv::GcovCool(), k, j, i)) >
+             net_heat_threshold); // checks that in the gain region
         bool is_outside_pns = (v(b, fluid_prim::entropy(), k, j, i) >
                                outside_pns_threshold); // checks that outside PNS
         const auto &coords = v.GetCoordinates(b);

diff --git a/src/progenitor/progenitordata.cpp b/src/progenitor/progenitordata.cpp
@@ -90,11 +90,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
                         2.42e-5); // corresponds to entropy > 3 kb/baryon
   Real inside_pns_threshold = pin->GetOrAddReal("progenitor", "inside_pns_threshold",
                                                 0.008); // corresponds to r < 80 km
+  Real net_heat_threshold = pin->GetOrAddReal("progenitor", "inside_pns_threshold",
+                                              1e-8); // corresponds to r < 80 km
   UnitConversions units(pin);
   Real LengthCGSToCode = units.GetLengthCGSToCode();
-  Real rc = 400e5 * LengthCGSToCode;
   auto mdot_radii = pin->GetOrAddVector("progenitor", "mdot_radii",
-                                        std::vector<Real>{rc}); // default 400km
+                                        std::vector<Real>{400}); // default 400km
 
   // Add Params
   params.Add("mass_density", mass_density);
@@ -123,6 +124,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
   params.Add("outside_pns_threshold", outside_pns_threshold);
   params.Add("inside_pns_threshold", inside_pns_threshold);
+  params.Add("net_heat_threshold", net_heat_threshold);
   params.Add("mdot_radii", mdot_radii);
 
   // Reductions
@@ -140,12 +142,12 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
            ReduceInGain<internal_variables::GcovCool>(md, 0, 0);
   };
   for (auto rc : mdot_radii) {
-    auto Mdot = [rc, LengthCGSToCode](MeshData<Real> *md) {
-      return History::CalculateMdot(md, rc, 0);
+    auto rc_code = rc * 1e5 * LengthCGSToCode;
+    auto Mdot = [rc_code](MeshData<Real> *md) {
+      return History::CalculateMdot(md, rc_code, 0);
     };
-    hst_vars.emplace_back(HistoryOutputVar(
-        HstSum, Mdot,
-        "Mdot at r = " + std::to_string(int(rc / LengthCGSToCode / 1.e5)) + "km"));
+    hst_vars.emplace_back(
+        HistoryOutputVar(HstSum, Mdot, "Mdot at r = " + std::to_string(int(rc)) + "km"));
   }
 
   Real x1max = pin->GetReal("parthenon/mesh", "x1max");