diff --git a/src/extern/include/gridfire/extern/gridfire_context.h b/src/extern/include/gridfire/extern/gridfire_context.h index a33904e6..72fc2651 100644 --- a/src/extern/include/gridfire/extern/gridfire_context.h +++ b/src/extern/include/gridfire/extern/gridfire_context.h @@ -7,19 +7,32 @@ #include #include -struct GridFireContext { +enum class GFContextType { + POINT, + GRID +}; + +struct GFContext { + virtual ~GFContext() = default; + std::unique_ptr policy; - gridfire::engine::DynamicEngine* engine; - std::unique_ptr solver; - + const gridfire::engine::DynamicEngine* engine; + std::unique_ptr engine_ctx; std::vector speciesList; - fourdst::composition::Composition working_comp; - void init_species_map(const std::vector& species_names); - void init_engine_from_policy(const std::string& policy_name, const double *abundances, size_t num_species); - void init_solver_from_engine(const std::string& solver_name); + virtual void init_species_map(const std::vector& species_names); + virtual void init_engine_from_policy(const std::string& policy_name, const double *abundances, size_t num_species); + virtual void init_solver_from_engine() = 0; - void init_composition_from_abundance_vector(const double* abundances, size_t num_species); + fourdst::composition::Composition init_composition_from_abundance_vector(const std::vector &abundances, size_t num_species) const; + std::string last_error; +}; + +struct GFPointContext final: GFContext{ + std::unique_ptr solver; + std::unique_ptr solver_ctx; + + void init_solver_from_engine() override; int evolve( const double* Y_in, @@ -35,9 +48,45 @@ struct GridFireContext { double& specific_neutrino_energy_loss, double& specific_neutrino_flux, double& mass_lost - ); + ) const; - std::string last_error; }; +struct GFGridContext final : GFContext { + std::unique_ptr local_solver; + std::unique_ptr solver; + std::unique_ptr solver_ctx; + + void init_solver_from_engine() override; + + size_t zones; + + void set_zones(const size_t num_zones) { + zones = num_zones; + } + + [[nodiscard]] size_t get_zones() const { + return zones; + } + + int evolve( + const double* Y_in, + size_t num_species, + const double* T, + const double* rho, + double tMax, + double dt0, + double* Y_out, + double* energy_out, + double* dEps_dT, + double* dEps_dRho, + double* specific_neutrino_energy_loss, + double* specific_neutrino_flux, + double* mass_lost + ) const; + +}; + +std::unique_ptr make_gf_context(const GFContextType& type); + #endif \ No newline at end of file diff --git a/src/extern/include/gridfire/extern/gridfire_extern.h b/src/extern/include/gridfire/extern/gridfire_extern.h index 7fa9e753..6799f7cf 100644 --- a/src/extern/include/gridfire/extern/gridfire_extern.h +++ b/src/extern/include/gridfire/extern/gridfire_extern.h @@ -6,6 +6,12 @@ #ifdef __cplusplus extern "C" { #endif + enum GF_TYPE { + SINGLE_ZONE = 0, + MULTI_ZONE = 1 + }; + + enum FDSSE_ERROR_CODES { FDSSE_NON_4DSTAR_ERROR = -102, FDSSE_UNKNOWN_ERROR = -101, @@ -54,37 +60,49 @@ extern "C" { GF_DEBUG_ERROR = 30, GF_GRIDFIRE_ERROR = 31, + GF_UNINITIALIZED_INPUT_MEMORY_ERROR = 32, + GF_UNINITIALIZED_OUTPUT_MEMORY_ERROR = 33, + + GF_INVALID_NUM_SPECIES = 34, + + GF_INVALID_TIMESTEPS = 35, + GF_UNKNOWN_FREE_TYPE = 36, + + GF_INVALID_TYPE = 37, }; char* gf_get_last_error_message(void* ptr); char* gf_error_code_to_string(int error_code); - void* gf_init(); + void* gf_init(const enum GF_TYPE type); - void gf_free(void* ctx); + int gf_free(const enum GF_TYPE type, void *ctx); + + int gf_set_num_zones(const enum GF_TYPE type, void* ptr, const size_t num_zones); int gf_register_species(void* ptr, const int num_species, const char** species_names); int gf_construct_engine_from_policy(void* ptr, const char* policy_name, const double *abundances, size_t num_species); - int gf_construct_solver_from_engine(void* ptr, const char* solver_name); + int gf_construct_solver_from_engine(void* ptr); int gf_evolve( + enum GF_TYPE type, void* ptr, - const double* Y_in, + const void* Y_in, size_t num_species, - double T, - double rho, + const void* T, + const void* rho, double tMax, double dt0, - double* Y_out, - double* energy_out, - double* dEps_dT, - double* dEps_dRho, - double* specific_neutrino_energy_loss, - double* specific_neutrino_flux, - double* mass_lost + void* Y_out, + void* energy_out, + void* dEps_dT, + void* dEps_dRho, + void* specific_neutrino_energy_loss, + void* specific_neutrino_flux, + void* mass_lost ); #ifdef __cplusplus diff --git a/src/extern/lib/gridfire_context.cpp b/src/extern/lib/gridfire_context.cpp index 41f81fc1..3c74f695 100644 --- a/src/extern/lib/gridfire_context.cpp +++ b/src/extern/lib/gridfire_context.cpp @@ -3,11 +3,10 @@ #include "fourdst/atomic/species.h" #include "fourdst/composition/exceptions/exceptions_composition.h" +#include "gridfire/exceptions/error_policy.h" +#include "gridfire/utils/logging.h" -void GridFireContext::init_species_map(const std::vector &species_names) { - for (const auto& name: species_names) { - working_comp.registerSymbol(name); - } +void GFContext::init_species_map(const std::vector &species_names) { this->speciesList.clear(); this->speciesList.reserve(species_names.size()); @@ -24,8 +23,9 @@ void GridFireContext::init_species_map(const std::vector &species_n } -void GridFireContext::init_engine_from_policy(const std::string &policy_name, const double *abundances, const size_t num_species) { - init_composition_from_abundance_vector(abundances, num_species); +void GFContext::init_engine_from_policy(const std::string &policy_name, const double *abundances, const size_t num_species) { + const std::vector Y_scratch(abundances, abundances + num_species); + fourdst::composition::Composition comp = init_composition_from_abundance_vector(Y_scratch, num_species); enum class EnginePolicy { MAIN_SEQUENCE_POLICY @@ -47,71 +47,53 @@ void GridFireContext::init_engine_from_policy(const std::string &policy_name, co switch (engine_map.at(policy_name)) { case EnginePolicy::MAIN_SEQUENCE_POLICY: { - this->policy = std::make_unique(this->working_comp); - this->engine = &policy->construct(); + this->policy = std::make_unique(comp); + const auto& [e, ctx] = policy->construct(); + this->engine = &e; + this->engine_ctx = ctx->clone_structure(); + break; } default: throw gridfire::exceptions::PolicyError( - "Unhandled engine policy in GridFireContext::init_engine_from_policy" + "Unhandled engine policy in GFPointContext::init_engine_from_policy" ); } } -void GridFireContext::init_solver_from_engine(const std::string &solver_name) { - enum class SolverType { - CVODE - }; - - static const std::unordered_map solver_map = { - {"CVODE", SolverType::CVODE} - }; - - if (!solver_map.contains(solver_name)) { - throw gridfire::exceptions::SolverError( - std::format( - "Solver {} is not recognized. Valid solvers are: {}", - solver_name, - gridfire::utils::iterable_to_delimited_string(solver_map, ", ", [](const auto& pair){ return pair.first; }) - ) - ); - } - - switch (solver_map.at(solver_name)) { - case SolverType::CVODE: { - this->solver = std::make_unique(*this->engine); - break; - } - default: - throw gridfire::exceptions::SolverError( - "Unhandled solver type in GridFireContext::init_solver_from_engine" - ); - } - -} - -void GridFireContext::init_composition_from_abundance_vector(const double *abundances, size_t num_species) { +fourdst::composition::Composition GFContext::init_composition_from_abundance_vector(const std::vector &abundances, size_t num_species) const { if (num_species == 0) { throw fourdst::composition::exceptions::InvalidCompositionError("Cannot initialize composition with zero species."); } - if (num_species != working_comp.size()) { + if (num_species != speciesList.size()) { throw fourdst::composition::exceptions::InvalidCompositionError( std::format( "Number of species provided ({}) does not match the registered species count ({}).", num_species, - working_comp.size() + speciesList.size() ) ); } + fourdst::composition::Composition comp; for (size_t i = 0; i < num_species; i++) { - this->working_comp.setMolarAbundance(this->speciesList[i], abundances[i]); + comp.registerSpecies(this->speciesList[i]); + comp.setMolarAbundance(this->speciesList[i], abundances[i]); } + + return comp; } -int GridFireContext::evolve( + +void GFPointContext::init_solver_from_engine() { + this->solver = std::make_unique(*this->engine); + this->solver_ctx = std::make_unique(*engine_ctx); +} + + +int GFPointContext::evolve( const double* Y_in, const size_t num_species, const double T, @@ -125,17 +107,19 @@ int GridFireContext::evolve( double& specific_neutrino_energy_loss, double& specific_neutrino_flux, double& mass_lost -) { - init_composition_from_abundance_vector(Y_in, num_species); +) const { + + const std::vector Y_scratch(Y_in, Y_in + num_species); + const fourdst::composition::Composition comp = init_composition_from_abundance_vector(Y_scratch, num_species); gridfire::NetIn netIn; netIn.temperature = T; netIn.density = rho; netIn.dt0 = dt0; netIn.tMax = tMax; - netIn.composition = this->working_comp; + netIn.composition = comp; - const gridfire::NetOut result = this->solver->evaluate(netIn); + const gridfire::NetOut result = this->solver->evaluate(*solver_ctx, netIn); energy_out = result.energy; dEps_dT = result.dEps_dT; @@ -158,4 +142,100 @@ int GridFireContext::evolve( } return 0; -} \ No newline at end of file +} + +void GFGridContext::init_solver_from_engine() { + this->local_solver = std::make_unique(*this->engine); + this->solver = std::make_unique(*this->engine, *this->local_solver); + this->solver_ctx = std::make_unique(*engine_ctx); +} + +int GFGridContext::evolve( + const double* Y_in, + size_t num_species, + const double *T, + const double *rho, + double tMax, + double dt0, + double *Y_out, + double *energy_out, + double *dEps_dT, + double *dEps_dRho, + double *specific_neutrino_energy_loss, + double *specific_neutrino_flux, + double *mass_lost +) const { + if (this->get_zones() == 0) { + throw gridfire::exceptions::GridFireError("GFGridContext has zero zones configured for evolution."); + } + + if (Y_in == nullptr || T == nullptr || rho == nullptr) { + throw gridfire::exceptions::GridFireError("Input abundance, temperature, or density pointers are null."); + } + + if (Y_out == nullptr) { + throw gridfire::exceptions::GridFireError("Output abundance pointer is null."); + } + + std::vector zone_compositions; + + zone_compositions.reserve(this->get_zones()); + + std::vector Y_scratch; + Y_scratch.resize(num_species); + + for (size_t i = 0; i < this->get_zones(); ++i) { + for (size_t j = 0; j < num_species; ++j) { + Y_scratch[j] = Y_in[i * num_species + j]; + } + zone_compositions.push_back(init_composition_from_abundance_vector(Y_scratch, num_species)); + } + + std::vector netIns; + netIns.reserve(this->get_zones()); + for (size_t i = 0; i < this->get_zones(); ++i) { + gridfire::NetIn netIn; + netIn.temperature = T[i]; + netIn.density = rho[i]; + netIn.dt0 = dt0; + netIn.tMax = tMax; + netIn.composition = zone_compositions[i]; + netIns.push_back(netIn); + } + + std::vector results = this->solver->evaluate(*this->solver_ctx, netIns); + + for (size_t zone_idx = 0; zone_idx < this->get_zones(); ++zone_idx) { + const gridfire::NetOut& netOut = results[zone_idx]; + energy_out[zone_idx] = netOut.energy; + dEps_dT[zone_idx] = netOut.dEps_dT;; + dEps_dRho[zone_idx] = netOut.dEps_dRho; + specific_neutrino_energy_loss[zone_idx] = netOut.specific_neutrino_energy_loss; + specific_neutrino_flux[zone_idx] = netOut.specific_neutrino_flux; + + std::set seen_species; + for (size_t i = 0; i < num_species; ++i) { + fourdst::atomic::Species species = this->speciesList[i]; + Y_out[zone_idx * num_species + i] = netOut.composition.getMolarAbundance(species); + seen_species.insert(species); + } + mass_lost[zone_idx] = 0.0; + for (const auto& species : netOut.composition.getRegisteredSpecies()) { + if (!seen_species.contains(species)) { + mass_lost[zone_idx] += species.mass() * netOut.composition.getMolarAbundance(species); + } + } + } + return 0; +} + +std::unique_ptr make_gf_context(const GFContextType &type) { + switch (type) { + case GFContextType::POINT: + return std::make_unique(); + case GFContextType::GRID: + return std::make_unique(); + default: + throw gridfire::exceptions::GridFireError("Unhandled GFContextType in make_gf_context"); + } +} diff --git a/src/extern/lib/gridfire_extern.cpp b/src/extern/lib/gridfire_extern.cpp index 54aaef22..bec417f4 100644 --- a/src/extern/lib/gridfire_extern.cpp +++ b/src/extern/lib/gridfire_extern.cpp @@ -3,120 +3,18 @@ #include "gridfire/extern/gridfire_context.h" #include "gridfire/extern/gridfire_extern.h" -extern "C" { - void* gf_init() { - return new GridFireContext(); - } +namespace { - void gf_free(void* ctx) { - delete static_cast(ctx); - } + template + concept ErrorTrackable = requires(T a) { + { a.last_error } -> std::convertible_to; + }; - int gf_register_species(void* ptr, const int num_species, const char** species_names) { - auto* ctx = static_cast(ptr); + template + int execute_guarded(Context* ctx, Func&& action) { try { - std::vector names; - for(int i=0; iinit_species_map(names); - return FDSSE_SUCCESS; - } catch (const fourdst::composition::exceptions::UnknownSymbolError& e) { - ctx->last_error = e.what(); - return FDSSE_UNKNOWN_SYMBOL_ERROR; - } catch (const fourdst::composition::exceptions::SpeciesError& e) { - ctx->last_error = e.what(); - return FDSSE_SPECIES_ERROR; - } catch (const std::exception& e) { - ctx->last_error = e.what(); - return FDSSE_NON_4DSTAR_ERROR; - } catch (...) { - ctx->last_error = "Unknown error occurred during species registration."; - return FDSSE_UNKNOWN_ERROR; - } - } + const int result = action(); - int gf_construct_engine_from_policy( - void* ptr, - const char* policy_name, - const double *abundances, - const size_t num_species - ) { - auto* ctx = static_cast(ptr); - try { - ctx->init_engine_from_policy(std::string(policy_name), abundances, num_species); - return GF_SUCCESS; - } catch (const gridfire::exceptions::MissingBaseReactionError& e) { - ctx->last_error = e.what(); - return GF_MISSING_BASE_REACTION_ERROR; - } catch (const gridfire::exceptions::MissingSeedSpeciesError& e) { - ctx->last_error = e.what(); - return GF_MISSING_SEED_SPECIES_ERROR; - } catch (const gridfire::exceptions::MissingKeyReactionError& e) { - ctx->last_error = e.what(); - return GF_MISSING_KEY_REACTION_ERROR; - } catch (const gridfire::exceptions::PolicyError& e) { - ctx->last_error = e.what(); - return GF_POLICY_ERROR; - } catch (std::exception& e) { - ctx->last_error = e.what(); - return GF_NON_GRIDFIRE_ERROR; - } catch (...) { - ctx->last_error = "Unknown error occurred during engine construction."; - return GF_UNKNOWN_ERROR; - } - } - - int gf_construct_solver_from_engine( - void* ptr, - const char* solver_name - ) { - auto* ctx = static_cast(ptr); - try { - ctx->init_solver_from_engine(std::string(solver_name)); - return GF_SUCCESS; - } catch (std::exception& e) { - ctx->last_error = e.what(); - return GF_NON_GRIDFIRE_ERROR; - } catch (...) { - ctx->last_error = "Unknown error occurred during solver construction."; - return GF_UNKNOWN_ERROR; - } - } - - int gf_evolve( - void* ptr, - const double* Y_in, - const size_t num_species, - const double T, - const double rho, - const double tMax, - const double dt0, - double* Y_out, - double* energy_out, - double* dEps_dT, - 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_in, - num_species, - T, - rho, - tMax, - dt0, - Y_out, - *energy_out, - *dEps_dT, - *dEps_dRho, - *specific_neutrino_energy_loss, - *specific_neutrino_flux, - *mass_lost - ); if (result != 0) { return result; } @@ -211,8 +109,218 @@ extern "C" { } } +} + +extern "C" { + + + void* gf_init(const enum GF_TYPE type) { + if (type == MULTI_ZONE) { + return new GFGridContext(); + } + if (type == SINGLE_ZONE) { + return new GFPointContext(); + } + return nullptr; + } + + int gf_free(const enum GF_TYPE type, void *ctx) { + if (!ctx) { + return GF_UNINITIALIZED_INPUT_MEMORY_ERROR; + } + if (type == MULTI_ZONE) { + delete static_cast(ctx); + return GF_SUCCESS; + } + if (type == SINGLE_ZONE) { + delete static_cast(ctx); + return GF_SUCCESS; + } + return GF_UNKNOWN_FREE_TYPE; + } + + int gf_set_num_zones(const enum GF_TYPE type, void* ptr, const size_t num_zones) { + if (type != MULTI_ZONE) { + return GF_INVALID_TYPE; + } + + if (!ptr) { + return GF_UNINITIALIZED_INPUT_MEMORY_ERROR; + } + + auto* ctx = static_cast(ptr); + return execute_guarded(ctx, [&]() { + ctx->set_zones(num_zones); + return GF_SUCCESS; + }); + } + + int gf_register_species(void* ptr, const int num_species, const char** species_names) { + if (num_species < 0) return GF_INVALID_NUM_SPECIES; + + if (num_species == 0) return GF_SUCCESS; + + if (!ptr || !species_names) { + return GF_UNINITIALIZED_INPUT_MEMORY_ERROR; + } + + for (int i=0; i < num_species; ++i) { + if (!species_names[i]) { + return GF_UNINITIALIZED_INPUT_MEMORY_ERROR; + } + } + + auto* ctx = static_cast(ptr); + return execute_guarded(ctx, [&]() { + std::vector names; + for (int i=0; iinit_species_map(names); + return FDSSE_SUCCESS; + }); + } + + int gf_construct_engine_from_policy( + void* ptr, + const char* policy_name, + const double *abundances, + const size_t num_species + ) { + auto* ctx = static_cast(ptr); + return execute_guarded(ctx, [&]() { + ctx->init_engine_from_policy(std::string(policy_name), abundances, num_species); + return GF_SUCCESS; + }); + } + + int gf_construct_solver_from_engine( + void* ptr + ) { + auto* ctx = static_cast(ptr); + return execute_guarded(ctx, [&]() { + ctx->init_solver_from_engine(); + return GF_SUCCESS; + }); + } + + int gf_evolve( + const enum GF_TYPE type, + void* ptr, + const void* Y_in, + const size_t num_species, + const void* T, + const void* rho, + const double tMax, + const double dt0, + void* Y_out, + void* energy_out, + void* dEps_dT, + void* dEps_dRho, + void* specific_neutrino_energy_loss, + void* specific_neutrino_flux, + void* mass_lost + ) { + + if (!ptr || !Y_in || !T || !rho) { + return GF_UNINITIALIZED_INPUT_MEMORY_ERROR; + } + + if (!Y_out || !energy_out || !dEps_dT || !dEps_dRho || !specific_neutrino_energy_loss || !specific_neutrino_flux || !mass_lost) { + return GF_UNINITIALIZED_OUTPUT_MEMORY_ERROR; + } + + if (tMax <= 0 || dt0 <= 0) { + return GF_INVALID_TIMESTEPS; + } + + if (num_species <= 0) { + return GF_INVALID_NUM_SPECIES; + } + + switch (type) { + case SINGLE_ZONE : { + auto* ctx = static_cast(ptr); + const auto T_ptr = static_cast(T); + const auto *rho_ptr = static_cast(rho); + + auto* Y_out_local = static_cast(Y_out); + auto* energy_out_local = static_cast(energy_out); + auto* dEps_dT_local = static_cast(dEps_dT); + auto* dEps_dRho_local = static_cast(dEps_dRho); + auto* specific_neutrino_energy_loss_local = static_cast(specific_neutrino_energy_loss); + auto* specific_neutrino_flux_local = static_cast(specific_neutrino_flux); + auto* mass_lost_local = static_cast(mass_lost); + + return execute_guarded(ctx, [&]() { + return ctx->evolve( + static_cast(Y_in), + num_species, + *T_ptr, + *rho_ptr, + tMax, + dt0, + Y_out_local, + *energy_out_local, + *dEps_dT_local, + *dEps_dRho_local, + *specific_neutrino_energy_loss_local, + *specific_neutrino_flux_local, + *mass_lost_local + ); + }); + } + case MULTI_ZONE : { + auto* ctx = static_cast(ptr); + const auto *T_ptr = static_cast(T); + const auto *rho_ptr = static_cast(rho); + + auto* Y_out_local = static_cast(Y_out); + auto* energy_out_local = static_cast(energy_out); + auto* dEps_dT_local = static_cast(dEps_dT); + auto* dEps_dRho_local = static_cast(dEps_dRho); + auto* specific_neutrino_energy_loss_local = static_cast(specific_neutrino_energy_loss); + auto* specific_neutrino_flux_local = static_cast(specific_neutrino_flux); + auto* mass_lost_local = static_cast(mass_lost); + + // for (size_t i = 0; i < ctx->get_zones(); ++i) { + // if (!Y_out_local[i]) { + // std::cerr << "Uninitialized memory for Y_out at zone " << i << std::endl; + // return GF_UNINITIALIZED_OUTPUT_MEMORY_ERROR; + // } + // } + + return execute_guarded(ctx, [&]() { + return ctx->evolve( + static_cast(Y_in), + num_species, + T_ptr, // T pointer + rho_ptr, // rho pointer + tMax, + dt0, + Y_out_local, + energy_out_local, + dEps_dT_local, + dEps_dRho_local, + specific_neutrino_energy_loss_local, + specific_neutrino_flux_local, + mass_lost_local + ); + }); + } + default : + return GF_UNKNOWN_ERROR; + } + + + + } + char* gf_get_last_error_message(void* ptr) { - const auto* ctx = static_cast(ptr); + if (!ptr) { + return const_cast("GF_UNINITIALIZED_INPUT_MEMORY_ERROR"); + } + const auto* ctx = static_cast(ptr); return const_cast(ctx->last_error.c_str()); } @@ -278,6 +386,18 @@ extern "C" { return const_cast("GF_DEBUG_ERROR"); case GF_GRIDFIRE_ERROR: return const_cast("GF_GRIDFIRE_ERROR"); + case GF_UNINITIALIZED_INPUT_MEMORY_ERROR: + return const_cast("GF_UNINITIALIZED_INPUT_MEMORY_ERROR"); + case GF_UNINITIALIZED_OUTPUT_MEMORY_ERROR: + return const_cast("GF_UNINITIALIZED_OUTPUT_MEMORY_ERROR"); + case GF_INVALID_NUM_SPECIES: + return const_cast("GF_INVALID_NUM_SPECIES"); + case GF_INVALID_TIMESTEPS: + return const_cast("GF_INVALID_TIMESTEPS"); + case GF_UNKNOWN_FREE_TYPE: + return const_cast("GF_UNKNOWN_FREE_TYPE"); + case GF_INVALID_TYPE: + return const_cast("GF_INVALID_TYPE"); case FDSSE_NON_4DSTAR_ERROR: return const_cast("FDSSE_NON_4DSTAR_ERROR"); case FDSSE_UNKNOWN_ERROR: diff --git a/tests/extern/C/gridfire_evolve_multi.c b/tests/extern/C/gridfire_evolve_multi.c new file mode 100644 index 00000000..6a414594 --- /dev/null +++ b/tests/extern/C/gridfire_evolve_multi.c @@ -0,0 +1,123 @@ +#include "gridfire/extern/gridfire_extern.h" +#include + +#define NUM_SPECIES 8 +#define ZONES 100 +#define ZONE_POLICY MULTI_ZONE + +// Define a macro to check return codes +#define CHECK_RET_CODE(ret, ctx, msg) \ + if ((ret) != 0 && (ret) != 1) { \ + printf("Error %s: %s (No. %d) [%s]\n", msg, gf_error_code_to_string(ret), ret, gf_get_last_error_message(ctx)); \ + gf_free(ZONE_POLICY, ctx); \ + return 1; \ + } + +int main() { + printf("Testing GridFireEvolve Multi\n"); + printf(" Number of zones: %d\n", ZONES); + printf(" Number of species: %d\n", NUM_SPECIES); + printf(" Initializing.."); + void* ctx = gf_init(ZONE_POLICY); + printf(" Done.\n"); + printf(" Setting number of zones to %d..", ZONES); + gf_set_num_zones(ZONE_POLICY, ctx, ZONES); + printf(" Done.\n"); + + const char* species_names[NUM_SPECIES]; + species_names[0] = "H-1"; + species_names[1] = "He-3"; + species_names[2] = "He-4"; + species_names[3] = "C-12"; + species_names[4] = "N-14"; + species_names[5] = "O-16"; + species_names[6] = "Ne-20"; + species_names[7] = "Mg-24"; + const double abundance_root[NUM_SPECIES] = { + 0.702616602672027, + 9.74791583949078e-06, + 0.06895512307276903, + 0.00025, + 7.855418029399437e-05, + 0.0006014411598306529, + 8.103062886768109e-05, + 2.151340851063217e-05 + }; + + double abundances[ZONES][NUM_SPECIES]; + for (size_t zone = 0; zone < ZONES; zone++) { + for (size_t i = 0; i < NUM_SPECIES; i++) { + abundances[zone][i] = abundance_root[i]; + } + } + + double Temps[ZONES]; + double Rhos[ZONES]; + + for (size_t zone = 0; zone < ZONES; zone++) { + Temps[zone] = 1.0e7; + Rhos[zone] = 1.5e2; + } + + printf(" Registering species..."); + int ret = gf_register_species(ctx, NUM_SPECIES, species_names); + CHECK_RET_CODE(ret, ctx, "SPECIES"); + printf(" Done.\n"); + + printf(" Constructing engine from policy..."); + ret = gf_construct_engine_from_policy(ctx, "MAIN_SEQUENCE_POLICY", abundance_root, NUM_SPECIES); + CHECK_RET_CODE(ret, ctx, "MAIN_SEQUENCE_POLICY"); + printf(" Done.\n"); + + printf(" Constructing solver from engine..."); + ret = gf_construct_solver_from_engine(ctx); + CHECK_RET_CODE(ret, ctx, "CVODE"); + printf(" Done.\n"); + + double Y_out[ZONES][NUM_SPECIES]; + double energy_out[ZONES]; + double dEps_dT[ZONES]; + double dEps_dRho[ZONES]; + double specific_neutrino_energy_loss[ZONES]; + double specific_neutrino_flux[ZONES]; + double mass_lost[ZONES]; + + printf(" Evolving...\n"); + ret = gf_evolve( + ZONE_POLICY, + ctx, + abundances, + NUM_SPECIES, + Temps, // Temperature in K + Rhos, // Density in g/cm^3 + 3e17, // Time step in seconds + 1e-12, + Y_out, + energy_out, + dEps_dT, + dEps_dRho, + specific_neutrino_energy_loss, + specific_neutrino_flux, &mass_lost + ); + CHECK_RET_CODE(ret, ctx, "EVOLUTION"); + printf(" Done.\n"); + + + printf("Evolved abundances:\n"); + for (size_t zone = 0; zone < ZONES; zone++) { + printf("=== Zone %zu ===\n", zone); + for (size_t i = 0; i < NUM_SPECIES; i++) { + printf(" Species %s: %e\n", species_names[i], Y_out[zone][i]); + } + printf(" Energy output: %e\n", energy_out[zone]); + printf(" dEps/dT: %e\n", dEps_dT[zone]); + printf(" dEps/dRho: %e\n", dEps_dRho[zone]); + printf(" Specific neutrino energy loss: %e\n", specific_neutrino_energy_loss[zone]); + printf(" Specific neutrino flux: %e\n", specific_neutrino_flux[zone]); + printf(" Mass lost: %e\n", mass_lost[zone]); + } + + gf_free(ZONE_POLICY, ctx); + + return 0; +} \ No newline at end of file diff --git a/tests/extern/C/gridfire_evolve.c b/tests/extern/C/gridfire_evolve_single.c similarity index 78% rename from tests/extern/C/gridfire_evolve.c rename to tests/extern/C/gridfire_evolve_single.c index 7834dabf..f7efb9a0 100644 --- a/tests/extern/C/gridfire_evolve.c +++ b/tests/extern/C/gridfire_evolve_single.c @@ -4,14 +4,14 @@ // Define a macro to check return codes #define CHECK_RET_CODE(ret, ctx, msg) \ - if (ret != 0 && ret != 1) { \ - printf("Error %s: %s\n", msg, gf_get_last_error_message(ctx)); \ - gf_free(ctx); \ + if ((ret) != 0 && (ret) != 1) { \ + printf("Error %s: %s [%s]\n", msg, gf_get_last_error_message(ctx), gf_error_code_to_string(ret)); \ + gf_free(SINGLE_ZONE, ctx); \ return 1; \ } int main() { - void* ctx = gf_init(); + void* ctx = gf_init(SINGLE_ZONE); const char* species_names[NUM_SPECIES]; species_names[0] = "H-1"; @@ -39,7 +39,7 @@ int main() { ret = gf_construct_engine_from_policy(ctx, "MAIN_SEQUENCE_POLICY", abundances, NUM_SPECIES); CHECK_RET_CODE(ret, ctx, "MAIN_SEQUENCE_POLICY"); - ret = gf_construct_solver_from_engine(ctx, "CVODE"); + ret = gf_construct_solver_from_engine(ctx); CHECK_RET_CODE(ret, ctx, "CVODE"); double Y_out[NUM_SPECIES]; @@ -50,21 +50,24 @@ int main() { double specific_neutrino_flux; double mass_lost; + const double T_in = 1.5e7; // Temperature in K + const double rho_in = 1.5e2; // Density in g/cm^3 + ret = gf_evolve( + SINGLE_ZONE, ctx, abundances, NUM_SPECIES, - 1.5e7, // Temperature in K - 1.5e2, // Density in g/cm^3 + &T_in, // Temperature in K + &rho_in, // Density in g/cm^3 3e17, // Time step in seconds - 1e-12, // Initial time step in seconds + 1e-12, Y_out, &energy_out, &dEps_dT, &dEps_dRho, &specific_neutrino_energy_loss, - &specific_neutrino_flux, - &mass_lost + &specific_neutrino_flux, &mass_lost ); CHECK_RET_CODE(ret, ctx, "EVOLUTION"); @@ -81,7 +84,7 @@ int main() { printf("Specific neutrino flux: %e\n", specific_neutrino_flux); printf("Mass lost: %e\n", mass_lost); - gf_free(ctx); + gf_free(SINGLE_ZONE, ctx); return 0; } \ No newline at end of file diff --git a/tests/extern/C/meson.build b/tests/extern/C/meson.build index 80db9f58..98f44f5e 100644 --- a/tests/extern/C/meson.build +++ b/tests/extern/C/meson.build @@ -1,5 +1,11 @@ -executable('test_c_extern', 'gridfire_evolve.c', +executable('gf_extern_c_single', 'gridfire_evolve_single.c', install: false, c_args: ['-Wall', '-Wextra'], dependencies: [gridfire_extern_dep] -) \ No newline at end of file +) + +executable('gf_extern_c_multi', 'gridfire_evolve_multi.c', + install: false, + c_args: ['-Wall', '-Wextra'], + dependencies: [gridfire_extern_dep] +)