tboudreaux/GridFire
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #8 from tboudreaux/feature/extern

Browse Source 
Stable C API and Fortran Bindings
This commit is contained in:
Emily Boudreaux
2025-11-27 11:33:48 -05:00
committed by GitHub
parent a46bd846a5 737f432b94
commit 6adff5b670
28 changed files with 1976 additions and 24809 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Doxyfile
View File

@@ -48,7 +48,7 @@ PROJECT_NAME = GridFire
# could be handy for archiving the generated documentation or if some version
# control system is used.
PROJECT_NUMBER = v0.7.0_rc1
PROJECT_NUMBER = v0.7.0_rc2
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewers a

218
README.md
View File

@@ -9,7 +9,7 @@
![GitHub License](https://img.shields.io/github/license/4D-STAR/GridFire?style=for-the-badge)
![ERC](https://img.shields.io/badge/Funded%20by-ERC-blue?style=for-the-badge&logo=europeancommission)
![Dynamic Regex Badge](https://img.shields.io/badge/dynamic/regex?url=https%3A%2F%2Fgithub.com%2F4D-STAR%2FGridFire%2Fblob%2Fmain%2Fmeson.build&search=version%3A%20'(%5B0-9a-zA-Z%5C.%5D%2B)'&style=for-the-badge&label=GitHub%20Main%20Branch)
![Dynamic Regex Badge](https://img.shields.io/badge/dynamic/regex?url=https%3A%2F%2Fraw.githubusercontent.com%2F4D-STAR%2FGridFire%2Frefs%2Fheads%2Fmain%2Fmeson.build&search=version%3A%20'(.%2B)'%2C&style=for-the-badge&label=Main%20Branch%20Version)
![GitHub commit activity](https://img.shields.io/github/commit-activity/w/4D-STAR/GridFire?style=for-the-badge)
@@ -896,6 +896,220 @@ if __name__ == "__main__":
```
# External Usage
C++ does not have a stable ABI nor does it make any strong guarantees about stl container layouts between compiler versions.
Therefore, GridFire includes a set of stable C bindings which can be used to interface with a limited subset of GridFire functionality
from other languages.
> **Note:** These bindings are not intended to allow GridFire to be extended from other languages; rather, they are intended to allow GridFire to be used as a
> black-box library from other languages.
> **Note:** One assumption for external usage is that the ordering of the species list will not change. That is to say that whatever order the array
> used to register the species is will be assumed to always be the order used when passing abundance arrays to and from GridFire.
> **Note:** Because the C API does not pass the general Composition object a `mass_lost`
> output parameter has been added to the evolve calls, this tracks the total mass in species which have not been registered with the C API GridFire by the caller
## C API Overview
In general when using the C API the workflow is to
1. create a `gf_context` pointer. This object holds the state of GridFire so that it does not need to be re-initialized for each call.
2. call initialization routines on the context to set up the engine and solver you wish to use.
3. call the `gf_evolve` function to evolve a network over some time.
4. At each state check the ret code of the function to ensure that no errors occurred. Valid ret-codes are 0 and 1. All other ret codes indicate an error.
5. Finally, call `gf_free` to free the context and all associated memory.
### C Example
```c++
#include "gridfire/extern/gridfire_extern.h"
#include <stdio.h>
#define NUM_SPECIES 8
// Define a macro to check return codes
#define GF_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); \
return ret; \
}
int main() {
void* gf_context = gf_init();
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 abundances[NUM_SPECIES] = {0.702616602672027, 9.74791583949078e-06, 0.06895512307276903, 0.00025, 7.855418029399437e-05, 0.0006014411598306529, 8.103062886768109e-05, 2.151340851063217e-05};
int ret = gf_register_species(gf_context, NUM_SPECIES, species_names);
GF_CHECK_RET_CODE(ret, gf_context, "Species Registration");
ret = gf_construct_engine_from_policy(gf_context, "MAIN_SEQUENCE_POLICY", abundances, NUM_SPECIES);
GF_CHECK_RET_CODE(ret, gf_context, "Policy and Engine Construction");
ret = gf_construct_solver_from_engine(gf_context, "CVODE");
GF_CHECK_RET_CODE(ret, gf_context, "Solver Construction");
// When using the C API it is assumed that the caller will ensure that the output arrays are large enough to hold the results.
double Y_out[NUM_SPECIES];
double energy_out;
double dEps_dT;
double dEps_dRho;
double mass_lost;
ret = gf_evolve(
gf_context,
abundances,
NUM_SPECIES,
1.5e7, // Temperature in K
1.5e2, // Density in g/cm^3
3e17, // Time step in seconds
1e-12, // Initial time step in seconds
Y_out,
&energy_out,
&dEps_dT,
&dEps_dRho, &mass_lost
);
GF_CHECK_RET_CODE(ret, gf_context, "Evolution");
printf("Evolved abundances:\n");
for (size_t i = 0; i < NUM_SPECIES; i++) {
printf("Species %s: %e\n", species_names[i], Y_out[i]);
}
printf("Energy output: %e\n", energy_out);
printf("dEps/dT: %e\n", dEps_dT);
printf("dEps/dRho: %e\n", dEps_dRho);
printf("Mass lost: %e\n", mass_lost);
gf_free(gf_context);
return 0;
}
```
## Fortran API Overview
GridFire makes use of the stable C API and Fortran 2003's `iso_c_bindings` to provide a Fortran interface for legacy
code. The fortran interface is designed to be very similar to the C API and exposes the same functionality.
1. `GridFire%gff_init`: Initializes a GridFire context and returns a handle to it.
2. `GridFire%register_species`: Registers species with the GridFire context.
3. `GridFire%setup_policy`: Configures the engine using a specified policy and initial abundances.
4. `GridFire%setup_solver`: Sets up the solver for the engine.
5. `GridFire%evolve`: Evolves the network over a specified time step.
6. `GridFire%get_last_error`: Retrieves the last error message from the GridFire context.
7. `GridFire%gff_free`: Frees the GridFire context and associated resources.
> **Note:** You must instantiate a `GridFire` type object to access these methods.
> **Note:** free and init have had the `gff_` prefix (GridFire Fortran) to avoid name clashes with common Fortran functions.
When building GridFire a fortran module file `gridfire_mod.mod` is generated which contains all the necessary
bindings to use GridFire from Fortran. You must also link your code against the C API library `libgridfire_extern`.
### Fortran Example
```fortran
program main
use iso_c_binding
use gridfire_mod
implicit none
type(GridFire) :: net
integer(c_int) :: ierr
integer :: i
! --- 1. Define Species and Initial Conditions ---
! Note: String lengths must match or exceed the longest name.
! We pad with spaces, which 'trim' handles inside the module.
character(len=5), dimension(8) :: species_names = [ &
"H-1 ", &
"He-3 ", &
"He-4 ", &
"C-12 ", &
"N-14 ", &
"O-16 ", &
"Ne-20", &
"Mg-24" &
]
! Initial Mass Fractions (converted to Molar Abundances Y = X/A)
! Standard solar-ish composition
real(c_double), dimension(8) :: Y_in = [ &
0.702616602672027, &
9.74791583949078e-06, &
0.06895512307276903, &
0.00025, &
7.855418029399437e-05, &
0.0006014411598306529, &
8.103062886768109e-05, &
2.151340851063217e-05 &
]
! Output buffers
real(c_double), dimension(8) :: Y_out
real(c_double) :: energy_out, dedt, dedrho, dmass
! Thermodynamic Conditions (Solar Core-ish)
real(c_double) :: T = 1.5e7 ! 15 Million K
real(c_double) :: rho = 150.0e0 ! 150 g/cm^3
real(c_double) :: dt = 3.1536e17 ! ~10 Gyr timestep
! --- 2. Initialize GridFire ---
print *, "Initializing GridFire..."
call net%gff_init()
! --- 3. Register Species ---
print *, "Registering species..."
call net%register_species(species_names)
! --- 4. Configure Engine & Solver ---
print *, "Setting up Main Sequence Policy..."
call net%setup_policy("MAIN_SEQUENCE_POLICY", Y_in)
print *, "Setting up CVODE Solver..."
call net%setup_solver("CVODE")
! --- 5. Evolve ---
print *, "Evolving system (dt =", dt, "s)..."
call net%evolve(Y_in, T, rho, dt, Y_out, energy_out, dedt, dedrho, dmass, ierr)
if (ierr /= 0) then
print *, "Evolution Failed with error code: ", ierr
print *, "Error Message: ", net%get_last_error()
call net%gff_free() ! Always cleanup
stop
end if
! --- 6. Report Results ---
print *, ""
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)', "Mass Change: ", dmass
print *, ""
print *, "Abundances:"
do i = 1, size(species_names)
print '(A, " : ", ES12.5, " -> ", ES12.5)', &
trim(species_names(i)), Y_in(i), Y_out(i)
end do
! --- 7. Cleanup ---
call net%gff_free()
end program main
```
# Related Projects
@@ -910,3 +1124,5 @@ GridFire integrates with and builds upon several key 4D-STAR libraries:
utilities.
- [liblogging](https://github.com/4D-STAR/liblogging): Flexible logging framework.
- [libconstants](https://github.com/4D-STAR/libconstants): Physical constants
- [libplugin](https://github.com/4D-STAR/libplugin): Dynamically loadable plugin
framework.

218
docs/static/mainpage.md vendored
View File

@@ -9,7 +9,7 @@
![GitHub License](https://img.shields.io/github/license/4D-STAR/GridFire?style=for-the-badge)
![ERC](https://img.shields.io/badge/Funded%20by-ERC-blue?style=for-the-badge&logo=europeancommission)
![Dynamic Regex Badge](https://img.shields.io/badge/dynamic/regex?url=https%3A%2F%2Fgithub.com%2F4D-STAR%2FGridFire%2Fblob%2Fmain%2Fmeson.build&search=version%3A%20'(%5B0-9a-zA-Z%5C.%5D%2B)'&style=for-the-badge&label=GitHub%20Main%20Branch)
![Dynamic Regex Badge](https://img.shields.io/badge/dynamic/regex?url=https%3A%2F%2Fraw.githubusercontent.com%2F4D-STAR%2FGridFire%2Frefs%2Fheads%2Fmain%2Fmeson.build&search=version%3A%20'(.%2B)'%2C&style=for-the-badge&label=Main%20Branch%20Version)
![GitHub commit activity](https://img.shields.io/github/commit-activity/w/4D-STAR/GridFire?style=for-the-badge)
@@ -896,6 +896,220 @@ if __name__ == "__main__":
```
# External Usage
C++ does not have a stable ABI nor does it make any strong guarantees about stl container layouts between compiler versions.
Therefore, GridFire includes a set of stable C bindings which can be used to interface with a limited subset of GridFire functionality
from other languages.
> **Note:** These bindings are not intended to allow GridFire to be extended from other languages; rather, they are intended to allow GridFire to be used as a
> black-box library from other languages.
> **Note:** One assumption for external usage is that the ordering of the species list will not change. That is to say that whatever order the array
> used to register the species is will be assumed to always be the order used when passing abundance arrays to and from GridFire.
> **Note:** Because the C API does not pass the general Composition object a `mass_lost`
> output parameter has been added to the evolve calls, this tracks the total mass in species which have not been registered with the C API GridFire by the caller
## C API Overview
In general when using the C API the workflow is to
1. create a `gf_context` pointer. This object holds the state of GridFire so that it does not need to be re-initialized for each call.
2. call initialization routines on the context to set up the engine and solver you wish to use.
3. call the `gf_evolve` function to evolve a network over some time.
4. At each state check the ret code of the function to ensure that no errors occurred. Valid ret-codes are 0 and 1. All other ret codes indicate an error.
5. Finally, call `gf_free` to free the context and all associated memory.
### C Example
```c++
#include "gridfire/extern/gridfire_extern.h"
#include <stdio.h>
#define NUM_SPECIES 8
// Define a macro to check return codes
#define GF_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); \
return ret; \
}
int main() {
void* gf_context = gf_init();
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 abundances[NUM_SPECIES] = {0.702616602672027, 9.74791583949078e-06, 0.06895512307276903, 0.00025, 7.855418029399437e-05, 0.0006014411598306529, 8.103062886768109e-05, 2.151340851063217e-05};
int ret = gf_register_species(gf_context, NUM_SPECIES, species_names);
GF_CHECK_RET_CODE(ret, gf_context, "Species Registration");
ret = gf_construct_engine_from_policy(gf_context, "MAIN_SEQUENCE_POLICY", abundances, NUM_SPECIES);
GF_CHECK_RET_CODE(ret, gf_context, "Policy and Engine Construction");
ret = gf_construct_solver_from_engine(gf_context, "CVODE");
GF_CHECK_RET_CODE(ret, gf_context, "Solver Construction");
// When using the C API it is assumed that the caller will ensure that the output arrays are large enough to hold the results.
double Y_out[NUM_SPECIES];
double energy_out;
double dEps_dT;
double dEps_dRho;
double mass_lost;
ret = gf_evolve(
gf_context,
abundances,
NUM_SPECIES,
1.5e7, // Temperature in K
1.5e2, // Density in g/cm^3
3e17, // Time step in seconds
1e-12, // Initial time step in seconds
Y_out,
&energy_out,
&dEps_dT,
&dEps_dRho, &mass_lost
);
GF_CHECK_RET_CODE(ret, gf_context, "Evolution");
printf("Evolved abundances:\n");
for (size_t i = 0; i < NUM_SPECIES; i++) {
printf("Species %s: %e\n", species_names[i], Y_out[i]);
}
printf("Energy output: %e\n", energy_out);
printf("dEps/dT: %e\n", dEps_dT);
printf("dEps/dRho: %e\n", dEps_dRho);
printf("Mass lost: %e\n", mass_lost);
gf_free(gf_context);
return 0;
}
```
## Fortran API Overview
GridFire makes use of the stable C API and Fortran 2003's `iso_c_bindings` to provide a Fortran interface for legacy
code. The fortran interface is designed to be very similar to the C API and exposes the same functionality.
1. `GridFire%gff_init`: Initializes a GridFire context and returns a handle to it.
2. `GridFire%register_species`: Registers species with the GridFire context.
3. `GridFire%setup_policy`: Configures the engine using a specified policy and initial abundances.
4. `GridFire%setup_solver`: Sets up the solver for the engine.
5. `GridFire%evolve`: Evolves the network over a specified time step.
6. `GridFire%get_last_error`: Retrieves the last error message from the GridFire context.
7. `GridFire%gff_free`: Frees the GridFire context and associated resources.
> **Note:** You must instantiate a `GridFire` type object to access these methods.
> **Note:** free and init have had the `gff_` prefix (GridFire Fortran) to avoid name clashes with common Fortran functions.
When building GridFire a fortran module file `gridfire_mod.mod` is generated which contains all the necessary
bindings to use GridFire from Fortran. You must also link your code against the C API library `libgridfire_extern`.
### Fortran Example
```fortran
program main
use iso_c_binding
use gridfire_mod
implicit none
type(GridFire) :: net
integer(c_int) :: ierr
integer :: i
! --- 1. Define Species and Initial Conditions ---
! Note: String lengths must match or exceed the longest name.
! We pad with spaces, which 'trim' handles inside the module.
character(len=5), dimension(8) :: species_names = [ &
"H-1 ", &
"He-3 ", &
"He-4 ", &
"C-12 ", &
"N-14 ", &
"O-16 ", &
"Ne-20", &
"Mg-24" &
]
! Initial Mass Fractions (converted to Molar Abundances Y = X/A)
! Standard solar-ish composition
real(c_double), dimension(8) :: Y_in = [ &
0.702616602672027, &
9.74791583949078e-06, &
0.06895512307276903, &
0.00025, &
7.855418029399437e-05, &
0.0006014411598306529, &
8.103062886768109e-05, &
2.151340851063217e-05 &
]
! Output buffers
real(c_double), dimension(8) :: Y_out
real(c_double) :: energy_out, dedt, dedrho, dmass
! Thermodynamic Conditions (Solar Core-ish)
real(c_double) :: T = 1.5e7 ! 15 Million K
real(c_double) :: rho = 150.0e0 ! 150 g/cm^3
real(c_double) :: dt = 3.1536e17 ! ~10 Gyr timestep
! --- 2. Initialize GridFire ---
print *, "Initializing GridFire..."
call net%gff_init()
! --- 3. Register Species ---
print *, "Registering species..."
call net%register_species(species_names)
! --- 4. Configure Engine & Solver ---
print *, "Setting up Main Sequence Policy..."
call net%setup_policy("MAIN_SEQUENCE_POLICY", Y_in)
print *, "Setting up CVODE Solver..."
call net%setup_solver("CVODE")
! --- 5. Evolve ---
print *, "Evolving system (dt =", dt, "s)..."
call net%evolve(Y_in, T, rho, dt, Y_out, energy_out, dedt, dedrho, dmass, ierr)
if (ierr /= 0) then
print *, "Evolution Failed with error code: ", ierr
print *, "Error Message: ", net%get_last_error()
call net%gff_free() ! Always cleanup
stop
end if
! --- 6. Report Results ---
print *, ""
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)', "Mass Change: ", dmass
print *, ""
print *, "Abundances:"
do i = 1, size(species_names)
print '(A, " : ", ES12.5, " -> ", ES12.5)', &
trim(species_names(i)), Y_in(i), Y_out(i)
end do
! --- 7. Cleanup ---
call net%gff_free()
end program main
```
# Related Projects
@@ -910,3 +1124,5 @@ GridFire integrates with and builds upon several key 4D-STAR libraries:
utilities.
- [liblogging](https://github.com/4D-STAR/liblogging): Flexible logging framework.
- [libconstants](https://github.com/4D-STAR/libconstants): Physical constants
- [libplugin](https://github.com/4D-STAR/libplugin): Dynamically loadable plugin
framework.

6
meson.build
View File

@@ -18,12 +18,16 @@
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#
# *********************************************************************** #
project('GridFire', 'cpp', version: 'v0.7.0_rc1', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
project('GridFire', ['c', 'cpp', 'fortran'], version: 'v0.7.0_rc2', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
# Add default visibility for all C++ targets
add_project_arguments('-fvisibility=default', language: 'cpp')
message('Found CXX compiler: ' + meson.get_compiler('cpp').get_id())
message('Found FORTRAN compiler: ' + meson.get_compiler('fortran').get_id())
message('C++ standard set to: ' + get_option('cpp_std'))
message('Fortran standard set to: ' + get_option('fortran_std'))
if meson.get_compiler('cpp').get_id() == 'clang'
# We disable these because of CppAD

258
src/extern/fortran/gridfire_mod.f90 vendored Normal file
View File

@@ -0,0 +1,258 @@
module gridfire_mod
use iso_c_binding
implicit none
enum, bind (C)
enumerator :: FDSSE_NON_4DSTAR_ERROR = -102
enumerator :: FDSSE_UNKNOWN_ERROR = -101
enumerator :: FDSSE_SUCCESS = 1
enumerator :: FDSSE_UNKNOWN_SYMBOL_ERROR = 100
enumerator :: FDSSE_SPECIES_ERROR = 101
enumerator :: FDSSE_INVALID_COMPOSITION_ERROR = 102
enumerator :: FDSSE_COMPOSITION_ERROR = 103
enumerator :: GF_NON_GRIDFIRE_ERROR = -2
enumerator :: GF_UNKNOWN_ERROR = -1
enumerator :: GF_SUCCESS = 0
enumerator :: GF_INVALID_QSE_SOLUTION_ERROR = 5
enumerator :: GF_FAILED_TO_PARTITION_ERROR = 6
enumerator :: GF_NETWORK_RESIZED_ERROR = 7
enumerator :: GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR = 8
enumerator :: GF_BAD_COLLECTION_ERROR = 9
enumerator :: GF_BAD_RHS_ENIGNE_ERROR = 10
enumerator :: GF_STALE_JACOBIAN_ERROR = 11
enumerator :: GF_UNINITIALIZED_JACOBIAN_ERROR = 12
enumerator :: GF_UNKNONWN_JACOBIAN_ERROR = 13
enumerator :: GF_JACOBIAN_ERROR = 14
enumerator :: GF_ENGINE_ERROR = 15
enumerator :: GF_MISSING_BASE_REACTION_ERROR = 16
enumerator :: GF_MISSING_SEED_SPECIES_ERROR = 17
enumerator :: GF_MISSING_KEY_REACTION_ERROR = 18
enumerator :: GF_POLICY_ERROR = 19
enumerator :: GF_REACTION_PARSING_ERROR = 20
enumerator :: GF_REACTOION_ERROR = 21
enumerator :: GF_SINGULAR_JACOBIAN_ERROR = 22
enumerator :: GF_ILL_CONDITIONED_JACOBIAN_ERROR = 23
enumerator :: GF_CVODE_SOLVER_FAILURE_ERROR = 24
enumerator :: GF_KINSOL_SOLVER_FAILURE_ERROR = 25
enumerator :: GF_SUNDIALS_ERROR = 26
enumerator :: GF_SOLVER_ERROR = 27
enumerator :: GF_HASHING_ERROR = 28
enumerator :: GF_UTILITY_ERROR = 29
enumerator :: GF_DEBUG_ERRROR = 30
enumerator :: GF_GRIDFIRE_ERROR = 31
end enum
interface
! void* gf_init()
function gf_init() bind(C, name="gf_init")
import :: c_ptr
type(c_ptr) :: gf_init
end function gf_init
! void gf_free(void* gf)
subroutine gf_free(gf) bind(C, name="gf_free")
import :: c_ptr
type(c_ptr), value :: gf
end subroutine gf_free
! char* gf_get_last_error_message(void* ptr);
function gf_get_last_error_message(ptr) result(c_msg) bind(C, name="gf_get_last_error_message")
import
type(c_ptr), value :: ptr
type(c_ptr) :: c_msg
end function
! int gf_register_species(void* ptr, const int num_species, const char** species_names);
function gf_register_species(ptr, num_species, species_names) result(ierr) bind(C, name="gf_register_species")
import
type(c_ptr), value :: ptr
integer(c_int), value :: num_species
type(c_ptr), dimension(*), intent(in) :: species_names ! Array of C pointers
integer(c_int) :: ierr
end function
! int gf_construct_engine_from_policy(void* ptr, const char* policy_name, const double *abundances, size_t num_species);
function gf_construct_engine_from_policy(ptr, policy_name, abundances, num_species) result(ierr) &
bind(C, name="gf_construct_engine_from_policy")
import
type(c_ptr), value :: ptr
character(kind=c_char), dimension(*), intent(in) :: policy_name
real(c_double), dimension(*), intent(in) :: abundances
integer(c_size_t), value :: num_species
integer(c_int) :: ierr
end function
! int gf_construct_solver_from_engine(void* ptr, const char* solver_name);
function gf_construct_solver_from_engine(ptr, solver_name) result(ierr) &
bind(C, name="gf_construct_solver_from_engine")
import
type(c_ptr), value :: ptr
character(kind=c_char), dimension(*), intent(in) :: solver_name
integer(c_int) :: ierr
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) &
bind(C, name="gf_evolve")
import
type(c_ptr), value :: ptr
real(c_double), dimension(*), intent(in) :: Y_in
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
integer(c_int) :: ierr
end function
end interface
type :: GridFire
type(c_ptr) :: ctx = c_null_ptr
integer(c_size_t) :: num_species = 0
contains
procedure :: gff_init
procedure :: gff_free
procedure :: register_species
procedure :: setup_policy
procedure :: setup_solver
procedure :: evolve
procedure :: get_last_error
end type GridFire
contains
subroutine gff_init(self)
class(GridFire), intent(out) :: self
self%ctx = gf_init()
end subroutine gff_init
subroutine gff_free(self)
class(GridFire), intent(inout) :: self
if (c_associated(self%ctx)) then
call gf_free(self%ctx)
self%ctx = c_null_ptr
end if
end subroutine gff_free
function get_last_error(self) result(msg)
class(GridFire), intent(in) :: self
character(len=:), allocatable :: msg
type(c_ptr) :: c_msg_ptr
character(kind=c_char), pointer :: char_ptr(:)
integer :: i, len_str
c_msg_ptr = gf_get_last_error_message(self%ctx)
if (.not. c_associated(c_msg_ptr)) then
msg = "GridFire: Unknown Error (Null Pointer returned)"
return
end if
call c_f_pointer(c_msg_ptr, char_ptr, [1024])
len_str = 0
do i = 1, 1024
if (char_ptr(i) == c_null_char) exit
len_str = len_str + 1
end do
msg = repeat(' ', len_str+10)
msg(1:10) = "GridFire: "
do i = 1, len_str
msg(i+10:i+10) = char_ptr(i)
end do
end function get_last_error
subroutine register_species(self, species_list)
class(GridFire), intent(inout) :: self
character(len=*), dimension(:), intent(in) :: species_list
type(c_ptr), allocatable, dimension(:) :: c_ptrs
character(kind=c_char, len=:), allocatable, target :: temp_strs(:)
integer :: i, n, ierr
print *, "Registering ", size(species_list), " species."
n = size(species_list)
self%num_species = int(n, c_size_t)
allocate(c_ptrs(n))
allocate(character(len=len(species_list(1))+1) :: temp_strs(n)) ! +1 for null terminator
do i = 1, n
temp_strs(i) = trim(species_list(i)) // c_null_char
c_ptrs(i) = c_loc(temp_strs(i))
end do
print *, "Calling gf_register_species..."
ierr = gf_register_species(self%ctx, int(n, c_int), c_ptrs)
print *, "gf_register_species returned with code: ", ierr
if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then
print *, "GridFire: ", self%get_last_error()
error stop
end if
end subroutine register_species
subroutine setup_policy(self, policy_name, abundances)
class(GridFire), intent(in) :: self
character(len=*), intent(in) :: policy_name
real(c_double), dimension(:), intent(in) :: abundances
integer(c_int) :: ierr
if (size(abundances) /= self%num_species) then
print *, "GridFire Error: Abundance array size mismatch."
error stop
end if
ierr = gf_construct_engine_from_policy(self%ctx, &
trim(policy_name) // c_null_char, &
abundances, &
self%num_species)
if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then
print *, "GridFire Policy Error: ", self%get_last_error()
error stop
end if
end subroutine setup_policy
subroutine setup_solver(self, solver_name)
class(GridFire), intent(in) :: self
character(len=*), intent(in) :: solver_name
integer(c_int) :: ierr
ierr = gf_construct_solver_from_engine(self%ctx, trim(solver_name) // c_null_char)
if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then
print *, "GridFire Solver Error: ", self%get_last_error()
error stop
end if
end subroutine setup_solver
subroutine evolve(self, Y_in, T, rho, dt, Y_out, energy, dedt, dedrho, 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
integer, intent(out) :: ierr
integer(c_int) :: c_ierr
c_ierr = gf_evolve(self%ctx, &
Y_in, self%num_species, &
T, rho, dt, &
Y_out, &
energy, dedt, dedrho, mass_lost)
ierr = int(c_ierr)
if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then
print *, "GridFire Evolve Error: ", self%get_last_error()
end if
end subroutine evolve
end module gridfire_mod

11
src/extern/fortran/meson.build vendored Normal file
View File

@@ -0,0 +1,11 @@
gridfire_fortran_lib = library('gridfire_fortran',
'gridfire_mod.f90',
link_with: libgridfire_extern,
install: true,
install_dir: get_option('libdir')
)
gridfire_fortran_dep = declare_dependency(
link_with: [gridfire_fortran_lib, libgridfire_extern],
include_directories: include_directories('.')
)

40
src/extern/include/gridfire/extern/gridfire_context.h vendored Normal file
View File

@@ -0,0 +1,40 @@
#ifndef GF_GRIDFIRE_CONTEXT_H
#define GF_GRIDFIRE_CONTEXT_H
#include "gridfire/gridfire.h"
#include "fourdst/atomic/atomicSpecies.h"
#include <memory>
#include <vector>
struct GridFireContext {
std::unique_ptr<gridfire::policy::NetworkPolicy> policy;
gridfire::engine::DynamicEngine* engine;
std::unique_ptr<gridfire::solver::DynamicNetworkSolverStrategy> solver;
std::vector<fourdst::atomic::Species> speciesList;
fourdst::composition::Composition working_comp;
void init_species_map(const std::vector<std::string>& 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);
void init_composition_from_abundance_vector(const double* abundances, size_t num_species);
int evolve(
const double* Y_in,
size_t num_species,
double T,
double rho,
double tMax,
double dt0,
double* Y_out,
double& energy_out,
double& dEps_dT,
double& dEps_dRho, double& mass_lost
);
std::string last_error;
};
#endif

91
src/extern/include/gridfire/extern/gridfire_extern.h vendored Normal file
View File

@@ -0,0 +1,91 @@
#ifndef GF_GRIDFIRE_EXTERN_H
#define GF_GRIDFIRE_EXTERN_H
#include <stddef.h>
#ifdef __cplusplus
extern "C" {
#endif
enum FDSSE_ERROR_CODES {
FDSSE_NON_4DSTAR_ERROR = -102,
FDSSE_UNKNOWN_ERROR = -101,
FDSSE_SUCCESS = 1,
FDSSE_UNKNOWN_SYMBOL_ERROR = 100,
FDSSE_SPECIES_ERROR = 101,
FDSSE_INVALID_COMPOSITION_ERROR = 102,
FDSSE_COMPOSITION_ERROR = 103
};
enum GF_ERROR_CODES {
GF_NON_GRIDFIRE_ERROR = -2,
GF_UNKNOWN_ERROR = -1,
GF_SUCCESS = 0,
GF_INVALID_QSE_SOLUTION_ERROR = 5,
GF_FAILED_TO_PARTITION_ENGINE_ERROR = 6,
GF_NETWORK_RESIZED_ERROR = 7,
GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR = 8,
GF_BAD_COLLECTION_ERROR = 9,
GF_BAD_RHS_ENGINE_ERROR = 10,
GF_STALE_JACOBIAN_ERROR = 11,
GF_UNINITIALIZED_JACOBIAN_ERROR = 12,
GF_UNKNOWN_JACOBIAN_ERROR = 13,
GF_JACOBIAN_ERROR = 14,
GF_ENGINE_ERROR = 15,
GF_MISSING_BASE_REACTION_ERROR = 16,
GF_MISSING_SEED_SPECIES_ERROR = 17,
GF_MISSING_KEY_REACTION_ERROR = 18,
GF_POLICY_ERROR = 19,
GF_REACTION_PARSING_ERROR = 20,
GF_REACTION_ERROR = 21,
GF_SINGULAR_JACOBIAN_ERROR = 22,
GF_ILL_CONDITIONED_JACOBIAN_ERROR = 23,
GF_CVODE_SOLVER_FAILURE_ERROR = 24,
GF_KINSOL_SOLVER_FAILURE_ERROR = 25,
GF_SUNDIALS_ERROR = 26,
GF_SOLVER_ERROR = 27,
GF_HASHING_ERROR = 28,
GF_UTILITY_ERROR = 29,
GF_DEBUG_ERROR = 30,
GF_GRIDFIRE_ERROR = 31,
};
char* gf_get_last_error_message(void* ptr);
char* gf_error_code_to_string(int error_code);
void* gf_init();
void gf_free(void* ctx);
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_evolve(
void* ptr,
const double* Y_in,
size_t num_species,
double T,
double rho,
double tMax,
double dt0,
double* Y_out,
double* energy_out,
double* dEps_dT,
double* dEps_dRho, double* mass_lost
);
#ifdef __cplusplus
}
#endif
#endif

156
src/extern/lib/gridfire_context.cpp vendored Normal file
View File

@@ -0,0 +1,156 @@
#include "gridfire/extern/gridfire_context.h"
#include "fourdst/atomic/species.h"
#include "fourdst/composition/exceptions/exceptions_composition.h"
void GridFireContext::init_species_map(const std::vector<std::string> &species_names) {
for (const auto& name: species_names) {
working_comp.registerSymbol(name);
}
this->speciesList.clear();
this->speciesList.reserve(species_names.size());
auto resolve_species_name = [](const std::string& name) -> fourdst::atomic::Species {
if (fourdst::atomic::species.contains(name)) {
return fourdst::atomic::species.at(name);
}
throw fourdst::composition::exceptions::UnknownSymbolError("Species " + name + " is not recognized in the atomic species database.");
};
for (const auto& name: species_names) {
this->speciesList.push_back(resolve_species_name(name));
}
}
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);
enum class EnginePolicy {
MAIN_SEQUENCE_POLICY
};
static const std::unordered_map<std::string, EnginePolicy> engine_map = {
{"MAIN_SEQUENCE_POLICY", EnginePolicy::MAIN_SEQUENCE_POLICY}
};
if (!engine_map.contains(policy_name)) {
throw gridfire::exceptions::PolicyError(
std::format(
"Engine Policy {} is not recognized. Valid policies are: {}",
policy_name,
gridfire::utils::iterable_to_delimited_string(engine_map, ", ", [](const auto& pair){ return pair.first; })
)
);
}
switch (engine_map.at(policy_name)) {
case EnginePolicy::MAIN_SEQUENCE_POLICY: {
this->policy = std::make_unique<gridfire::policy::MainSequencePolicy>(this->working_comp);
this->engine = &policy->construct();
break;
}
default:
throw gridfire::exceptions::PolicyError(
"Unhandled engine policy in GridFireContext::init_engine_from_policy"
);
}
}
void GridFireContext::init_solver_from_engine(const std::string &solver_name) {
enum class SolverType {
CVODE
};
static const std::unordered_map<std::string, SolverType> 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<gridfire::solver::CVODESolverStrategy>(*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) {
if (num_species == 0) {
throw fourdst::composition::exceptions::InvalidCompositionError("Cannot initialize composition with zero species.");
}
if (num_species != working_comp.size()) {
throw fourdst::composition::exceptions::InvalidCompositionError(
std::format(
"Number of species provided ({}) does not match the registered species count ({}).",
num_species,
working_comp.size()
)
);
}
for (size_t i = 0; i < num_species; i++) {
this->working_comp.setMolarAbundance(this->speciesList[i], abundances[i]);
}
}
int GridFireContext::evolve(
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& mass_lost
) {
init_composition_from_abundance_vector(Y_in, num_species);
gridfire::NetIn netIn;
netIn.temperature = T;
netIn.density = rho;
netIn.dt0 = dt0;
netIn.tMax = tMax;
netIn.composition = this->working_comp;
const gridfire::NetOut result = this->solver->evaluate(netIn);
energy_out = result.energy;
dEps_dT = result.dEps_dT;
dEps_dRho = result.dEps_dRho;
std::set<fourdst::atomic::Species> seen_species;
for (size_t i = 0; i < num_species; i++) {
fourdst::atomic::Species species = this->speciesList[i];
Y_out[i] = result.composition.getMolarAbundance(species);
seen_species.insert(species);
}
mass_lost = 0.0;
for (const auto& species : result.composition.getRegisteredSpecies()) {
if (!seen_species.contains(species)) {
mass_lost += species.mass() * result.composition.getMolarAbundance(species);
}
}
return 0;
}

293
src/extern/lib/gridfire_extern.cpp vendored Normal file
View File

@@ -0,0 +1,293 @@
#include "gridfire/gridfire.h"
#include "fourdst/composition/exceptions/exceptions_composition.h"
#include "gridfire/extern/gridfire_context.h"
#include "gridfire/extern/gridfire_extern.h"
extern "C" {
void* gf_init() {
return new GridFireContext();
}
void gf_free(void* ctx) {
delete static_cast<GridFireContext*>(ctx);
}
int gf_register_species(void* ptr, const int num_species, const char** species_names) {
auto* ctx = static_cast<GridFireContext*>(ptr);
try {
std::vector<std::string> names;
for(int i=0; i<num_species; ++i) {
names.emplace_back(species_names[i]);
}
ctx->init_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;
}
}
int gf_construct_engine_from_policy(
void* ptr,
const char* policy_name,
const double *abundances,
const size_t num_species
) {
auto* ctx = static_cast<GridFireContext*>(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<GridFireContext*>(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,
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* mass_lost
) {
auto* ctx = static_cast<GridFireContext*>(ptr);
try {
const int result = ctx->evolve(
Y,
num_species,
T,
rho,
tMax,
dt0,
Y_out,
*energy_out,
*dEps_dT,
*dEps_dRho, *mass_lost
);
if (result != 0) {
return result;
}
return GF_SUCCESS;
} catch (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 fourdst::composition::exceptions::InvalidCompositionError& e) {
ctx->last_error = e.what();
return FDSSE_INVALID_COMPOSITION_ERROR;
} catch (const fourdst::composition::exceptions::CompositionError& e) {
ctx->last_error = e.what();
return FDSSE_COMPOSITION_ERROR;
} catch (const gridfire::exceptions::InvalidQSESolutionError& e) {
ctx->last_error = e.what();
return GF_INVALID_QSE_SOLUTION_ERROR;
} catch (const gridfire::exceptions::FailedToPartitionEngineError& e) {
ctx->last_error = e.what();
return GF_FAILED_TO_PARTITION_ENGINE_ERROR;
} catch (const gridfire::exceptions::NetworkResizedError& e) {
ctx->last_error = e.what();
return GF_NETWORK_RESIZED_ERROR;
} catch (const gridfire::exceptions::UnableToSetNetworkReactionsError& e) {
ctx->last_error = e.what();
return GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR;
} catch (const gridfire::exceptions::BadCollectionError& e) {
ctx->last_error = e.what();
return GF_BAD_COLLECTION_ERROR;
} catch (const gridfire::exceptions::BadRHSEngineError& e) {
ctx->last_error = e.what();
return GF_BAD_RHS_ENGINE_ERROR;
} catch (const gridfire::exceptions::StaleJacobianError& e) {
ctx->last_error = e.what();
return GF_STALE_JACOBIAN_ERROR;
} catch (const gridfire::exceptions::UninitializedJacobianError& e) {
ctx->last_error = e.what();
return GF_UNINITIALIZED_JACOBIAN_ERROR;
} catch (const gridfire::exceptions::UnknownJacobianError& e) {
ctx->last_error = e.what();
return GF_UNKNOWN_JACOBIAN_ERROR;
} catch (const gridfire::exceptions::JacobianError& e) {
ctx->last_error = e.what();
return GF_JACOBIAN_ERROR;
} catch (const gridfire::exceptions::EngineError& e) {
ctx->last_error = e.what();
return GF_ENGINE_ERROR;
} catch (const gridfire::exceptions::ReactionParsingError& e) {
ctx->last_error = e.what();
return GF_REACTION_PARSING_ERROR;
} catch (const gridfire::exceptions::ReactionError& e) {
ctx->last_error = e.what();
return GF_REACTION_ERROR;
} catch (const gridfire::exceptions::SingularJacobianError& e) {
ctx->last_error = e.what();
return GF_SINGULAR_JACOBIAN_ERROR;
} catch (const gridfire::exceptions::IllConditionedJacobianError& e) {
ctx->last_error = e.what();
return GF_ILL_CONDITIONED_JACOBIAN_ERROR;
} catch (const gridfire::exceptions::CVODESolverFailureError& e) {
ctx->last_error = e.what();
return GF_CVODE_SOLVER_FAILURE_ERROR;
} catch (const gridfire::exceptions::KINSolSolverFailureError& e) {
ctx->last_error = e.what();
return GF_KINSOL_SOLVER_FAILURE_ERROR;
} catch (const gridfire::exceptions::SUNDIALSError& e) {
ctx->last_error = e.what();
return GF_SUNDIALS_ERROR;
} catch (const gridfire::exceptions::SolverError& e) {
ctx->last_error = e.what();
return GF_SOLVER_ERROR;
} catch (const gridfire::exceptions::HashingError& e) {
ctx->last_error = e.what();
return GF_HASHING_ERROR;
} catch (const gridfire::exceptions::UtilityError& e) {
ctx->last_error = e.what();
return GF_UTILITY_ERROR;
} catch (const gridfire::exceptions::DebugException& e) {
ctx->last_error = e.what();
return GF_DEBUG_ERROR;
} catch (const gridfire::exceptions::GridFireError& e) {
ctx->last_error = e.what();
return GF_GRIDFIRE_ERROR;
} catch (std::exception& e) {
ctx->last_error = e.what();
return GF_NON_GRIDFIRE_ERROR;
} catch (...) {
ctx->last_error = "Unknown error occurred during evolution.";
return GF_UNKNOWN_ERROR;
}
}
char* gf_get_last_error_message(void* ptr) {
const auto* ctx = static_cast<GridFireContext*>(ptr);
return const_cast<char*>(ctx->last_error.c_str());
}
char* gf_error_code_to_string(const int error_code) {
switch (error_code) {
case GF_SUCCESS:
return const_cast<char*>("GF_SUCCESS");
case GF_UNKNOWN_ERROR:
return const_cast<char*>("GF_UNKNOWN_ERROR");
case GF_NON_GRIDFIRE_ERROR:
return const_cast<char*>("GF_NON_GRIDFIRE_ERROR");
case GF_INVALID_QSE_SOLUTION_ERROR:
return const_cast<char*>("GF_INVALID_QSE_SOLUTION_ERROR");
case GF_FAILED_TO_PARTITION_ENGINE_ERROR:
return const_cast<char*>("GF_FAILED_TO_PARTITION_ENGINE_ERROR");
case GF_NETWORK_RESIZED_ERROR:
return const_cast<char*>("GF_NETWORK_RESIZED_ERROR");
case GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR:
return const_cast<char*>("GF_UNABLE_TO_SET_NETWORK_REACTIONS_ERROR");
case GF_BAD_COLLECTION_ERROR:
return const_cast<char*>("GF_BAD_COLLECTION_ERROR");
case GF_BAD_RHS_ENGINE_ERROR:
return const_cast<char*>("GF_BAD_RHS_ENGINE_ERROR");
case GF_STALE_JACOBIAN_ERROR:
return const_cast<char*>("GF_STALE_JACOBIAN_ERROR");
case GF_UNINITIALIZED_JACOBIAN_ERROR:
return const_cast<char*>("GF_UNINITIALIZED_JACOBIAN_ERROR");
case GF_UNKNOWN_JACOBIAN_ERROR:
return const_cast<char*>("GF_UNKNOWN_JACOBIAN_ERROR");
case GF_JACOBIAN_ERROR:
return const_cast<char*>("GF_JACOBIAN_ERROR");
case GF_ENGINE_ERROR:
return const_cast<char*>("GF_ENGINE_ERROR");
case GF_MISSING_BASE_REACTION_ERROR:
return const_cast<char*>("GF_MISSING_BASE_REACTION_ERROR");
case GF_MISSING_SEED_SPECIES_ERROR:
return const_cast<char*>("GF_MISSING_SEED_SPECIES_ERROR");
case GF_MISSING_KEY_REACTION_ERROR:
return const_cast<char*>("GF_MISSING_KEY_REACTION_ERROR");
case GF_POLICY_ERROR:
return const_cast<char*>("GF_POLICY_ERROR");
case GF_REACTION_PARSING_ERROR:
return const_cast<char*>("GF_REACTION_PARSING_ERROR");
case GF_REACTION_ERROR:
return const_cast<char*>("GF_REACTION_ERROR");
case GF_SINGULAR_JACOBIAN_ERROR:
return const_cast<char*>("GF_SINGULAR_JACOBIAN_ERROR");
case GF_ILL_CONDITIONED_JACOBIAN_ERROR:
return const_cast<char*>("GF_ILL_CONDITIONED_JACOBIAN_ERROR");
case GF_CVODE_SOLVER_FAILURE_ERROR:
return const_cast<char*>("GF_CVODE_SOLVER_FAILURE_ERROR");
case GF_KINSOL_SOLVER_FAILURE_ERROR:
return const_cast<char*>("GF_KINSOL_SOLVER_FAILURE_ERROR");
case GF_SUNDIALS_ERROR:
return const_cast<char*>("GF_SUNDIALS_ERROR");
case GF_SOLVER_ERROR:
return const_cast<char*>("GF_SOLVER_ERROR");
case GF_HASHING_ERROR:
return const_cast<char*>("GF_HASHING_ERROR");
case GF_UTILITY_ERROR:
return const_cast<char*>("GF_UTILITY_ERROR");
case GF_DEBUG_ERROR:
return const_cast<char*>("GF_DEBUG_ERROR");
case GF_GRIDFIRE_ERROR:
return const_cast<char*>("GF_GRIDFIRE_ERROR");
case FDSSE_NON_4DSTAR_ERROR:
return const_cast<char*>("FDSSE_NON_4DSTAR_ERROR");
case FDSSE_UNKNOWN_ERROR:
return const_cast<char*>("FDSSE_UNKNOWN_ERROR");
case FDSSE_SUCCESS:
return const_cast<char*>("FDSSE_SUCCESS");
case FDSSE_UNKNOWN_SYMBOL_ERROR:
return const_cast<char*>("FDSSE_UNKNOWN_SYMBOL_ERROR");
case FDSSE_SPECIES_ERROR:
return const_cast<char*>("FDSSE_SPECIES_ERROR");
case FDSSE_INVALID_COMPOSITION_ERROR:
return const_cast<char*>("FDSSE_INVALID_COMPOSITION_ERROR");
case FDSSE_COMPOSITION_ERROR:
return const_cast<char*>("FDSSE_COMPOSITION_ERROR");
default:
return const_cast<char*>("GF_UNRECOGNIZED_ERROR_CODE");
}
}
}

26
src/extern/meson.build vendored Normal file
View File

@@ -0,0 +1,26 @@
gridfire_extern_sources = files(
'lib/gridfire_context.cpp',
'lib/gridfire_extern.cpp',
)
gridfire_extern_dependencies = [
gridfire_dep
]
libgridfire_extern = library('gridfire_extern',
gridfire_extern_sources,
include_directories: include_directories('include'),
dependencies: gridfire_extern_dependencies,
install : true
)
gridfire_extern_dep = declare_dependency(
include_directories: include_directories('include'),
link_with: libgridfire_extern,
sources: gridfire_extern_sources,
dependencies: gridfire_extern_dependencies,
)
install_subdir('include/gridfire', install_dir: get_option('includedir'))
subdir('fortran')

2
src/meson.build
View File

@@ -60,6 +60,8 @@ gridfire_dep = declare_dependency(
install_subdir('include/gridfire', install_dir: get_option('includedir'))
subdir('extern')
if get_option('build-python')
message('Configuring Python bindings...')
subdir('python')

2
subprojects/fourdst.wrap
View File

@@ -1,4 +1,4 @@
[wrap-git]
url = https://github.com/4D-STAR/fourdst
revision = v0.9.6
revision = v0.9.7
depth = 1

80
tests/extern/C/gridfire_evolve.c vendored Normal file
View File

@@ -0,0 +1,80 @@
#include "gridfire/extern/gridfire_extern.h"
#include <stdio.h>
#define NUM_SPECIES 8
// 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); \
return 1; \
}
int main() {
void* ctx = gf_init();
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 abundances[NUM_SPECIES] = {
0.702616602672027,
9.74791583949078e-06,
0.06895512307276903,
0.00025,
7.855418029399437e-05,
0.0006014411598306529,
8.103062886768109e-05,
2.151340851063217e-05
};
int ret = gf_register_species(ctx, NUM_SPECIES, species_names);
CHECK_RET_CODE(ret, ctx, "SPECIES");
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");
CHECK_RET_CODE(ret, ctx, "CVODE");
double Y_out[NUM_SPECIES];
double energy_out;
double dEps_dT;
double dEps_dRho;
double mass_lost;
ret = gf_evolve(
ctx,
abundances,
NUM_SPECIES,
1.5e7, // Temperature in K
1.5e2, // Density in g/cm^3
3e17, // Time step in seconds
1e-12, // Initial time step in seconds
Y_out,
&energy_out,
&dEps_dT,
&dEps_dRho, &mass_lost
);
CHECK_RET_CODE(ret, ctx, "EVOLUTION");
printf("Evolved abundances:\n");
for (size_t i = 0; i < NUM_SPECIES; i++) {
printf("Species %s: %e\n", species_names[i], Y_out[i]);
}
printf("Energy output: %e\n", energy_out);
printf("dEps/dT: %e\n", dEps_dT);
printf("dEps/dRho: %e\n", dEps_dRho);
printf("Mass lost: %e\n", mass_lost);
gf_free(ctx);
return 0;
}

5
tests/extern/C/meson.build vendored Normal file
View File

@@ -0,0 +1,5 @@
executable('test_c_extern', 'gridfire_evolve.c',
install: false,
c_args: ['-Wall', '-Wextra'],
dependencies: [gridfire_extern_dep]
)

89
tests/extern/fortran/gridfire_evolve.f90 vendored Normal file
View File

@@ -0,0 +1,89 @@
program main
use iso_c_binding
use gridfire_mod
implicit none
type(GridFire) :: net
integer(c_int) :: ierr
integer :: i
! --- 1. Define Species and Initial Conditions ---
! Note: String lengths must match or exceed the longest name.
! We pad with spaces, which 'trim' handles inside the module.
character(len=5), dimension(8) :: species_names = [ &
"H-1 ", &
"He-3 ", &
"He-4 ", &
"C-12 ", &
"N-14 ", &
"O-16 ", &
"Ne-20", &
"Mg-24" &
]
! Initial Mass Fractions (converted to Molar Abundances Y = X/A)
! Standard solar-ish composition
real(c_double), dimension(8) :: Y_in = [ &
0.702616602672027, &
9.74791583949078e-06, &
0.06895512307276903, &
0.00025, &
7.855418029399437e-05, &
0.0006014411598306529, &
8.103062886768109e-05, &
2.151340851063217e-05 &
]
! Output buffers
real(c_double), dimension(8) :: Y_out
real(c_double) :: energy_out, dedt, dedrho, dmass
! Thermodynamic Conditions (Solar Core-ish)
real(c_double) :: T = 1.5e7 ! 15 Million K
real(c_double) :: rho = 150.0e0 ! 150 g/cm^3
real(c_double) :: dt = 3.0e17 ! 1 second timestep
! --- 2. Initialize GridFire ---
print *, "Initializing GridFire..."
call net%gff_init()
! --- 3. Register Species ---
print *, "Registering species..."
call net%register_species(species_names)
! --- 4. Configure Engine & Solver ---
print *, "Setting up Main Sequence Policy..."
call net%setup_policy("MAIN_SEQUENCE_POLICY", Y_in)
print *, "Setting up CVODE Solver..."
call net%setup_solver("CVODE")
! --- 5. Evolve ---
print *, "Evolving system (dt =", dt, "s)..."
call net%evolve(Y_in, T, rho, dt, Y_out, energy_out, dedt, dedrho, dmass, ierr)
if (ierr /= 0) then
print *, "Evolution Failed with error code: ", ierr
print *, "Error Message: ", net%get_last_error()
call net%gff_free() ! Always cleanup
stop
end if
! --- 6. Report Results ---
print *, ""
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)', "Mass Change: ", dmass
print *, ""
print *, "Abundances:"
do i = 1, size(species_names)
print '(A, " : ", ES12.5, " -> ", ES12.5)', &
trim(species_names(i)), Y_in(i), Y_out(i)
end do
! --- 7. Cleanup ---
call net%gff_free()
end program main

5
tests/extern/fortran/meson.build vendored Normal file
View File

@@ -0,0 +1,5 @@
executable('test_fortran_extern', 'gridfire_evolve.f90',
install: false,
fortran_args: ['-Wall', '-Wextra'],
dependencies: [gridfire_fortran_dep]
)

2
tests/extern/meson.build vendored Normal file
View File

@@ -0,0 +1,2 @@
subdir('C')
subdir('fortran')

1
tests/meson.build
View File

@@ -5,3 +5,4 @@ gtest_nomain_dep = dependency('gtest', main: false, required : true)
# Subdirectories for unit and integration tests
subdir('graphnet_sandbox')
subdir('extern')

95
validation/pynucastro/GridFireEquiv/GridFireEvolve.py
View File

@@ -1,95 +0,0 @@
from fourdst.composition import Composition
from gridfire.type import NetIn
from gridfire.policy import MainSequencePolicy
from gridfire.solver import CVODESolverStrategy
from enum import Enum
from typing import Dict, Union, SupportsFloat
import json
import dicttoxml
def init_composition() -> Composition:
Y = [7.0262E-01, 9.7479E-06, 6.8955E-02, 2.5000E-04, 7.8554E-05, 6.0144E-04, 8.1031E-05, 2.1513E-05]
S = ["H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"]
return Composition(S, Y)
def init_netIn(temp: float, rho: float, time: float, comp: Composition) -> NetIn:
netIn = NetIn()
netIn.temperature = temp
netIn.density = rho
netIn.tMax = time
netIn.dt0 = 1e-12
netIn.composition = comp
return netIn
class StepData(Enum):
TIME = 0
DT = 1
COMP = 2
CONTRIB = 3
class StepLogger:
def __init__(self):
self.num_steps: int = 0
self.step_data: Dict[int, Dict[StepData, Union[SupportsFloat, Dict[str, SupportsFloat]]]] = {}
def log_step(self, context):
engine = context.engine
self.step_data[self.num_steps] = {}
self.step_data[self.num_steps][StepData.TIME] = context.t
self.step_data[self.num_steps][StepData.DT] = context.dt
comp_data: Dict[str, SupportsFloat] = {}
for species in engine.getNetworkSpecies():
sid = engine.getSpeciesIndex(species)
comp_data[species.name()] = context.state[sid]
self.step_data[self.num_steps][StepData.COMP] = comp_data
self.num_steps += 1
def to_json (self, filename: str):
serializable_data = {
stepNum: {
StepData.TIME.name: step[StepData.TIME],
StepData.DT.name: step[StepData.DT],
StepData.COMP.name: step[StepData.COMP],
}
for stepNum, step in self.step_data.items()
}
with open(filename, 'w') as f:
json.dump(serializable_data, f, indent=4)
def to_xml(self, filename: str):
serializable_data = {
stepNum: {
StepData.TIME.name: step[StepData.TIME],
StepData.DT.name: step[StepData.DT],
StepData.COMP.name: step[StepData.COMP],
}
for stepNum, step in self.step_data.items()
}
xml_data = dicttoxml.dicttoxml(serializable_data, custom_root='StepLog', attr_type=False)
with open(filename, 'wb') as f:
f.write(xml_data)
def main(temp: float, rho: float, time: float):
comp = init_composition()
netIn = init_netIn(temp, rho, time, comp)
policy = MainSequencePolicy(comp)
engine = policy.construct()
solver = CVODESolverStrategy(engine)
step_logger = StepLogger()
solver.set_callback(lambda context: step_logger.log_step(context))
solver.evaluate(netIn, False)
step_logger.to_xml("log_data.xml")
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description="Simple python example of GridFire usage")
parser.add_argument("-t", "--temp", type=float, help="Temperature in K", default=1.5e7)
parser.add_argument("-r", "--rho", type=float, help="Density in g/cm^3", default=1.5e2)
parser.add_argument("--tMax", type=float, help="Time in s", default=3.1536 * 1e17)
args = parser.parse_args()
main(args.temp, args.rho, args.tMax)

122
validation/pynucastro/GridFireEquiv/evolve.py
View File

@@ -1,122 +0,0 @@
# import numpy as np
# from scipy.integrate import solve_ivp
#
# import pp_chain_robust as network
#
#
# T = 1e7
# rho = 1.5e2
#
# Y0 = np.zeros(network.nnuc)
# Y0[network.jp] = 0.702583
# Y0[network.jhe3] = 9.74903e-6
# Y0[network.jhe4] = 0.068963
# Y0[network.jc12] = 0.000250029
# Y0[network.jn14] = 7.85632e-5
# Y0[network.jo16] = 0.00060151
# Y0[network.jne20] = 8.10399e-5
# Y0[network.jmg24] = 2.15159e-5
#
# t_start = 0.0
# t_end = 3.14e17
# t_span = (t_start, t_end)
#
# sol = solve_ivp(
# network.rhs,
# t_span=t_span,
# y0=Y0,
# method='Radau',
# jac=network.jacobian,
# args=(rho, T),
# dense_output=True,
# rtol=1e-8,
# atol=1e-8
# )
#
#
# with open("pynucastro_results.csv", 'w') as f:
# f.write("t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24\n")
# for (t,h1,h2,he3,he4,c12,n14,o16,ne20,mg24) in zip(sol.t, sol.y[network.jp, :], sol.y[network.jd, :], sol.y[network.jhe3, :], sol.y[network.jhe4, :], sol.y[network.jc12, :], sol.y[network.jn14, :], sol.y[network.jo16, :], sol.y[network.jne20, :], sol.y[network.jmg24, :]):
# f.write(f"{t},{h1},{h2},{he3},{he4},{c12},{n14},{o16},{ne20},{mg24}\n")
import sys
from gridfire.engine import GraphEngine, MultiscalePartitioningEngineView, AdaptiveEngineView
from gridfire.solver import CVODESolverStrategy
from gridfire.type import NetIn
from fourdst.composition import Composition
from fourdst.atomic import species, Species
from datetime import datetime
from typing import List, Dict, Set, Tuple
X : list[float] = [0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4]
symbols : list[str] = ["H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"]
comp = Composition()
comp.registerSymbol(symbols)
comp.setMassFraction(symbols, X)
comp.finalize(True)
print(f"Initial H-1 mass fraction {comp.getMassFraction("H-1")}")
netIn = NetIn()
netIn.composition = comp
netIn.temperature = 1.5e7
netIn.density = 1.6e2
netIn.tMax = 3e17
netIn.dt0 = 1e-12
baseEngine = GraphEngine(netIn.composition, 2)
baseEngine.setUseReverseReactions(False)
qseEngine = MultiscalePartitioningEngineView(baseEngine)
adaptiveEngine = AdaptiveEngineView(qseEngine)
solver = CVODESolverStrategy(adaptiveEngine)
data: List[Tuple[float, Dict[str, Tuple[float, float]]]] = []
def callback(context):
engine = context.engine
abundances: Dict[str, Tuple[float, float]] = {}
for species in engine.getNetworkSpecies():
sid = engine.getSpeciesIndex(species)
abundances[species.name()] = (species.mass(), context.state[sid])
data.append((context.t,abundances))
solver.set_callback(callback)
try:
results = solver.evaluate(netIn, False)
print(f"H-1 final molar abundance: {results.composition.getMolarAbundance("H-1"):0.3f}")
except:
print("Warning: solver did not converge", file=sys.stderr)
uniqueSpecies: Set[Tuple[str, float]] = set()
for t, timestep in data:
for (name, (mass, abundance)) in timestep.items():
uniqueSpecies.add((name, mass))
sortedUniqueSpecies = list(sorted(uniqueSpecies, key=lambda e: e[1]))
with open(f"gridfire_results_{datetime.now().strftime("%H-%M-%d_%m_%Y")}.csv", 'w') as f:
f.write('t,')
for i, (species,mass) in enumerate(sortedUniqueSpecies):
f.write(f"{species}")
if i < len(sortedUniqueSpecies)-1:
f.write(",")
f.write("\n")
for t, timestep in data:
f.write(f"{t},")
for i, (species, mass) in enumerate(sortedUniqueSpecies):
if timestep.get(species, None) is not None:
f.write(f"{timestep.get(species)[1]}")
else:
f.write(f"")
if i < len(sortedUniqueSpecies) - 1:
f.write(",")
f.write("\n")

46
validation/pynucastro/GridFireEquiv/generate_proton_proton.py
View File

@@ -1,46 +0,0 @@
import pynucastro as pyna
from gridfire.engine import GraphEngine, MultiscalePartitioningEngineView, AdaptiveEngineView
from gridfire.type import NetIn
from fourdst.composition import Composition
symbols : list[str] = ["H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"]
X : list[float] = [0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4]
comp = Composition()
comp.registerSymbol(symbols)
comp.setMassFraction(symbols, X)
comp.finalize(True)
print(f"Initial H-1 mass fraction {comp.getMassFraction("H-1")}")
netIn = NetIn()
netIn.composition = comp
netIn.temperature = 1e7
netIn.density = 1.6e2
netIn.tMax = 1e-9
netIn.dt0 = 1e-12
baseEngine = GraphEngine(netIn.composition, 2)
equiv_species = baseEngine.getNetworkSpecies()
equiv_species = [x.name().replace("-","").lower() for x in equiv_species]
equiv_species = [x if x != 'h1' else 'p' for x in equiv_species]
equiv_species = [x if x != 'n1' else 'n' for x in equiv_species]
rl = pyna.ReacLibLibrary()
# equiv_species = ['p', 'd', 'he3', 'he4', 'c12', 'n14', 'o16', 'ne20', 'mg24']
full_lib = rl.linking_nuclei(equiv_species)
print(f"\nFound {len(full_lib.get_rates())} rates.")
net = pyna.PythonNetwork(libraries=[full_lib])
output_filename = "pp_chain_robust.py"
net.write_network(output_filename)
print(f"\nSuccessfully wrote robust PP-chain network to: {output_filename}")

17677
validation/pynucastro/GridFireEquiv/pp_chain_robust.py
View File

File diff suppressed because it is too large Load Diff

6864
validation/pynucastro/GridFireEquiv/pynucastro.py
View File

File diff suppressed because it is too large Load Diff

102
validation/vv/GridFireValidationSuite.py Normal file
View File

@@ -0,0 +1,102 @@
from gridfire.policy import MainSequencePolicy, NetworkPolicy
from gridfire.engine import DynamicEngine, GraphEngine
from gridfire.type import NetIn
from fourdst.composition import Composition
from testsuite import TestSuite
from utils import init_netIn, init_composition, years_to_seconds
from enum import Enum
class SolarLikeStar_QSE_Suite(TestSuite):
def __init__(self):
initialComposition : Composition = init_composition()
super().__init__(
name="SolarLikeStar_QSE",
description="GridFire simulation of a roughly solar like star over 10Gyr with QSE enabled.",
temp=1.5e7,
density=1.5e2,
tMax=years_to_seconds(1e10),
composition=initialComposition,
notes="Thermodynamically Static, MultiscalePartitioning Engine View"
)
def __call__(self):
policy : MainSequencePolicy = MainSequencePolicy(self.composition)
engine : DynamicEngine = policy.construct()
netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition)
self.evolve(engine, netIn)
class MetalEnhancedSolarLikeStar_QSE_Suite(TestSuite):
def __init__(self):
initialComposition : Composition = init_composition(ZZs=1)
super().__init__(
name="MetalEnhancedSolarLikeStar_QSE",
description="GridFire simulation of a star with solar core temp and density but enhanced by 1 dex in Z.",
temp=0.8 * 1.5e7,
density=1.5e2,
tMax=years_to_seconds(1e10),
composition=initialComposition,
notes="Thermodynamically Static, MultiscalePartitioning Engine View, Z enhanced by 1 dex, temperature reduced to 80% of solar core"
)
def __call__(self):
policy : MainSequencePolicy = MainSequencePolicy(self.composition)
engine : GraphEngine = policy.construct()
netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition)
self.evolve(engine, netIn)
class MetalDepletedSolarLikeStar_QSE_Suite(TestSuite):
def __init__(self):
initialComposition : Composition = init_composition(ZZs=-1)
super().__init__(
name="MetalDepletedSolarLikeStar_QSE",
description="GridFire simulation of a star with solar core temp and density but depleted by 1 dex in Z.",
temp=1.2 * 1.5e7,
density=1.5e2,
tMax=years_to_seconds(1e10),
composition=initialComposition,
notes="Thermodynamically Static, MultiscalePartitioning Engine View, Z depleted by 1 dex, temperature increased to 120% of solar core"
)
def __call__(self):
policy : MainSequencePolicy = MainSequencePolicy(self.composition)
engine : GraphEngine = policy.construct()
netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition)
self.evolve(engine, netIn)
class SolarLikeStar_No_QSE_Suite(TestSuite):
def __init__(self):
initialComposition : Composition = init_composition()
super().__init__(
name="SolarLikeStar_No_QSE",
description="GridFire simulation of a roughly solar like star over 10Gyr with QSE disabled.",
temp=1.5e7,
density=1.5e2,
tMax=years_to_seconds(1e10),
composition=initialComposition,
notes="Thermodynamically Static, No MultiscalePartitioning Engine View"
)
def __call__(self):
engine : GraphEngine = GraphEngine(self.composition, 3)
netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition)
self.evolve(engine, netIn)
class ValidationSuites(Enum):
SolarLikeStar_QSE = SolarLikeStar_QSE_Suite
SolarLikeStar_No_QSE = SolarLikeStar_No_QSE_Suite
MetalDepletedSolarLikeStar_QSE = MetalDepletedSolarLikeStar_QSE_Suite
MetalEnhancedSolarLikeStar_QSE = MetalEnhancedSolarLikeStar_QSE_Suite
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description="Run some subset of the GridFire validation suite.")
parser.add_argument('--suite', type=str, choices=[suite.name for suite in ValidationSuites], nargs="+", help="The validation suite to run.")
args = parser.parse_args()
for suite_name in args.suite:
suite = ValidationSuites[suite_name]
instance : TestSuite = suite.value()
instance()

77
validation/vv/logger.py Normal file
View File

@@ -0,0 +1,77 @@
from enum import Enum
from typing import Dict, List, Any, SupportsFloat
import json
from datetime import datetime
import os
import sys
from gridfire.solver import CVODETimestepContext
import gridfire
class LogEntries(Enum):
Step = "Step"
t = "t"
dt = "dt"
eps = "eps"
Composition = "Composition"
ReactionContributions = "ReactionContributions"
class StepLogger:
def __init__(self):
self.num_steps : int = 0
self.steps : List[Dict[LogEntries, Any]] = []
def log_step(self, ctx : CVODETimestepContext):
comp_data: Dict[str, SupportsFloat] = {}
for species in ctx.engine.getNetworkSpecies():
sid = ctx.engine.getSpeciesIndex(species)
comp_data[species.name()] = ctx.state[sid]
entry : Dict[LogEntries, Any] = {
LogEntries.Step: ctx.num_steps,
LogEntries.t: ctx.t,
LogEntries.dt: ctx.dt,
LogEntries.Composition: comp_data,
}
self.steps.append(entry)
self.num_steps += 1
def to_json(self, filename: str, **kwargs):
serializable_steps : List[Dict[str, Any]] = [
{
LogEntries.Step.value: step[LogEntries.Step],
LogEntries.t.value: step[LogEntries.t],
LogEntries.dt.value: step[LogEntries.dt],
LogEntries.Composition.value: step[LogEntries.Composition],
}
for step in self.steps
]
out_data : Dict[str, Any] = {
"Metadata": {
"NumSteps": self.num_steps,
**kwargs,
"DateCreated": datetime.now().isoformat(),
"GridFireVersion": gridfire.__version__,
"Author": "Emily M. Boudreaux",
"OS": os.uname().sysname,
"ClangVersion": os.popen("clang --version").read().strip(),
"GccVersion": os.popen("gcc --version").read().strip(),
"PythonVersion": sys.version,
},
"Steps": serializable_steps
}
with open(filename, 'w') as f:
json.dump(out_data, f, indent=4)
def summary(self) -> Dict[str, Any]:
if not self.steps:
return {}
final_step = self.steps[-1]
summary_data : Dict[str, Any] = {
"TotalSteps": self.num_steps,
"FinalTime": final_step[LogEntries.t],
"FinalComposition": final_step[LogEntries.Composition],
}
return summary_data

241
validation/vv/testsuite.py Normal file
View File

@@ -0,0 +1,241 @@
from abc import ABC, abstractmethod
import fourdst.atomic
import scipy.integrate
from fourdst.composition import Composition
from gridfire.engine import DynamicEngine, GraphEngine
from gridfire.type import NetIn, NetOut
from gridfire.exceptions import GridFireError
from gridfire.solver import CVODESolverStrategy
from logger import StepLogger
from typing import List
import re
from typing import Dict, Tuple, Any, Union
from datetime import datetime
import pynucastro as pyna
import os
import importlib.util
import sys
import numpy as np
import json
import time
def load_network_module(filepath):
module_name = os.path.basename(filepath).replace(".py", "")
if module_name in sys.modules: # clear any existing module with the same name
del sys.modules[module_name]
spec = importlib.util.spec_from_file_location(module_name, filepath)
if spec is None:
raise FileNotFoundError(f"Could not find module at {filepath}")
network_module = importlib.util.module_from_spec(spec)
sys.modules[module_name] = network_module
spec.loader.exec_module(network_module)
return network_module
def get_pyna_rate(my_rate_str, library):
match = re.match(r"([a-zA-Z0-9]+)\(([^,]+),([^)]*)\)(.*)", my_rate_str)
if not match:
print(f"Could not parse string format: {my_rate_str}")
return None
target = match.group(1)
projectile = match.group(2)
ejectiles = match.group(3)
product = match.group(4)
def expand_species(s_str):
if not s_str or s_str.strip() == "":
return []
# Split by space (handling "p a" or "2p a")
parts = s_str.split()
expanded = []
for p in parts:
# Check for multipliers like 2p, 3a
mult_match = re.match(r"(\d+)([a-zA-Z0-9]+)", p)
if mult_match:
count = int(mult_match.group(1))
spec = mult_match.group(2)
# Map common aliases if necessary (though pyna handles most)
if spec == 'a': spec = 'he4'
expanded.extend([spec] * count)
else:
spec = p
if spec == 'a': spec = 'he4'
expanded.append(spec)
return expanded
reactants_str = [target] + expand_species(projectile)
products_str = expand_species(ejectiles) + [product]
# Convert strings to pyna.Nucleus objects
try:
r_nuc = [pyna.Nucleus(r) for r in reactants_str]
p_nuc = [pyna.Nucleus(p) for p in products_str]
except Exception as e:
print(f"Error converting nuclei for {my_rate_str}: {e}")
return None
rates = library.get_rate_by_nuclei(r_nuc, p_nuc)
if rates:
if isinstance(rates, list):
return rates[0] # Return the first match
return rates
else:
return None
class TestSuite(ABC):
def __init__(self, name: str, description: str, temp: float, density: float, tMax: float, composition: Composition, notes: str = ""):
self.name : str = name
self.description : str = description
self.temperature : float = temp
self.density : float = density
self.tMax : float = tMax
self.composition : Composition = composition
self.notes : str = notes
def evolve_pynucastro(self, engine: GraphEngine):
print("Evolution complete. Now building equivalent pynucastro network...")
# Build equivalent pynucastro network for comparison
reaclib_library : pyna.ReacLibLibrary = pyna.ReacLibLibrary()
rate_names = [r.id().replace("e+","").replace("e-","").replace(", ", ",") for r in engine.getNetworkReactions()]
goodRates : List[pyna.rates.reaclib_rate.ReacLibRate] = []
missingRates = []
for r_str in rate_names:
# Try the exact name match first (fastest)
try:
pyna_rate = reaclib_library.get_rate_by_name(r_str)
if isinstance(pyna_rate, list):
goodRates.append(pyna_rate[0])
else:
goodRates.append(pyna_rate)
except:
# Fallback to the smart parser
pyna_rate = get_pyna_rate(r_str, reaclib_library)
if pyna_rate:
goodRates.append(pyna_rate)
else:
missingRates.append(r_str)
pynet : pyna.PythonNetwork = pyna.PythonNetwork(rates=goodRates)
pynet.write_network(f"{self.name}_pynucastro_network.py")
net = load_network_module(f"{self.name}_pynucastro_network.py")
Y0 = np.zeros(net.nnuc)
Y0[net.jp] = self.composition.getMolarAbundance("H-1")
Y0[net.jhe3] = self.composition.getMolarAbundance("He-3")
Y0[net.jhe4] = self.composition.getMolarAbundance("He-4")
Y0[net.jc12] = self.composition.getMolarAbundance("C-12")
Y0[net.jn14] = self.composition.getMolarAbundance("N-14")
Y0[net.jo16] = self.composition.getMolarAbundance("O-16")
Y0[net.jne20] = self.composition.getMolarAbundance("Ne-20")
Y0[net.jmg24] = self.composition.getMolarAbundance("Mg-24")
print("Starting pynucastro integration...")
startTime = time.time()
sol = scipy.integrate.solve_ivp(
net.rhs,
[0, self.tMax],
Y0,
args=(self.density, self.temperature),
method="BDF",
jac=net.jacobian,
rtol=1e-5,
atol=1e-8
)
endTime = time.time()
print("Pynucastro integration complete. Writing results to JSON...")
data: List[Dict[str, Union[float, Dict[str, float]]]] = []
for time_step, t in enumerate(sol.t):
data.append({"t": t, "Composition": {}})
for j in range(net.nnuc):
A = net.A[j]
Z = net.Z[j]
species: str
try:
species = fourdst.atomic.az_to_species(A, Z).name()
except:
species = f"SP-A_{A}_Z_{Z}"
data[-1]["Composition"][species] = sol.y[j, time_step]
pynucastro_json : Dict[str, Any] = {
"Metadata": {
"Name": f"{self.name}_pynucastro",
"Description": f"pynucastro simulation equivalent to GridFire validation suite: {self.description}",
"Status": "Success",
"Notes": self.notes,
"Temperature": self.temperature,
"Density": self.density,
"tMax": self.tMax,
"ElapsedTime": endTime - startTime,
"DateCreated": datetime.now().isoformat()
},
"Steps": data
}
with open(f"GridFireValidationSuite_{self.name}_pynucastro.json", "w") as f:
json.dump(pynucastro_json, f, indent=4)
def evolve(self, engine: GraphEngine, netIn: NetIn, pynucastro_compare: bool = True):
solver : CVODESolverStrategy = CVODESolverStrategy(engine)
stepLogger : StepLogger = StepLogger()
solver.set_callback(lambda ctx: stepLogger.log_step(ctx))
startTime = time.time()
try:
startTime = time.time()
netOut : NetOut = solver.evaluate(netIn)
endTime = time.time()
stepLogger.to_json(
f"GridFireValidationSuite_{self.name}_OKAY.json",
Name = f"{self.name}_Success",
Description=self.description,
Status="Success",
Notes=self.notes,
Temperature=netIn.temperature,
Density=netIn.density,
tMax=netIn.tMax,
FinalEps = netOut.energy,
FinaldEpsdT = netOut.dEps_dT,
FinaldEpsdRho = netOut.dEps_dRho,
ElapsedTime = endTime - startTime
)
except GridFireError as e:
endTime = time.time()
stepLogger.to_json(
f"GridFireValidationSuite_{self.name}_FAIL.json",
Name = f"{self.name}_Failure",
Description=self.description,
Status=f"Error",
ErrorMessage=str(e),
Notes=self.notes,
Temperature=netIn.temperature,
Density=netIn.density,
tMax=netIn.tMax,
ElapsedTime = endTime - startTime
)
if pynucastro_compare:
self.evolve_pynucastro(engine)
@abstractmethod
def __call__(self):
pass

56
validation/vv/utils.py Normal file
View File

@@ -0,0 +1,56 @@
from fourdst.composition import Composition
from fourdst.composition import CanonicalComposition
from fourdst.atomic import Species
from gridfire.type import NetIn
def rescale_composition(comp_ref : Composition, ZZs : float, Y_primordial : float = 0.248) -> Composition:
CC : CanonicalComposition = comp_ref.getCanonicalComposition()
dY_dZ = (CC.Y - Y_primordial) / CC.Z
Z_new = CC.Z * (10**ZZs)
Y_bulk_new = Y_primordial + (dY_dZ * Z_new)
X_new = 1.0 - Z_new - Y_bulk_new
if X_new < 0: raise ValueError(f"ZZs={ZZs} yields unphysical composition (X < 0)")
ratio_H = X_new / CC.X if CC.X > 0 else 0
ratio_He = Y_bulk_new / CC.Y if CC.Y > 0 else 0
ratio_Z = Z_new / CC.Z if CC.Z > 0 else 0
Y_new_list = []
newComp : Composition = Composition()
s: Species
for s in comp_ref.getRegisteredSpecies():
Xi_ref = comp_ref.getMassFraction(s)
if s.el() == "H":
Xi_new = Xi_ref * ratio_H
elif s.el() == "He":
Xi_new = Xi_ref * ratio_He
else:
Xi_new = Xi_ref * ratio_Z
Y = Xi_new / s.mass()
newComp.registerSpecies(s)
newComp.setMolarAbundance(s, Y)
return newComp
def init_composition(ZZs : float = 0) -> Composition:
Y_solar = [7.0262E-01, 9.7479E-06, 6.8955E-02, 2.5000E-04, 7.8554E-05, 6.0144E-04, 8.1031E-05, 2.1513E-05]
S = ["H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"]
return rescale_composition(Composition(S, Y_solar), ZZs)
def init_netIn(temp: float, rho: float, time: float, comp: Composition) -> NetIn:
n : NetIn = NetIn()
n.temperature = temp
n.density = rho
n.tMax = time
n.dt0 = 1e-12
n.composition = comp
return n
def years_to_seconds(years: float) -> float:
return years * 3.1536e7