diff --git a/src/extern/fortran/gridfire_mod.f90 b/src/extern/fortran/gridfire_mod.f90 index fa02f1ce..1c4c3253 100644 --- a/src/extern/fortran/gridfire_mod.f90 +++ b/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 diff --git a/src/extern/include/gridfire/extern/gridfire_context.h b/src/extern/include/gridfire/extern/gridfire_context.h index 5c432820..a33904e6 100644 --- a/src/extern/include/gridfire/extern/gridfire_context.h +++ b/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; diff --git a/src/extern/include/gridfire/extern/gridfire_extern.h b/src/extern/include/gridfire/extern/gridfire_extern.h index 7ae2edb5..7fa9e753 100644 --- a/src/extern/include/gridfire/extern/gridfire_extern.h +++ b/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 diff --git a/src/extern/lib/gridfire_context.cpp b/src/extern/lib/gridfire_context.cpp index 4df5c68f..41f81fc1 100644 --- a/src/extern/lib/gridfire_context.cpp +++ b/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 seen_species; for (size_t i = 0; i < num_species; i++) { diff --git a/src/extern/lib/gridfire_extern.cpp b/src/extern/lib/gridfire_extern.cpp index a8c7fe7d..54aaef22 100644 --- a/src/extern/lib/gridfire_extern.cpp +++ b/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(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; diff --git a/src/include/gridfire/reaction/reaclib.h b/src/include/gridfire/reaction/reaclib.h index b44183a8..199ca5cd 100644 --- a/src/include/gridfire/reaction/reaclib.h +++ b/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 diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h index 35dd23fd..d60d0281 100644 --- a/src/include/gridfire/reaction/reaction.h +++ b/src/include/gridfire/reaction/reaction.h @@ -10,7 +10,6 @@ #include #include - #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> 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. diff --git a/src/include/gridfire/types/types.h b/src/include/gridfire/types/types.h index d8b2a1e9..64ef6f3a 100644 --- a/src/include/gridfire/types/types.h +++ b/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 total_neutrino_flux; ///< Total neutrino flux in neutrinos/g/s + double specific_neutrino_energy_loss; ///< Neutrino energy loss rate in ergs/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 << ")"; diff --git a/src/lib/engine/procedures/construction.cpp b/src/lib/engine/procedures/construction.cpp index dee405d1..2e98fd31 100644 --- a/src/lib/engine/procedures/construction.cpp +++ b/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; diff --git a/src/lib/reaction/reaclib.cpp b/src/lib/reaction/reaclib.cpp index efbc930b..0157e353 100644 --- a/src/lib/reaction/reaclib.cpp +++ b/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& reactants = reaction.reactants(); - const std::vector& 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 \ No newline at end of file diff --git a/src/lib/reaction/reaction.cpp b/src/lib/reaction/reaction.cpp index 94211d7b..debee6f6 100644 --- a/src/lib/reaction/reaction.cpp +++ b/src/lib/reaction/reaction.cpp @@ -231,7 +231,7 @@ namespace gridfire::reaction { std::unordered_set 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& reactants = reaction.reactants(); + const std::vector& 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; + } + } diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index 5fa97af6..7d430481 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/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.total_neutrino_flux = accumulated_total_neutrino_flux; + netOut.specific_neutrino_energy_loss = accumulated_neutrino_energy_loss; + netOut.specific_neutrino_flux = accumulated_total_neutrino_flux; LOG_TRACE_L2(m_logger, "Output data built!"); LOG_TRACE_L2(m_logger, "Solver evaluation complete!."); diff --git a/src/python/types/bindings.cpp b/src/python/types/bindings.cpp index 5e611c69..6938234a 100644 --- a/src/python/types/bindings.cpp +++ b/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ρ = {})", diff --git a/tests/extern/C/gridfire_evolve.c b/tests/extern/C/gridfire_evolve.c index 9055f601..7834dabf 100644 --- a/tests/extern/C/gridfire_evolve.c +++ b/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); diff --git a/tests/extern/fortran/gridfire_evolve.f90 b/tests/extern/fortran/gridfire_evolve.f90 index 592a01f9..ac0b9c4d 100644 --- a/tests/extern/fortran/gridfire_evolve.f90 +++ b/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 *, "" diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index 874e54b3..4455606f 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/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); - delta.push_back(netOut.neutrino_energy_loss_rate); + final.push_back(netOut.specific_neutrino_energy_loss); + 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); - delta.push_back(netOut.total_neutrino_flux); + final.push_back(netOut.specific_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 arg) { - return std::format("{:5}: {:10.4E}", arg.first.name(), arg.second); - })); policy::MainSequencePolicy stellarPolicy(netIn.composition); stellarPolicy.construct();