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

@@ -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)')

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")
+        function gf_init(ctx_type) bind(C, name="gf_init")
-            import :: c_ptr
+            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)
+        ! int gf_free(void* gf)
-        subroutine gf_free(gf) bind(C, name="gf_free")
+        function gf_free(ctx_type, ptr) result(c_res) bind(C, name="gf_free")
-            import :: c_ptr
+            import :: c_ptr, c_int
-            type(c_ptr), value :: gf
+            type(c_ptr), value :: ptr
-        end subroutine gf_free
+            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 :: gff_register_species
-        procedure :: setup_policy
+        procedure :: gff_setup_policy
-        procedure :: setup_solver
+        procedure :: gff_setup_solver
-        procedure :: evolve
+        procedure :: gff_get_last_error
-        procedure :: 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
+        end subroutine gff_evolve_single
-                print *, "GridFire Evolve Error: ", self%get_last_error()
+
-            end if
+        subroutine gff_evolve_multi(self, Y_in, T, rho, tMax, dt0, Y_out, energy, dedt, dedrho, nu_e_loss, nu_flux, mass_lost, ierr)
-        end subroutine evolve
+            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

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

@@ -7,8 +7,8 @@
 extern "C" {
 #endif
     enum GF_TYPE {
-        SINGLE_ZONE = 0,
+        SINGLE_ZONE = 1001,
-        MULTI_ZONE = 1
+        MULTI_ZONE = 1002
     };
 
 

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<const double*>(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<const double*>(T)));
+        printf("In C rho value: %e\n", *(static_cast<const double*>(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<double*>(specific_neutrino_flux);
                 auto* mass_lost_local = static_cast<double*>(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<const double*>(Y_in),
@@ -283,12 +300,7 @@ extern "C" {
                 auto* specific_neutrino_flux_local = static_cast<double*>(specific_neutrino_flux);
                 auto* mass_lost_local = static_cast<double*>(mass_lost);
 
-                // for (size_t i = 0; i < ctx->get_zones(); ++i) {
+                printf("Evolving multi zone for tMax = %e and dt0 = %e\n", tMax, dt0);
-                //     if (!Y_out_local[i]) {
-                //         std::cerr << "Uninitialized memory for Y_out at zone " << i << std::endl;
-                //         return GF_UNINITIALIZED_OUTPUT_MEMORY_ERROR;
-                //     }
-                // }
 
                 return execute_guarded(ctx, [&]() {
                     return ctx->evolve(

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);

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(

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

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)..."
+    print *, "Evolving system (t = ", tMax, "s dt =", dt0, "s)..."
-    call net%evolve(Y_in, T, rho, dt, Y_out, energy_out, dedt, dedrho, snu_e_loss, snu_flux, dmass, ierr)
+    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

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]
-)
+)
+
+executable('gf_fortran_multi_zone_test', 'gridfire_evolve_multi.f90',
+           install: false,
+           fortran_args: ['-Wall', '-Wextra'],
+           dependencies: [gridfire_fortran_dep]
+)

tools/cli/meson.build

@@ -1 +0,0 @@
-subdir('gf_quick')

tools/gf_multi/main.cpp

@@ -0,0 +1,393 @@
+// ReSharper disable CppUnusedIncludeDirective
+#include <iostream>
+#include <fstream>
+#include <chrono>
+#include <thread>
+#include <format>
+
+#include "gridfire/gridfire.h"
+#include <cppad/utility/thread_alloc.hpp> // 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 <clocale>
+
+#include "gridfire/utils/gf_omp.h"
+
+static std::terminate_handler g_previousHandler = nullptr;
+static std::vector<std::pair<double, std::unordered_map<std::string, std::pair<double, double>>>> g_callbackHistory;
+static bool s_wrote_abundance_history = false;
+void quill_terminate_handler();
+
+enum class ScalingTypes {
+    LINEAR,
+    LOG,
+    GEOM
+};
+
+std::map<std::string, ScalingTypes> scaling_type_map = {
+    {"linear", ScalingTypes::LINEAR},
+    {"log",    ScalingTypes::LOG},
+    {"geom",   ScalingTypes::GEOM}
+};
+
+template<typename T>
+concept IsOrderableLinear = std::floating_point<T>;
+
+template<typename T>
+concept IsOrderableLog = std::floating_point<T>;
+
+template<IsOrderableLinear T>
+constexpr std::vector<T> linspace(T start, T end, size_t N) {
+    if (N == 0) {
+        return {};
+    }
+    if (N == 1) {
+        return {start};
+    }
+
+    std::vector<T> result{};
+    result.resize(N);
+
+    for (std::size_t i = 0; i < N; ++i) {
+        const T t = static_cast<T>(i) / static_cast<T>(N - 1);
+        result[i] = std::lerp(start, end, t);
+    }
+
+    return result;
+}
+
+template<IsOrderableLog T>
+std::vector<T> logspace(T start, T end, size_t N) {
+    if (N == 0) {
+        return {};
+    }
+    if (N == 1) {
+        return {start};
+    }
+
+    std::vector<T> 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<T>(i) / static_cast<T>(N - 1);
+        const T exponent = std::lerp(log_start, log_end, t);
+        result[i] = std::pow(static_cast<T>(10), exponent);
+    }
+
+    return result;
+}
+
+template<IsOrderableLog T>
+std::vector<T> 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<T> result = logspace<T>(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<double> X            = {0.7081145999999999, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
+    const std::vector<std::string> 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<fourdst::atomic::Species> 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<double> initial;
+    std::vector<double> final;
+    std::vector<double> delta;
+    std::vector<double> 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<std::string> rowLabels = [&]() -> std::vector<std::string> {
+        std::vector<std::string> 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<std::string> paramCol("Parameter", rowLabels);
+    gridfire::utils::Column<double> initialCol("Initial", initial);
+    gridfire::utils::Column<double> finalCol  ("Final", final);
+    gridfire::utils::Column<double> deltaCol  ("δ", delta);
+    gridfire::utils::Column<double> percentCol("% Change", fractional);
+
+    std::vector<std::unique_ptr<gridfire::utils::ColumnBase>> columns;
+    columns.push_back(std::make_unique<gridfire::utils::Column<std::string>>(paramCol));
+    columns.push_back(std::make_unique<gridfire::utils::Column<double>>(initialCol));
+    columns.push_back(std::make_unique<gridfire::utils::Column<double>>(finalCol));
+    columns.push_back(std::make_unique<gridfire::utils::Column<double>>(deltaCol));
+    columns.push_back(std::make_unique<gridfire::utils::Column<double>>(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<std::string, std::pair<double, double>> abundances;
+    std::vector<double> 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<std::string, std::pair<double, double>> 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<std::string> 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<double> temps;
+    std::vector<double> densities;
+
+    switch (tempScaling) {
+        case ScalingTypes::LINEAR:
+            temps = linspace<double>(tMin, tMax, static_cast<size_t>(shells));
+            break;
+        case ScalingTypes::LOG:
+            temps = logspace<double>(std::log10(tMin), std::log10(tMax), static_cast<size_t>(shells));
+            break;
+        case ScalingTypes::GEOM:
+            temps = geomspace<double>(tMin, tMax, static_cast<size_t>(shells));
+            break;
+    }
+    switch (rhoScaling) {
+        case ScalingTypes::LINEAR:
+            densities = linspace<double>(rhoMin, rhoMax, static_cast<size_t>(shells));
+            break;
+        case ScalingTypes::LOG:
+            densities = logspace<double>(std::log10(rhoMin), std::log10(rhoMax), static_cast<size_t>(shells));
+            break;
+        case ScalingTypes::GEOM:
+            densities = geomspace<double>(rhoMin, rhoMax, static_cast<size_t>(shells));
+            break;
+    }
+
+    std::vector<NetIn> netIns;
+    netIns.reserve(shells);
+    for (size_t i = 0; i < static_cast<size_t>(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<NetOut> 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]);
+    }
+}

tools/gf_multi/meson.build

@@ -0,0 +1 @@
+executable('gf_multi', 'main.cpp', dependencies: [gridfire_dep, cli11_dep])

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);

tools/meson.build

@@ -1,4 +1,5 @@
 if get_option('build_tools')
-    subdir('config')
+    subdir('gf_config')
-    subdir('cli')
+    subdir('gf_quick')
+    subdir('gf_multi')
 endif