From d65c237b26cf45550095328d57b92fb8fa3a391a Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 19 Dec 2025 09:58:47 -0500 Subject: [PATCH] feat(fortran): Fortran interface can now use multi-zone Fortran interface uses the new C api ability to call the naieve multi-zone solver. This allows fortran calling code to make use of in build parellaism for solving multiple zones --- build-check/FC/meson.build | 7 +- src/extern/fortran/gridfire_mod.f90 | 201 +++++++-- .../include/gridfire/extern/gridfire_extern.h | 4 +- src/extern/lib/gridfire_extern.cpp | 24 +- tests/extern/C/gridfire_evolve_multi.c | 4 +- tests/extern/C/gridfire_evolve_single.c | 2 +- .../extern/fortran/gridfire_evolve_multi.f90 | 134 ++++++ ..._evolve.f90 => gridfire_evolve_single.f90} | 17 +- tests/extern/fortran/meson.build | 10 +- tools/cli/meson.build | 1 - .../generate_config_files.cpp | 0 tools/{config => gf_config}/meson.build | 0 tools/gf_multi/main.cpp | 393 ++++++++++++++++++ tools/gf_multi/meson.build | 1 + tools/{cli => }/gf_quick/main.cpp | 4 + tools/{cli => }/gf_quick/meson.build | 0 tools/meson.build | 5 +- 17 files changed, 738 insertions(+), 69 deletions(-) create mode 100644 tests/extern/fortran/gridfire_evolve_multi.f90 rename tests/extern/fortran/{gridfire_evolve.f90 => gridfire_evolve_single.f90} (82%) delete mode 100644 tools/cli/meson.build rename tools/{config => gf_config}/generate_config_files.cpp (100%) rename tools/{config => gf_config}/meson.build (100%) create mode 100644 tools/gf_multi/main.cpp create mode 100644 tools/gf_multi/meson.build rename tools/{cli => }/gf_quick/main.cpp (98%) rename tools/{cli => }/gf_quick/meson.build (100%) diff --git a/build-check/FC/meson.build b/build-check/FC/meson.build index 5d1cc40e..6366845d 100644 --- a/build-check/FC/meson.build +++ b/build-check/FC/meson.build @@ -1,5 +1,10 @@ if get_option('build_fortran') - add_languages('fortran', native: true) + found_fortran = add_languages('fortran') + if not found_fortran + error('Fortran compiler not found, but build_fortran option is enabled.') + else + message('Fortran compiler found.') + endif message('Found FORTRAN compiler: ' + meson.get_compiler('fortran').get_id()) message('Fortran standard set to: ' + get_option('fortran_std')) message('Building fortran module (gridfire_mod.mod)') diff --git a/src/extern/fortran/gridfire_mod.f90 b/src/extern/fortran/gridfire_mod.f90 index b9fabf53..38d0889b 100644 --- a/src/extern/fortran/gridfire_mod.f90 +++ b/src/extern/fortran/gridfire_mod.f90 @@ -2,6 +2,14 @@ module gridfire_mod use iso_c_binding implicit none + type, public :: GF_TYPE + integer(c_int) :: value + end type GF_TYPE + + type(GF_TYPE), parameter, public :: & + SINGLE_ZONE = GF_TYPE(1001), & + MULTI_ZONE = GF_TYPE(1002) + enum, bind (C) enumerator :: FDSSE_NON_4DSTAR_ERROR = -102 enumerator :: FDSSE_UNKNOWN_ERROR = -101 @@ -50,24 +58,46 @@ module gridfire_mod enumerator :: GF_DEBUG_ERRROR = 30 enumerator :: GF_GRIDFIRE_ERROR = 31 + enumerator :: GF_UNINITIALIZED_INPUT_MEMORY_ERROR = 32 + enumerator :: GF_UNINITIALIZED_OUTPUT_MEMORY_ERROR = 33 + + enumerator :: GF_INVALD_NUM_SPECIES = 34 + enumerator :: GF_INVALID_TIMESTEPS = 35 + enumerator :: GF_UNKNONWN_FREE_TYPE = 36 + + enumerator :: GF_INVALID_TYPE = 37 + + enumerator :: GF_SINGLE_ZONE = 1001 + enumerator :: GF_MULTI_ZONE = 1002 end enum interface ! void* gf_init() - function gf_init() bind(C, name="gf_init") - import :: c_ptr + function gf_init(ctx_type) bind(C, name="gf_init") + import :: c_ptr, c_int type(c_ptr) :: gf_init + integer(c_int), value :: ctx_type 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 + ! int gf_free(void* gf) + function gf_free(ctx_type, ptr) result(c_res) bind(C, name="gf_free") + import :: c_ptr, c_int + type(c_ptr), value :: ptr + integer(c_int), value :: ctx_type + integer(c_int) :: c_res + end function gf_free + + function gf_set_num_zones(ctx_type, ptr, num_zones) result(c_res) bind(C, name="gf_set_num_zones") + import :: c_ptr, c_int, c_size_t + type(c_ptr), value :: ptr + integer(c_int), value :: ctx_type + integer(c_size_t), value :: num_zones + integer(c_int) :: c_res + end function gf_set_num_zones ! 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 + import :: c_ptr, c_int type(c_ptr), value :: ptr type(c_ptr) :: c_msg end function @@ -102,49 +132,116 @@ module gridfire_mod end function ! int gf_evolve(...) - function gf_evolve(ptr, Y_in, num_species, T, rho, dt, Y_out, energy_out, dEps_dT, dEps_dRho, specific_neutrino_loss, specific_neutrino_flux, mass_lost) result(ierr) & + function gf_evolve_c_scalar(ctx_type, ptr, Y_in, num_species, T, rho, tMax, dt0, & + Y_out, energy, dedt, dedrho, & + nue_loss, nu_flux, mass_lost) result(ierr) & bind(C, name="gf_evolve") - import + import :: c_ptr, c_int, c_double, c_size_t type(c_ptr), value :: ptr - real(c_double), dimension(*), intent(in) :: Y_in + integer(c_int), value :: ctx_type integer(c_size_t), value :: num_species - real(c_double), value :: T, rho, dt + + ! Arrays + real(c_double), dimension(*), intent(in) :: Y_in real(c_double), dimension(*), intent(out) :: Y_out - real(c_double), intent(out) :: energy_out, dEps_dT, dEps_dRho, specific_neutrino_loss, specific_neutrino_flux, mass_lost + + ! Scalars (Passed by Reference -> matches void*) + real(c_double), intent(in) :: T, rho + real(c_double), intent(out) :: energy, dedt, dedrho, nue_loss, nu_flux, mass_lost + + ! Scalars (Passed by Value) + real(c_double), value :: tMax, dt0 + + integer(c_int) :: ierr + end function + + ! 2. Interface for Multi Zone (Arrays) + function gf_evolve_c_array(ctx_type, ptr, Y_in, num_species, T, rho, tMax, dt0, & + Y_out, energy, dedt, dedrho, & + nue_loss, nu_flux, mass_lost) result(ierr) & + bind(C, name="gf_evolve") + import :: c_ptr, c_int, c_double, c_size_t + type(c_ptr), value :: ptr + integer(c_int), value :: ctx_type + integer(c_size_t), value :: num_species + + ! All Arrays (dimension(*)) + real(c_double), dimension(*), intent(in) :: Y_in + real(c_double), dimension(*), intent(in) :: T, rho + + real(c_double), dimension(*), intent(out) :: Y_out + real(c_double), dimension(*), intent(out) :: energy, dedt, dedrho, nue_loss, nu_flux, mass_lost + + ! Scalars (Passed by Value) + real(c_double), value :: tMax, dt0 + integer(c_int) :: ierr end function end interface type :: GridFire type(c_ptr) :: ctx = c_null_ptr + integer(c_int) :: ctx_type = SINGLE_ZONE%value integer(c_size_t) :: num_species = 0 + integer(c_size_t) :: num_zones = 1 contains procedure :: gff_init procedure :: gff_free - procedure :: register_species - procedure :: setup_policy - procedure :: setup_solver - procedure :: evolve - procedure :: get_last_error + procedure :: gff_register_species + procedure :: gff_setup_policy + procedure :: gff_setup_solver + procedure :: gff_get_last_error + + procedure :: gff_evolve_single + procedure :: gff_evolve_multi + + generic :: gff_evolve => gff_evolve_single, gff_evolve_multi end type GridFire contains - subroutine gff_init(self) + subroutine gff_init(self, type, zones) class(GridFire), intent(out) :: self + type(GF_TYPE), intent(in) :: type + integer(c_size_t), intent(in), optional :: zones + integer(c_int) :: ierr - self%ctx = gf_init() + if (type%value==1002) then + if (.not. present(zones)) then + print *, "GridFire Error: Multi-zone type requires number of zones to be specficied in the GridFire init method (i.e. GridFire(MULTI_ZONE, 10) for 10 zones)." + error stop + end if + + self%num_zones = zones + end if + + self%ctx_type = type%value + + self%ctx = gf_init(self%ctx_type) + + if (type%value==1002) then + ierr = gf_set_num_zones(self%ctx_type, self%ctx, self%num_zones) + if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then + print *, "GridFire Multi-Zone Error: ", self%gff_get_last_error() + error stop + end if + end if end subroutine gff_init subroutine gff_free(self) class(GridFire), intent(inout) :: self + integer(c_int) :: ierr if (c_associated(self%ctx)) then - call gf_free(self%ctx) + ierr = gf_free(self%ctx_type, self%ctx) + if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then + print *, "GridFire Free Error: ", self%gff_get_last_error() + error stop + end if self%ctx = c_null_ptr end if end subroutine gff_free - function get_last_error(self) result(msg) + function gff_get_last_error(self) result(msg) class(GridFire), intent(in) :: self character(len=:), allocatable :: msg type(c_ptr) :: c_msg_ptr @@ -169,9 +266,9 @@ module gridfire_mod do i = 1, len_str msg(i+10:i+10) = char_ptr(i) end do - end function get_last_error + end function gff_get_last_error - subroutine register_species(self, species_list) + subroutine gff_register_species(self, species_list) class(GridFire), intent(inout) :: self character(len=*), dimension(:), intent(in) :: species_list @@ -179,7 +276,6 @@ module gridfire_mod 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) @@ -191,17 +287,14 @@ module gridfire_mod 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() + print *, "GridFire: ", self%gff_get_last_error() error stop end if - end subroutine register_species + end subroutine gff_register_species - subroutine setup_policy(self, policy_name, abundances) + subroutine gff_setup_policy(self, policy_name, abundances) class(GridFire), intent(in) :: self character(len=*), intent(in) :: policy_name real(c_double), dimension(:), intent(in) :: abundances @@ -218,41 +311,59 @@ module gridfire_mod self%num_species) if (ierr /= GF_SUCCESS .AND. ierr /= FDSSE_SUCCESS) then - print *, "GridFire Policy Error: ", self%get_last_error() + print *, "GridFire Policy Error: ", self%gff_get_last_error() error stop end if - end subroutine setup_policy + end subroutine gff_setup_policy - subroutine setup_solver(self, solver_name) + subroutine gff_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() + print *, "GridFire Solver Error: ", self%gff_get_last_error() error stop end if - end subroutine setup_solver + end subroutine gff_setup_solver - subroutine evolve(self, Y_in, T, rho, dt, Y_out, energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost, ierr) + subroutine gff_evolve_single(self, Y_in, T, rho, tMax, dt0, Y_out, energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost, ierr) class(GridFire), intent(in) :: self real(c_double), dimension(:), intent(in) :: Y_in - real(c_double), value :: T, rho, dt + real(c_double), intent(in) :: T, rho + real(c_double), value :: tMax, dt0 + real(c_double), dimension(:), intent(out) :: Y_out real(c_double), intent(out) :: energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost integer, intent(out) :: ierr integer(c_int) :: c_ierr - c_ierr = gf_evolve(self%ctx, & + c_ierr = gf_evolve_c_scalar(self%ctx_type, self%ctx, & Y_in, self%num_species, & - T, rho, dt, & + T, rho, tMax, dt0, & Y_out, & energy, dedt, dedrho, nu_e_loss, nu_flux, 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 subroutine gff_evolve_single + + subroutine gff_evolve_multi(self, Y_in, T, rho, tMax, dt0, Y_out, energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost, ierr) + class(GridFire), intent(in) :: self + real(c_double), dimension(:,:), intent(in) :: Y_in + real(c_double), dimension(:), intent(in) :: T, rho + real(c_double), value :: tMax, dt0 + + real(c_double), dimension(:,:), intent(out) :: Y_out + real(c_double), dimension(:), intent(out) :: energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost + integer, intent(out) :: ierr + integer(c_int) :: c_ierr + + c_ierr = gf_evolve_c_array(self%ctx_type, self%ctx, & + Y_in, self%num_species, & + T, rho, tMax, dt0, & + Y_out, & + energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost) + ierr = int(c_ierr) + end subroutine gff_evolve_multi + end module gridfire_mod diff --git a/src/extern/include/gridfire/extern/gridfire_extern.h b/src/extern/include/gridfire/extern/gridfire_extern.h index 6799f7cf..0f58d38e 100644 --- a/src/extern/include/gridfire/extern/gridfire_extern.h +++ b/src/extern/include/gridfire/extern/gridfire_extern.h @@ -7,8 +7,8 @@ extern "C" { #endif enum GF_TYPE { - SINGLE_ZONE = 0, - MULTI_ZONE = 1 + SINGLE_ZONE = 1001, + MULTI_ZONE = 1002 }; diff --git a/src/extern/lib/gridfire_extern.cpp b/src/extern/lib/gridfire_extern.cpp index bec417f4..4bc4be64 100644 --- a/src/extern/lib/gridfire_extern.cpp +++ b/src/extern/lib/gridfire_extern.cpp @@ -222,6 +222,21 @@ extern "C" { void* mass_lost ) { + printf("In C Starting gf_evolve with type %d\n", type); + printf("In C num_species: %zu, tMax: %e, dt0: %e\n", num_species, tMax, dt0); + printf("In C Y_in ptr: %p, T ptr: %p, rho ptr: %p\n", Y_in, T, rho); + // values + printf("In C Y_in first 5 values: "); + const auto* Y_in_ptr = static_cast(Y_in); + for (size_t i = 0; i < std::min(num_species, size_t(5)); ++i) { + printf("%e ", Y_in_ptr[i]); + } + printf("\n"); + printf("In C T value: %e\n", *(static_cast(T))); + printf("In C rho value: %e\n", *(static_cast(rho))); + printf("In C tMax value: %e\n", tMax); + printf("In C dt0 value: %e\n", dt0); + if (!ptr || !Y_in || !T || !rho) { return GF_UNINITIALIZED_INPUT_MEMORY_ERROR; } @@ -252,6 +267,8 @@ extern "C" { auto* specific_neutrino_flux_local = static_cast(specific_neutrino_flux); auto* mass_lost_local = static_cast(mass_lost); + printf("Evolving single zone with T = %e, rho = %e for tMax = %e and dt0 = %e\n", *T_ptr, *rho_ptr, tMax, dt0); + return execute_guarded(ctx, [&]() { return ctx->evolve( static_cast(Y_in), @@ -283,12 +300,7 @@ extern "C" { 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; - // } - // } + printf("Evolving multi zone for tMax = %e and dt0 = %e\n", tMax, dt0); return execute_guarded(ctx, [&]() { return ctx->evolve( diff --git a/tests/extern/C/gridfire_evolve_multi.c b/tests/extern/C/gridfire_evolve_multi.c index 6a414594..bc6609c9 100644 --- a/tests/extern/C/gridfire_evolve_multi.c +++ b/tests/extern/C/gridfire_evolve_multi.c @@ -55,9 +55,11 @@ int main() { double Rhos[ZONES]; for (size_t zone = 0; zone < ZONES; zone++) { - Temps[zone] = 1.0e7; + Temps[zone] = 1.0e7 + (double)zone * 1.0e5; // From 10 million K to 20 million K Rhos[zone] = 1.5e2; + printf("Zone %zu - Temp: %e K, Rho: %e g/cm^3\n", zone, Temps[zone], Rhos[zone]); } + return 0; printf(" Registering species..."); int ret = gf_register_species(ctx, NUM_SPECIES, species_names); diff --git a/tests/extern/C/gridfire_evolve_single.c b/tests/extern/C/gridfire_evolve_single.c index f7efb9a0..8ee076e3 100644 --- a/tests/extern/C/gridfire_evolve_single.c +++ b/tests/extern/C/gridfire_evolve_single.c @@ -50,7 +50,7 @@ int main() { double specific_neutrino_flux; double mass_lost; - const double T_in = 1.5e7; // Temperature in K + const double T_in = 1e7; // Temperature in K const double rho_in = 1.5e2; // Density in g/cm^3 ret = gf_evolve( diff --git a/tests/extern/fortran/gridfire_evolve_multi.f90 b/tests/extern/fortran/gridfire_evolve_multi.f90 new file mode 100644 index 00000000..83381788 --- /dev/null +++ b/tests/extern/fortran/gridfire_evolve_multi.f90 @@ -0,0 +1,134 @@ +program main_multi + use iso_c_binding + use gridfire_mod + implicit none + + ! --- Constants --- + integer, parameter :: NUM_SPECIES = 8 + integer, parameter :: ZONES = 100 + + type(GridFire) :: net + integer(c_int) :: ierr + integer :: i, z + + ! --- 1. Define Species --- + character(len=5), dimension(NUM_SPECIES) :: 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 template + real(c_double), dimension(NUM_SPECIES) :: abundance_root = [ & + 0.702616602672027d0, & + 9.74791583949078d-06, & + 0.06895512307276903d0, & + 0.00025d0, & + 7.855418029399437d-05, & + 0.0006014411598306529d0, & + 8.103062886768109d-05, & + 2.151340851063217d-05 & + ] + + ! --- Multi-Zone Arrays --- + ! Fortran is Column-Major. To match C's Row-Major [ZONES][SPECIES], + ! we dimension as (SPECIES, ZONES) so that Species are contiguous for each Zone. + real(c_double), dimension(NUM_SPECIES, ZONES) :: Y_in, Y_out + + ! Thermodynamic Conditions Arrays + real(c_double), dimension(ZONES) :: T_arr + real(c_double), dimension(ZONES) :: rho_arr + + ! Output Arrays + real(c_double), dimension(ZONES) :: energy_out + real(c_double), dimension(ZONES) :: dedt + real(c_double), dimension(ZONES) :: dedrho + real(c_double), dimension(ZONES) :: snu_e_loss + real(c_double), dimension(ZONES) :: snu_flux + real(c_double), dimension(ZONES) :: dmass + + ! Time settings + real(c_double) :: tMax = 3.0e17 ! 10 Gyr total time + real(c_double) :: dt0 = 1e-12 ! Starting Timestep + + ! --- 2. Setup Data --- + print *, "Testing GridFireEvolve Multi (Fortran)" + print *, " Number of zones: ", ZONES + print *, " Number of species: ", NUM_SPECIES + + do z = 1, ZONES + ! Initialize Abundances (Copy root to every zone) + Y_in(:, z) = abundance_root + + ! Initialize T and Rho gradient + ! T: 1.0e7 -> 2.0e7 in steps of 1.0e5 + T_arr(z) = 1.0d7 + dble(z-1) * 1.0d5 + rho_arr(z) = 1.5d2 + + ! Debug print for first few zones + if (z <= 3) then + print '(A, I0, A, ES12.5, A, ES12.5, A)', & + " Zone ", z-1, " - Temp: ", T_arr(z), " K, Rho: ", rho_arr(z), " g/cm^3" + end if + end do + + ! --- 3. Initialize GridFire Multi-Zone --- + print *, "Initializing GridFire..." + ! Note: Pass integer(c_size_t) for zone count + call net%gff_init(MULTI_ZONE, int(ZONES, c_size_t)) + + ! --- 4. Register Species --- + print *, "Registering species..." + call net%gff_register_species(species_names) + + ! --- 5. Configure Engine & Solver --- + print *, "Setting up Main Sequence Policy..." + call net%gff_setup_policy("MAIN_SEQUENCE_POLICY", abundance_root) + + print *, "Setting up CVODE Solver..." + call net%gff_setup_solver("CVODE") + + ! --- 6. Evolve --- + print *, "Evolving system..." + + ! Note: We pass the arrays T_arr and rho_arr. + ! Ensure your interface change (removing 'value' attribute) is applied + ! so these are passed by reference (address). + call net%gff_evolve(Y_in, T_arr, rho_arr, tMax, dt0, & + Y_out, energy_out, dedt, dedrho, & + snu_e_loss, snu_flux, dmass, ierr) + + if (ierr /= 0) then + print *, "Evolution Failed with error code: ", ierr + print *, "Error Message: ", net%gff_get_last_error() + call net%gff_free() + stop + end if + + ! --- 7. Report Results --- + print *, "" + print *, "--- Results (First 3 Zones) ---" + + do z = 1, 3 + print *, "=== Zone ", z-1, " ===" + print '(A, ES12.5, A)', " Energy Gen: ", energy_out(z), " erg/g/s" + print '(A, ES12.5)', " Mass Lost: ", dmass(z) + print '(A, ES12.5)', " T: ", T_arr(z) + + print *, " Key Abundances (H-1, He-4, C-12):" + print '(3(ES12.5, 1X))', Y_out(1, z), Y_out(3, z), Y_out(4, z) + print *, "" + end do + + print *, "... (Zones 4-", ZONES, " omitted) ..." + + ! --- 8. Cleanup --- + call net%gff_free() + +end program main_multi \ No newline at end of file diff --git a/tests/extern/fortran/gridfire_evolve.f90 b/tests/extern/fortran/gridfire_evolve_single.f90 similarity index 82% rename from tests/extern/fortran/gridfire_evolve.f90 rename to tests/extern/fortran/gridfire_evolve_single.f90 index ac0b9c4d..0f38326a 100644 --- a/tests/extern/fortran/gridfire_evolve.f90 +++ b/tests/extern/fortran/gridfire_evolve_single.f90 @@ -41,30 +41,31 @@ program main ! 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 + real(c_double) :: tMax = 3.0e17 ! 10 Gyr total time + real(c_double) :: dt0 = 1e-12 ! Starting Timestep ! --- 2. Initialize GridFire --- print *, "Initializing GridFire..." - call net%gff_init() + call net%gff_init(SINGLE_ZONE) ! --- 3. Register Species --- print *, "Registering species..." - call net%register_species(species_names) + call net%gff_register_species(species_names) ! --- 4. Configure Engine & Solver --- print *, "Setting up Main Sequence Policy..." - call net%setup_policy("MAIN_SEQUENCE_POLICY", Y_in) + call net%gff_setup_policy("MAIN_SEQUENCE_POLICY", Y_in) print *, "Setting up CVODE Solver..." - call net%setup_solver("CVODE") + call net%gff_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, snu_e_loss, snu_flux, dmass, ierr) + print *, "Evolving system (t = ", tMax, "s dt =", dt0, "s)..." + call net%gff_evolve(Y_in, T, rho, tMax, dt0, Y_out, energy_out, dedt, dedrho, snu_e_loss, snu_flux, dmass, ierr) if (ierr /= 0) then print *, "Evolution Failed with error code: ", ierr - print *, "Error Message: ", net%get_last_error() + print *, "Error Message: ", net%gff_get_last_error() call net%gff_free() ! Always cleanup stop end if diff --git a/tests/extern/fortran/meson.build b/tests/extern/fortran/meson.build index 60601a00..6e7738c3 100644 --- a/tests/extern/fortran/meson.build +++ b/tests/extern/fortran/meson.build @@ -1,5 +1,11 @@ -executable('test_fortran_extern', 'gridfire_evolve.f90', +executable('gf_fortran_single_zone_test', 'gridfire_evolve_single.f90', install: false, fortran_args: ['-Wall', '-Wextra'], dependencies: [gridfire_fortran_dep] -) \ No newline at end of file +) + +executable('gf_fortran_multi_zone_test', 'gridfire_evolve_multi.f90', + install: false, + fortran_args: ['-Wall', '-Wextra'], + dependencies: [gridfire_fortran_dep] +) diff --git a/tools/cli/meson.build b/tools/cli/meson.build deleted file mode 100644 index 336435f3..00000000 --- a/tools/cli/meson.build +++ /dev/null @@ -1 +0,0 @@ -subdir('gf_quick') \ No newline at end of file diff --git a/tools/config/generate_config_files.cpp b/tools/gf_config/generate_config_files.cpp similarity index 100% rename from tools/config/generate_config_files.cpp rename to tools/gf_config/generate_config_files.cpp diff --git a/tools/config/meson.build b/tools/gf_config/meson.build similarity index 100% rename from tools/config/meson.build rename to tools/gf_config/meson.build diff --git a/tools/gf_multi/main.cpp b/tools/gf_multi/main.cpp new file mode 100644 index 00000000..41aa2e51 --- /dev/null +++ b/tools/gf_multi/main.cpp @@ -0,0 +1,393 @@ +// ReSharper disable CppUnusedIncludeDirective +#include +#include +#include +#include +#include + +#include "gridfire/gridfire.h" +#include // Required for parallel_setup + +#include "fourdst/composition/composition.h" +#include "fourdst/logging/logging.h" +#include "fourdst/atomic/species.h" +#include "fourdst/composition/utils.h" + +#include "quill/Logger.h" +#include "quill/Backend.h" +#include "CLI/CLI.hpp" + +#include + +#include "gridfire/utils/gf_omp.h" + +static std::terminate_handler g_previousHandler = nullptr; +static std::vector>>> g_callbackHistory; +static bool s_wrote_abundance_history = false; +void quill_terminate_handler(); + +enum class ScalingTypes { + LINEAR, + LOG, + GEOM +}; + +std::map scaling_type_map = { + {"linear", ScalingTypes::LINEAR}, + {"log", ScalingTypes::LOG}, + {"geom", ScalingTypes::GEOM} +}; + +template +concept IsOrderableLinear = std::floating_point; + +template +concept IsOrderableLog = std::floating_point; + +template +constexpr std::vector linspace(T start, T end, size_t N) { + if (N == 0) { + return {}; + } + if (N == 1) { + return {start}; + } + + std::vector result{}; + result.resize(N); + + for (std::size_t i = 0; i < N; ++i) { + const T t = static_cast(i) / static_cast(N - 1); + result[i] = std::lerp(start, end, t); + } + + return result; +} + +template +std::vector logspace(T start, T end, size_t N) { + if (N == 0) { + return {}; + } + if (N == 1) { + return {start}; + } + + std::vector result{}; + result.resize(N); + + const T log_start = start; + const T log_end = end; + + for (std::size_t i = 0; i < N; ++i) { + const T t = static_cast(i) / static_cast(N - 1); + const T exponent = std::lerp(log_start, log_end, t); + result[i] = std::pow(static_cast(10), exponent); + } + + return result; +} + +template +std::vector geomspace(T start, T end, size_t N) { + if (start <= 0 || end <= 0) { + throw std::domain_error("geomspace requires positive start/end values"); + } + + const T log_start = std::log10(start); + const T log_end = std::log10(end); + + std::vector result = logspace(log_start, log_end, N); + + result[0] = start; + result[N - 1] = end; + + return result; +} + + + +gridfire::NetIn init(const double temp, const double rho, const double tMax) { + std::setlocale(LC_ALL, ""); + g_previousHandler = std::set_terminate(quill_terminate_handler); + quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + logger->set_log_level(quill::LogLevel::Info); + + using namespace gridfire; + const std::vector X = {0.7081145999999999, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; + const std::vector symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"}; + + + const fourdst::composition::Composition composition = fourdst::composition::buildCompositionFromMassFractions(symbols, X); + + NetIn netIn; + netIn.composition = composition; + netIn.temperature = temp; + netIn.density = rho; + netIn.energy = 0; + + netIn.tMax = tMax; + netIn.dt0 = 1e-12; + + return netIn; +} + +void log_results(const gridfire::NetOut& netOut, const gridfire::NetIn& netIn) { + std::vector logSpecies = { + fourdst::atomic::H_1, + fourdst::atomic::He_3, + fourdst::atomic::He_4, + fourdst::atomic::C_12, + fourdst::atomic::N_14, + fourdst::atomic::O_16, + fourdst::atomic::Ne_20, + fourdst::atomic::Mg_24 + }; + + std::vector initial; + std::vector final; + std::vector delta; + std::vector fractional; + for (const auto& species : logSpecies) { + double initial_X = netIn.composition.getMassFraction(species); + double final_X = netOut.composition.getMassFraction(species); + double delta_X = final_X - initial_X; + double fractionalChange = (delta_X) / initial_X * 100.0; + + initial.push_back(initial_X); + final.push_back(final_X); + delta.push_back(delta_X); + fractional.push_back(fractionalChange); + } + + initial.push_back(0.0); // Placeholder for energy + final.push_back(netOut.energy); + delta.push_back(netOut.energy); + fractional.push_back(0.0); // Placeholder for energy + + initial.push_back(0.0); + final.push_back(netOut.dEps_dT); + delta.push_back(netOut.dEps_dT); + fractional.push_back(0.0); + + initial.push_back(0.0); + final.push_back(netOut.dEps_dRho); + delta.push_back(netOut.dEps_dRho); + fractional.push_back(0.0); + + initial.push_back(0.0); + final.push_back(netOut.specific_neutrino_energy_loss); + delta.push_back(netOut.specific_neutrino_energy_loss); + fractional.push_back(0.0); + + initial.push_back(0.0); + final.push_back(netOut.specific_neutrino_flux); + delta.push_back(netOut.specific_neutrino_flux); + fractional.push_back(0.0); + + initial.push_back(netIn.composition.getMeanParticleMass()); + final.push_back(netOut.composition.getMeanParticleMass()); + delta.push_back(final.back() - initial.back()); + fractional.push_back((final.back() - initial.back()) / initial.back() * 100.0); + + std::vector rowLabels = [&]() -> std::vector { + std::vector labels; + for (const auto& species : logSpecies) { + labels.emplace_back(species.name()); + } + labels.emplace_back("ε"); + labels.emplace_back("dε/dT"); + labels.emplace_back("dε/dρ"); + labels.emplace_back("Eν"); + labels.emplace_back("Fν"); + labels.emplace_back("<μ>"); + return labels; + }(); + + + gridfire::utils::Column paramCol("Parameter", rowLabels); + gridfire::utils::Column initialCol("Initial", initial); + gridfire::utils::Column finalCol ("Final", final); + gridfire::utils::Column deltaCol ("δ", delta); + gridfire::utils::Column percentCol("% Change", fractional); + + std::vector> columns; + columns.push_back(std::make_unique>(paramCol)); + columns.push_back(std::make_unique>(initialCol)); + columns.push_back(std::make_unique>(finalCol)); + columns.push_back(std::make_unique>(deltaCol)); + columns.push_back(std::make_unique>(percentCol)); + + + gridfire::utils::print_table("Simulation Results", columns); +} + + +void record_abundance_history_callback(const gridfire::solver::PointSolverTimestepContext& ctx) { + s_wrote_abundance_history = true; + const auto& engine = ctx.engine; + // std::unordered_map> abundances; + std::vector Y; + for (const auto& species : engine.getNetworkSpecies(ctx.state_ctx)) { + const size_t sid = engine.getSpeciesIndex(ctx.state_ctx, species); + double y = N_VGetArrayPointer(ctx.state)[sid]; + Y.push_back(y > 0.0 ? y : 0.0); // Regularize tiny negative abundances to zero + } + + fourdst::composition::Composition comp(engine.getNetworkSpecies(ctx.state_ctx), Y); + + + std::unordered_map> abundances; + for (const auto& sp : comp | std::views::keys) { + abundances.emplace(std::string(sp.name()), std::make_pair(sp.mass(), comp.getMolarAbundance(sp))); + } + g_callbackHistory.emplace_back(ctx.t, abundances); +} + + +void save_callback_data(const std::string_view filename) { + std::set unique_species; + for (const auto &abundances: g_callbackHistory | std::views::values) { + for (const auto &species_name: abundances | std::views::keys) { + unique_species.insert(species_name); + } + } + std::ofstream csvFile(filename.data(), std::ios::out); + csvFile << "t,"; + + size_t i = 0; + for (const auto& species_name : unique_species) { + csvFile << species_name; + if (i < unique_species.size() - 1) { + csvFile << ","; + } + i++; + } + + csvFile << "\n"; + + for (const auto& [time, data] : g_callbackHistory) { + csvFile << time << ","; + size_t j = 0; + for (const auto& species_name : unique_species) { + if (!data.contains(species_name)) { + csvFile << "0.0"; + } else { + csvFile << data.at(species_name).second; + } + if (j < unique_species.size() - 1) { + csvFile << ","; + } + ++j; + } + csvFile << "\n"; + } + + csvFile.close(); +} + +void log_callback_data(const double temp) { + if (s_wrote_abundance_history) { + std::cout << "Saving abundance history to abundance_history.csv" << std::endl; + save_callback_data("abundance_history_" + std::to_string(temp) + ".csv"); + } + +} + +void quill_terminate_handler() +{ + log_callback_data(1.5e7); + quill::Backend::stop(); + if (g_previousHandler) + g_previousHandler(); + else + std::abort(); +} + +void callback_main(const gridfire::solver::PointSolverTimestepContext& ctx) { + record_abundance_history_callback(ctx); +} + +int main(int argc, char** argv) { + GF_PAR_INIT(); + using namespace gridfire; + + double tMin = 1e7; + double tMax = 3e7; + ScalingTypes tempScaling = ScalingTypes::LINEAR; + + double rhoMin = 1e2; + double rhoMax = 1e3; + ScalingTypes rhoScaling = ScalingTypes::LOG; + + double evolveTime = 1e13; + size_t shells = 10; + + + CLI::App app("GridFire Quick CLI Test"); + // Add temp, rho, and tMax as options if desired + app.add_option("--tMin", tMin, "Initial Temperature")->default_val(std::format("{:5.2E}", tMin)); + app.add_option("--tMax", tMax, "Maximum Temperature")->default_val(std::format("{:5.2E}", tMax)); + app.add_option("--tempScaling", tempScaling, "Temperature Scaling Type (linear, log, geom)")->default_val(ScalingTypes::LINEAR)->transform(CLI::CheckedTransformer(scaling_type_map, CLI::ignore_case)); + + app.add_option("--rhoMin", rhoMin, "Initial Density")->default_val(std::format("{:5.2E}", rhoMin)); + app.add_option("--rhoMax", rhoMax, "Maximum Density")->default_val(std::format("{:5.2E}", rhoMax)); + app.add_option("--rhoScaling", rhoScaling, "Density Scaling Type (linear, log, geom)")->default_val(ScalingTypes::LOG)->transform(CLI::CheckedTransformer(scaling_type_map, CLI::ignore_case)); + + app.add_option("--shell", shells, "Number of Shells")->default_val(shells); + + app.add_option("--evolveTime", evolveTime, "Evolution Time")->default_val(std::format("{:5.2E}", evolveTime)); + + + CLI11_PARSE(app, argc, argv); + NetIn rootNetIn = init(tMin, tMax, evolveTime); + + std::vector temps; + std::vector densities; + + switch (tempScaling) { + case ScalingTypes::LINEAR: + temps = linspace(tMin, tMax, static_cast(shells)); + break; + case ScalingTypes::LOG: + temps = logspace(std::log10(tMin), std::log10(tMax), static_cast(shells)); + break; + case ScalingTypes::GEOM: + temps = geomspace(tMin, tMax, static_cast(shells)); + break; + } + switch (rhoScaling) { + case ScalingTypes::LINEAR: + densities = linspace(rhoMin, rhoMax, static_cast(shells)); + break; + case ScalingTypes::LOG: + densities = logspace(std::log10(rhoMin), std::log10(rhoMax), static_cast(shells)); + break; + case ScalingTypes::GEOM: + densities = geomspace(rhoMin, rhoMax, static_cast(shells)); + break; + } + + std::vector netIns; + netIns.reserve(shells); + for (size_t i = 0; i < static_cast(shells); ++i) { + NetIn netIn = rootNetIn; + netIn.temperature = temps[i]; + netIn.density = densities[i]; + netIns.emplace_back(netIn); + } + + policy::MainSequencePolicy stellarPolicy(rootNetIn.composition); + auto [engine, ctx_template] = stellarPolicy.construct(); + + solver::PointSolver point_solver(engine); + solver::GridSolver grid_solver(engine, point_solver); + solver::GridSolverContext solver_ctx(*ctx_template); + + + std::vector results = grid_solver.evaluate(solver_ctx, netIns); + for (size_t i = 0; i < results.size(); ++i) { + std::cout << "=== Shell " << i << " ===" << std::endl; + log_results(results[i], netIns[i]); + } +} \ No newline at end of file diff --git a/tools/gf_multi/meson.build b/tools/gf_multi/meson.build new file mode 100644 index 00000000..66b00b49 --- /dev/null +++ b/tools/gf_multi/meson.build @@ -0,0 +1 @@ +executable('gf_multi', 'main.cpp', dependencies: [gridfire_dep, cli11_dep]) \ No newline at end of file diff --git a/tools/cli/gf_quick/main.cpp b/tools/gf_quick/main.cpp similarity index 98% rename from tools/cli/gf_quick/main.cpp rename to tools/gf_quick/main.cpp index 333f4a44..88a59fb2 100644 --- a/tools/cli/gf_quick/main.cpp +++ b/tools/gf_quick/main.cpp @@ -360,6 +360,10 @@ int main(int argc, char** argv) { CLI11_PARSE(app, argc, argv); NetIn netIn = init(temp, rho, tMax); + for (const auto& [sp, y] : netIn.composition) { + std::println("Species: {}, Abundance: {}", sp.name(), y); + } + return 0; // netIn.composition = rescale(netIn.composition, X, Z); policy::MainSequencePolicy stellarPolicy(netIn.composition); diff --git a/tools/cli/gf_quick/meson.build b/tools/gf_quick/meson.build similarity index 100% rename from tools/cli/gf_quick/meson.build rename to tools/gf_quick/meson.build diff --git a/tools/meson.build b/tools/meson.build index 93d1e5a2..d4fae098 100644 --- a/tools/meson.build +++ b/tools/meson.build @@ -1,4 +1,5 @@ if get_option('build_tools') - subdir('config') - subdir('cli') + subdir('gf_config') + subdir('gf_quick') + subdir('gf_multi') endif \ No newline at end of file