feat(neutrino): Updated neutrino output

GridFire now reports neutrino loss for reaclib reactions. Note this
currently is only computed if precomputation is enabled.

src/extern/fortran/gridfire_mod.f90

@@ -102,7 +102,7 @@ module gridfire_mod
         end function
 
         ! int gf_evolve(...)
-        function gf_evolve(ptr, Y_in, num_species, T, rho, dt, Y_out, energy_out, dEps_dT, dEps_dRho, mass_lost) result(ierr) &
+        function gf_evolve(ptr, Y_in, num_species, T, rho, dt, Y_out, energy_out, dEps_dT, dEps_dRho, specific_neutrino_loss, specific_neutrino_flux, mass_lost) result(ierr) &
                 bind(C, name="gf_evolve")
             import
             type(c_ptr), value :: ptr
@@ -110,7 +110,7 @@ module gridfire_mod
             integer(c_size_t), value :: num_species
             real(c_double), value :: T, rho, dt
             real(c_double), dimension(*), intent(out) :: Y_out
-            real(c_double), intent(out) :: energy_out, dEps_dT, dEps_dRho, mass_lost
+            real(c_double), intent(out) :: energy_out, dEps_dT, dEps_dRho, specific_neutrino_loss, specific_neutrino_flux, mass_lost
             integer(c_int) :: ierr
         end function
     end interface
@@ -235,12 +235,12 @@ module gridfire_mod
             end if
         end subroutine setup_solver
 
-        subroutine evolve(self, Y_in, T, rho, dt, Y_out, energy, dedt, dedrho, mass_lost, ierr)
+        subroutine evolve(self, Y_in, T, rho, dt, Y_out, energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost, ierr)
             class(GridFire), intent(in) :: self
             real(c_double), dimension(:), intent(in) :: Y_in
             real(c_double), value :: T, rho, dt
             real(c_double), dimension(:), intent(out) :: Y_out
-            real(c_double), intent(out) :: energy, dedt, dedrho, mass_lost
+            real(c_double), intent(out) :: energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost
             integer, intent(out) :: ierr
             integer(c_int) :: c_ierr
 
@@ -248,7 +248,7 @@ module gridfire_mod
                     Y_in, self%num_species, &
                     T, rho, dt, &
                     Y_out, &
-                    energy, dedt, dedrho, mass_lost)
+                    energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost)
 
             ierr = int(c_ierr)
             if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then

src/extern/include/gridfire/extern/gridfire_context.h

@@ -31,7 +31,10 @@ struct GridFireContext {
         double* Y_out,
         double& energy_out,
         double& dEps_dT,
-        double& dEps_dRho, double& mass_lost
+        double& dEps_dRho,
+        double& specific_neutrino_energy_loss,
+        double& specific_neutrino_flux,
+        double& mass_lost
     );
 
     std::string last_error;

src/extern/include/gridfire/extern/gridfire_extern.h

@@ -81,7 +81,10 @@ extern "C" {
         double* Y_out,
         double* energy_out,
         double* dEps_dT,
-        double* dEps_dRho, double* mass_lost
+        double* dEps_dRho,
+        double* specific_neutrino_energy_loss,
+        double* specific_neutrino_flux,
+        double* mass_lost
     );
 
 #ifdef __cplusplus

src/extern/lib/gridfire_context.cpp

@@ -121,7 +121,10 @@ int GridFireContext::evolve(
     double* Y_out,
     double& energy_out,
     double& dEps_dT,
-    double& dEps_dRho, double& mass_lost
+    double& dEps_dRho,
+    double& specific_neutrino_energy_loss,
+    double& specific_neutrino_flux,
+    double& mass_lost
 ) {
     init_composition_from_abundance_vector(Y_in, num_species);
 
@@ -137,6 +140,8 @@ int GridFireContext::evolve(
     energy_out = result.energy;
     dEps_dT  = result.dEps_dT;
     dEps_dRho = result.dEps_dRho;
+    specific_neutrino_energy_loss = result.specific_neutrino_energy_loss;
+    specific_neutrino_flux = result.specific_neutrino_flux;
 
     std::set<fourdst::atomic::Species> seen_species;
     for (size_t i = 0; i < num_species; i++) {

src/extern/lib/gridfire_extern.cpp

@@ -86,7 +86,7 @@ extern "C" {
 
     int gf_evolve(
         void* ptr,
-        const double* Y,
+        const double* Y_in,
         const size_t num_species,
         const double T,
         const double rho,
@@ -95,12 +95,15 @@ extern "C" {
         double* Y_out,
         double* energy_out,
         double* dEps_dT,
-        double* dEps_dRho, double* mass_lost
+        double* dEps_dRho,
+        double* specific_neutrino_energy_loss,
+        double* specific_neutrino_flux,
+        double* mass_lost
     ) {
         auto* ctx = static_cast<GridFireContext*>(ptr);
         try {
             const int result = ctx->evolve(
-                Y,
+                Y_in,
                 num_species,
                 T,
                 rho,
@@ -109,7 +112,10 @@ extern "C" {
                 Y_out,
                 *energy_out,
                 *dEps_dT,
-                *dEps_dRho, *mass_lost
+                *dEps_dRho,
+                *specific_neutrino_energy_loss,
+                *specific_neutrino_flux,
+                *mass_lost
             );
             if (result != 0) {
                 return result;

src/include/gridfire/reaction/reaclib.h

@@ -15,13 +15,4 @@ namespace gridfire::reaclib {
      * @return A constant reference to the application-wide reaction set.
      */
     const reaction::ReactionSet &get_all_reaclib_reactions();
-
-    // Simple heuristic to check if a reaclib reaction is a strong or weak reaction
-    /*  A weak reaction is defined here as one where:
-        - The number of reactants is equal to the number of products
-        - There is only one reactant and one product
-        - The mass number (A) of the reactant is equal to the mass number (A) of the product
-    */
-    bool reaction_is_weak(const reaction::Reaction& reaction);
-
 } // namespace gridfire::reaclib

src/include/gridfire/reaction/reaction.h

@@ -10,7 +10,6 @@
 #include <vector>
 #include <unordered_set>
 
-
 #include "cppad/cppad.hpp"
 #include "fourdst/composition/composition.h"
 
@@ -48,6 +47,7 @@ namespace gridfire::reaction {
         {ReactionType::REACLIB_WEAK, "Weak"},
         {ReactionType::LOGICAL_REACLIB_WEAK, "Weak"},
     };
+
     /**
      * @struct RateCoefficientSet
      * @brief Holds the seven coefficients for the REACLIB rate equation.
@@ -358,6 +358,15 @@ namespace gridfire::reaction {
         [[nodiscard]] virtual std::optional<std::vector<RateCoefficientSet>> getRateCoefficients() const = 0;
 
     };
+
+    // Simple heuristic to check if a reaclib reaction is a strong or weak reaction
+    /*  A weak reaction is defined here as one where:
+        - The number of reactants is equal to the number of products
+        - There is only one reactant and one product
+        - The mass number (A) of the reactant is equal to the mass number (A) of the product
+    */
+    bool reaction_is_weak(const Reaction& reaction);
+
     class ReaclibReaction : public Reaction {
     public:
         ~ReaclibReaction() override = default;
@@ -448,7 +457,12 @@ namespace gridfire::reaction {
          */
         [[nodiscard]] std::string_view sourceLabel() const { return m_sourceLabel; }
 
-        [[nodiscard]] ReactionType type() const override { return ReactionType::REACLIB; }
+        [[nodiscard]] ReactionType type() const override {
+            if (reaction::reaction_is_weak(*this)) {
+                return ReactionType::REACLIB_WEAK;
+            }
+            return ReactionType::REACLIB;
+        }
 
         /**
          * @brief Gets the set of rate coefficients.

src/include/gridfire/types/types.h

@@ -40,8 +40,8 @@ namespace gridfire {
         double energy; ///< Energy in ergs after evaluation
         double dEps_dT; ///< Partial derivative of energy generation rate with respect to temperature
         double dEps_dRho; ///< Partial derivative of energy generation rate with respect to density
-        double neutrino_energy_loss_rate; ///< Neutrino energy loss rate in ergs/g/s
+        double specific_neutrino_energy_loss; ///< Neutrino energy loss rate in ergs/g/s
-        double total_neutrino_flux; ///< Total neutrino flux in neutrinos/g/s
+        double specific_neutrino_flux; ///< Total neutrino flux in neutrinos/g/s
 
         friend std::ostream& operator<<(std::ostream& os, const NetOut& netOut) {
             os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", ε=" << netOut.energy << ", dε/dT=" << netOut.dEps_dT << ", dε/dρ=" << netOut.dEps_dRho << ")";

src/lib/engine/procedures/construction.cpp

@@ -88,7 +88,7 @@ namespace {
         if (has_flag(reaction_types, gridfire::engine::NetworkConstructionFlags::REACLIB_STRONG)) {
             const auto& allReaclibReactions = gridfire::reaclib::get_all_reaclib_reactions();
             for (const auto& reaction : allReaclibReactions) {
-                const bool isWeakReaction = gridfire::reaclib::reaction_is_weak(*reaction);
+                const bool isWeakReaction = gridfire::reaction::reaction_is_weak(*reaction);
                 const bool okayToUseReaclibWeakReaction = has_flag(reaction_types, gridfire::engine::NetworkConstructionFlags::REACLIB_WEAK);
 
                 const bool reaclibWeakOkay = !isWeakReaction || okayToUseReaclibWeakReaction;

src/lib/reaction/reaclib.cpp

@@ -224,22 +224,4 @@ namespace gridfire::reaclib {
         return *s_all_reaclib_reactions_ptr;
     }
 
-    bool reaction_is_weak(const reaction::Reaction& reaction) {
-        const std::vector<fourdst::atomic::Species>& reactants = reaction.reactants();
-        const std::vector<fourdst::atomic::Species>& products = reaction.products();
-
-        if (reactants.size() != products.size()) {
-            return false;
-        }
-
-        if (reactants.size() != 1 || products.size() != 1) {
-            return false;
-        }
-
-        if (std::floor(reactants[0].a()) != std::floor(products[0].a())) {
-            return false;
-        }
-
-        return true;
-    }
 } // namespace gridfire::reaclib

src/lib/reaction/reaction.cpp

@@ -231,7 +231,7 @@ namespace gridfire::reaction {
 
         std::unordered_set<bool> reaction_weak_types;
         for (const auto& reaction : reactions) {
-            reaction_weak_types.insert(reaclib::reaction_is_weak(*reaction));
+            reaction_weak_types.insert(reaction::reaction_is_weak(*reaction));
         }
 
         if (reaction_weak_types.size() != 1) {
@@ -611,8 +611,28 @@ namespace gridfire::reaction {
             }
         }
 
-    return finalReactionSet;
+        return finalReactionSet;
-}
+    }
+
+    bool reaction_is_weak(const reaction::Reaction& reaction) {
+        const std::vector<fourdst::atomic::Species>& reactants = reaction.reactants();
+        const std::vector<fourdst::atomic::Species>& products = reaction.products();
+
+        if (reactants.size() != products.size()) {
+            return false;
+        }
+
+        if (reactants.size() != 1 || products.size() != 1) {
+            return false;
+        }
+
+        if (std::floor(reactants[0].a()) != std::floor(products[0].a())) {
+            return false;
+        }
+
+        return true;
+    }
+
 
 }
 

src/lib/solver/strategies/CVODE_solver_strategy.cpp

@@ -190,8 +190,11 @@ namespace gridfire::solver {
 
             sunrealtype* y_data = N_VGetArrayPointer(m_Y);
             const double current_energy = y_data[numSpecies]; // Specific energy rate
-            accumulated_neutrino_energy_loss += user_data.neutrino_energy_loss_rate;
+
-            accumulated_total_neutrino_flux += user_data.total_neutrino_flux;
+            // TODO: Accumulate neutrino loss through the state vector directly which will allow CVODE to properly integrate it
+            accumulated_neutrino_energy_loss += user_data.neutrino_energy_loss_rate * last_step_size;
+            accumulated_total_neutrino_flux += user_data.total_neutrino_flux * last_step_size;
+
             size_t iter_diff = (total_nonlinear_iterations + nliters) - prev_nonlinear_iterations;
             size_t convFail_diff = (total_convergence_failures + nlcfails) - prev_convergence_failures;
             if (m_stdout_logging_enabled) {
@@ -508,8 +511,8 @@ namespace gridfire::solver {
 
         netOut.dEps_dT = dEps_dT;
         netOut.dEps_dRho = dEps_dRho;
-        netOut.neutrino_energy_loss_rate = accumulated_neutrino_energy_loss;
+        netOut.specific_neutrino_energy_loss = accumulated_neutrino_energy_loss;
-        netOut.total_neutrino_flux = accumulated_total_neutrino_flux;
+        netOut.specific_neutrino_flux = accumulated_total_neutrino_flux;
         LOG_TRACE_L2(m_logger, "Output data built!");
         LOG_TRACE_L2(m_logger, "Solver evaluation complete!.");
 

src/python/types/bindings.cpp

@@ -35,6 +35,8 @@ void register_type_bindings(const pybind11::module &m) {
         .def_readonly("energy", &gridfire::NetOut::energy)
         .def_readonly("dEps_dT", &gridfire::NetOut::dEps_dT)
         .def_readonly("dEps_dRho", &gridfire::NetOut::dEps_dRho)
+        .def_readonly("specific_neutrino_flux", &gridfire::NetOut::specific_neutrino_flux)
+        .def_readonly("specific_neutrino_energy_loss", &gridfire::NetOut::specific_neutrino_energy_loss)
         .def("__repr__", [](const gridfire::NetOut &netOut) {
             std::string repr = std::format(
                 "NetOut(<μ> = {} steps = {}, ε = {}, dε/dT = {}, dε/dρ = {})",

tests/extern/C/gridfire_evolve.c

@@ -46,6 +46,8 @@ int main() {
     double energy_out;
     double dEps_dT;
     double dEps_dRho;
+    double specific_neutrino_energy_loss;
+    double specific_neutrino_flux;
     double mass_lost;
 
     ret = gf_evolve(
@@ -59,7 +61,10 @@ int main() {
         Y_out,
         &energy_out,
         &dEps_dT,
-        &dEps_dRho, &mass_lost
+        &dEps_dRho,
+        &specific_neutrino_energy_loss,
+        &specific_neutrino_flux,
+        &mass_lost
     );
 
     CHECK_RET_CODE(ret, ctx, "EVOLUTION");
@@ -72,6 +77,8 @@ int main() {
     printf("Energy output: %e\n", energy_out);
     printf("dEps/dT: %e\n", dEps_dT);
     printf("dEps/dRho: %e\n", dEps_dRho);
+    printf("Specific neutrino energy loss: %e\n", specific_neutrino_energy_loss);
+    printf("Specific neutrino flux: %e\n", specific_neutrino_flux);
     printf("Mass lost: %e\n", mass_lost);
 
     gf_free(ctx);

tests/extern/fortran/gridfire_evolve.f90

@@ -36,7 +36,7 @@ program main
 
     ! Output buffers
     real(c_double), dimension(8) :: Y_out
-    real(c_double) :: energy_out, dedt, dedrho, dmass
+    real(c_double) :: energy_out, dedt, dedrho, snu_e_loss, snu_flux, dmass
 
     ! Thermodynamic Conditions (Solar Core-ish)
     real(c_double) :: T = 1.5e7      ! 15 Million K
@@ -60,7 +60,7 @@ program main
 
     ! --- 5. Evolve ---
     print *, "Evolving system (dt =", dt, "s)..."
-    call net%evolve(Y_in, T, rho, dt, Y_out, energy_out, dedt, dedrho, dmass, ierr)
+    call net%evolve(Y_in, T, rho, dt, Y_out, energy_out, dedt, dedrho, snu_e_loss, snu_flux, dmass, ierr)
 
     if (ierr /= 0) then
         print *, "Evolution Failed with error code: ", ierr
@@ -74,6 +74,9 @@ program main
     print *, "--- Results ---"
     print '(A, ES12.5, A)', "Energy Generation: ", energy_out, " erg/g/s"
     print '(A, ES12.5)',    "dEps/dT:           ", dedt
+    print '(A, ES12.5)',    "dEps/drho:         ", dedrho
+    print '(A, ES12.5)',    "Neutrino Energy Loss:      ", snu_e_loss
+    print '(A, ES12.5)',    "Neutrino Flux:          ", snu_flux
     print '(A, ES12.5)',    "Mass Change:       ", dmass
 
     print *, ""

tests/graphnet_sandbox/main.cpp

@@ -96,13 +96,13 @@ void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) {
     fractional.push_back(0.0);
 
     initial.push_back(0.0);
-    final.push_back(netOut.neutrino_energy_loss_rate);
+    final.push_back(netOut.specific_neutrino_energy_loss);
-    delta.push_back(netOut.neutrino_energy_loss_rate);
+    delta.push_back(netOut.specific_neutrino_energy_loss);
     fractional.push_back(0.0);
 
     initial.push_back(0.0);
-    final.push_back(netOut.total_neutrino_flux);
+    final.push_back(netOut.specific_neutrino_flux);
-    delta.push_back(netOut.total_neutrino_flux);
+    delta.push_back(netOut.specific_neutrino_flux);
     fractional.push_back(0.0);
 
     initial.push_back(netIn.composition.getMeanParticleMass());
@@ -282,10 +282,6 @@ int main(int argc, char** argv) {
     CLI11_PARSE(app, argc, argv);
 
     const NetIn netIn = init(temp, rho, tMax);
-    std::println("Starting Integration with T = {} K, ρ = {} g/cm³, tMax = {} s", temp, rho, tMax);
-    std::println("Composition is: {}", utils::iterable_to_delimited_string(netIn.composition, ", ", [&netIn](std::pair<const fourdst::atomic::Species&, double> arg) {
-        return std::format("{:5}: {:10.4E}", arg.first.name(), arg.second);
-    }));
 
     policy::MainSequencePolicy stellarPolicy(netIn.composition);
     stellarPolicy.construct();