docs(readme): updated readme to v0.7.0-alpha

Doxyfile

@@ -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         = 0.6.0
+PROJECT_NUMBER         = v0.7.0-alpha
 
 # 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

README.md

@@ -24,6 +24,7 @@
     - [pypi](#pypi)
     - [source](#source)
     - [source for developers](#source-for-developers)
+    - [patching shared object files](#patching-shared-object-files)
   - [Automatic Build and Installation](#automatic-build-and-installation)
     - [Script Build and Installation Instructions](#script-build-and-installation-instructions)
     - [Currently, known good platforms](#currently-known-good-platforms)
@@ -62,6 +63,16 @@
     - [Python callbacks](#python-callbacks)
 - [Related Projects](#related-projects)
 
+# Version and Notes
+This repository is currently tracking GridFire version v0.7.0-alpha. Note that this
+is a development version which has known scientific and build issues:
+
+1. Over long runs (say 10Gyr) repartitioning stages can introduce discontinuities into abundances
+2. We do not currently produce He-4 at a rate consistent with literature values. This is a known issue and is being addressed.
+3. When using Weak Rate Library (WRL) weak reactions the network becomes pathologically stiff. Reaclib includes a limited set of reactions which can be used to close the CNO cycle. Network construction then defaults to using all reaclib reactions while we address pathological stiffness with WRL rates.
+4. WRL reactions do track energy loss due to neutrinos and neutrino flux; however, these are not currently reported to the user. They will be in the final v0.7.0 release.
+5. There is a current bug in meson-python which results in multiple duplicate LC_RPATH entries in any shared object files compiled for python linking. Recent versions of the macos dynamic loader (XCode command line tools versions >= 16) refuse to load shared object files with duplicat rpath entries. Because of this running `pip install`. in the root will result in a broken gridfire python install. Instead we have bundled a helper script `pip_install_mac_patch.sh` which should be used to install python bindings on macos for the time being. 
+
 # Introduction
 GridFire is a C++ library designed to perform general nuclear network
 evolution. It is part of the larger SERiF project within the 4D-STAR
@@ -109,7 +120,8 @@ are just that and should maintain nearly the same speed as the C++ code. End
 users are strongly encouraged to use the python module rather than the C++ code.
 
 ### pypi
-Installing from pip is as simple as
+Installing from pip is as simple as. Note that this will install gridfire v0.5.0, currently the latest version on pip. Once
+v0.7.0 is released this will be pushed to pip.
 ```bash
 pip install gridfire
 ```
@@ -150,6 +162,8 @@ pip install .
 > **Note:** that if you do not have all system dependencies installed this will
 > fail, the steps in further sections address these in more detail.
 
+> **Note:** If you are using macos you should use the included `pip_install_mac_patch.sh` script instead of `pip install .` as this will automatically patch the build shared object libraries such that they can be loaded by the macos dynamic loader.
+
 ### source for developers
 If you are a developer and would like an editable and incremental python
 install `meson-python` makes this very easy
@@ -166,6 +180,29 @@ of these it does still take a few seconds to recompile regardless of how small
 a source code change you have made). It is **strongly** recommended that
 developers use this approach and end users *do not*.
 
+#### Patching Shared Object Files
+If you need to patch shared object files generated by meson-python directly you should first located the shared object file
+these will be in the site-packages and site-packages/fourdst directories for your python enviroment.
+
+Look for files named
+
+- `site-packages/gridfire.cpython-3*-darwin.so`
+- `site-packages/fourdst/_phys.cpython-3*-darwin.so` 
+
+then, for each of these files, run
+
+```bash
+otool -l <Path/to/file> | grep RPATH -A2
+```
+
+count the number of occurences of duplicate RPATH entries (these should look like `@loaderpath/.gridfire.mesonpy.libs` or `@loaderpath/../.fourdst.mesonpy.libs`). Then use `install_name_tool` to remove **all but one of these** from each shared object file.
+
+If for example there are 4 occurrences of the path `@loader_path/../.fourdst.mesonpy.libs` in `_phys.cpython-3*-darwin.so` then you should run the following command 3 times
+```bash
+install_name_tool -delete_rpath @loader_path/../.fourdst.mesonpy.libs site-packages/fourdst/_phys.cpython-314-darwin.so
+```
+
+the same for the other shared object file (make sure to count the duplicate rpath entries for each separately as there may be a different number of duplicates in each shared object file).
 
 ## Automatic Build and Installation
 ### Script Build and Installation Instructions
@@ -351,8 +388,10 @@ include:
   2000](https://www.sciencedirect.com/science/article/pii/S0092640X00908349?via%3Dihub]))
   to weight reaction rates based on nuclear properties.
 - **Solver Module:** Defines numerical integration strategies (e.g.,
-  `DirectNetworkSolver`) for solving the stiff ODE systems arising from reaction
+  `CVODESolverStrategy`) for solving the stiff ODE systems arising from reaction
   networks.
+- **io Module:** Defines shared interface for parsing network data from files
+- **trigger Module:** Defines interface for complex trigger logic so that repartitioning can be followed.
 - **Python Interface:** Exposes *almost* all C++ functionality to Python,
   allowing users to define compositions, configure engines, and run simulations
   directly from Python scripts.
@@ -404,6 +443,7 @@ construction and rate evaluations:
       member of the `NetworkBuildDepth` enum or an integer.
     - `partition::PartitionFunction`: Partition function used when evaluating
       detailed balance for inverse rates.
+    - `NetworkConstructionFlags`: A bitwise flag telling the network how to construct itself. That is, what reaction types should be used in construction. For example one might use `NetworkConstructionFlags::STRONG | NetworkConstructionFlags::BETA_PLUS` to use all strong reactions and β+ decay. By Default this is set to use reaclib strong and reaclib weak (no WRL included by default due to current pathological stiffness issues).
 
 - **setPrecomputation(bool precompute):**
     - Enable/disable caching of reaction rates and stoichiometric data at initialization.
@@ -533,48 +573,23 @@ A `NetOut` struct contains
 >**Note:** Currently `GraphEngine` only considers energy due to nuclear mass
 >defect and not neutrino loss.
 
+### CVODESolverStrategy
 
-### DirectNetworkSolver (Implicit Rosenbrock Method)
+We use the CVODE module from [SUNDIALS](https://computing.llnl.gov/projects/sundials/cvode) as our primary numerical
+solver. Specifically we use the BDF linear multistep method from that which includes advanced adaptive timestepping.
 
-- **Integrator:** Implicit Rosenbrock4 scheme (order 4) via `Boost.Odeint`’s
+Further, we use a trigger system to periodically repartition the network as the state of the network changes. This
-  `rosenbrock4<double>`, optimized for stiff reaction networks with adaptive step
+keeps the stiffness of the network tractable. The algorithm we use for that is 
-  size control using configurable absolute and relative tolerances.
-- **Jacobian Assembly:** Asks the base engine for the Jacobian Matrix
-- **RHS Evaluation:** Asks the base engine for RHS of the abundance evolution
-  equations
-- **Linear Algebra:** Utilizes `Boost.uBLAS` for state vectors and dense Jacobian
-  matrices, with sparse access patterns supported via coordinate lists of nonzero
-  entries.
-- **Error Control and Logging:** Absolute and relative tolerance parameters
-  (`absTol`, `relTol`) are read from configuration; Quill loggers, which run in a
-  separate non blocking thread, capture integration diagnostics and step
-  statistics.
 
-### Algorithmic Workflow in DirectNetworkSolver
+1. Trigger every 1000th time that the simulation time exceeds the simulationTimeInterval
-1. **Initialization:** Convert input temperature to T9 units, retrieve
+2. OR if any off-diagonal Jacobian entry exceeds the offDiagonalThreshold
-   tolerances, and initialize state vector `Y` from equilibrated composition.
+3. OR every 10th time that the timestep growth exceeds the timestepGrowthThreshold (relative or absolute)
-2. **Integrator Setup:** Construct the controlled Rosenbrock4 stepper and bind
+4. OR if the number of convergence failures grows more than 100% from one step to the next or exceeds 5 at any given step.
-   `RHSManager` and `JacobianFunctor`.
-3. **Adaptive Integration Loop:**
-    - Perform `integrate_adaptive` advancing until `tMax`, catching any
-      `StaleEngineTrigger` to repartition the network and update composition.
-    - On each substep, observe states and log via `RHSManager::observe`.
-4. **Finalization:** Assemble final mass fractions, compute accumulated energy,
-   and populate `NetOut` with updated composition and diagnostics.
 
-### Future Solver Implementations
+Moreover, callback functions can be registered in either python or C++ which will take a `const CVODESolverStrategy::TimestepContext&` struct 
-- **Operator Splitting Solvers:** Strategies to decouple thermodynamics,
+as argument. This allows for more complex logging logic. Note that callbacks **do not** let you reach inside the 
-  screening, and reaction substeps for performance on stiff, multiscale
+solver and adjust the state of the network. They are only intended for investigation not extension of physics. If you
-  networks.
+wish to extend the physics this must be implemented at the engine or engine view level.
-- **GPU-Accelerated Solvers:** Planned use of CUDA/OpenCL backends for
-  large-scale network integration.
-- **Callback observer support:** Currently we use an observer built into our
-  `RHSManager` (`RHSManager::observe`); however, we intend to include support for
-  custom, user defined, observer method.
-
-These strategies can be developed by inheriting from `NetworkSolverStrategy`
-and registering against the same engine types without modifying existing engine
-code.
 
 ## Python Extensibility
 Through the Python bindings, users can subclass engine view classes directly in
@@ -667,7 +682,7 @@ int main(){
     gridfire::AdaptiveEngineView adaptView(msView);
 
     // 5. Construct implicit solver (handles remaining stiffness)
-    gridfire::DirectNetworkSolver solver(adaptView);
+    gridfire::CVODESolverStrategey solver(adaptView);
 
     // 6. Prepare input conditions
     NetIn input{
@@ -684,20 +699,6 @@ int main(){
 }
 ```
 
-#### Workflow Components and Effects
-- **GraphEngine** constructs the full reaction network, capturing all species
-  and reactions.
-- **MultiscalePartitioningEngineView** segregates reactions by characteristic
-  timescales (Hix & Thielemann), reducing the effective stiffness by treating
-  fast processes separately.
-- **AdaptiveEngineView** prunes low-flux species/reactions at runtime,
-  decreasing dimensionality and improving computational efficiency.
-- **DirectNetworkSolver** employs an implicit Rosenbrock method to stably
-  integrate the remaining stiff system with adaptive step control.
-
-This layered approach enhances stability for stiff networks while maintaining
-accuracy and performance.
-
 ### Callback Example
 Custom callback functions can be registered with any solver. Because it might make sense for each solver to provide
 different context to the callback function, you should use the struct `gridfire::solver::<SolverName>::TimestepContext`
@@ -712,7 +713,7 @@ the callback function.
 
 #include <iostream>
 
-void callback(const gridfire::solver::DirectNetworkSolver::TimestepContext& context) {
+void callback(const gridfire::solver::CVODESolverStrategy::TimestepContext& context) {
     int H1Index = context.engine.getSpeciesIndex(fourdst::atomic::H_1);
     int He4Index = context.engine.getSpeciesIndex(fourdst::atomic::He_4);
 
@@ -741,7 +742,7 @@ int main(){
     gridfire::AdaptiveEngineView adaptView(msView);
 
     // 5. Construct implicit solver (handles remaining stiffness)
-    gridfire::DirectNetworkSolver solver(adaptView);
+    gridfire::CVODESolverStrategy solver(adaptView);
     solver.set_callback(callback);
 
     // 6. Prepare input conditions
@@ -758,12 +759,29 @@ int main(){
     std::cout << "Final results are: " << output << std::endl;
 }
 ```
+>**Note:** If you want to see exactly why each repartitioning stage was triggered in a human readable manner add the flag True to `solver.evaluate` (`solver.evaluate(input, true)`).
 
 >**Note:** A fully detailed list of all available information in the TimestepContext struct is available in the API documentation.
 
 >**Note:** The order of species in the boost state vector (`ctx.state`) is **not guaranteed** to be any particular order run over run. Therefore, in order to reliably extract
 > values from it, you **must** use the `getSpeciesIndex` method of the engine to get the index of the species you are interested in (these will always be in the same order).
 
+If you wish to know what is provided by a solver context without investigating the code you can simply do
+
+```c++
+void callback(const gridfire::solver::SolverContextBase& context) {
+    for (const auto& [parameterName, description] : context.describe()) {
+        std::cout << parameterName << ": " << description << "\n";
+    }
+    std::cout << std::flush();
+    exit(0);
+}
+```
+
+If you set this as the callback (to any solver strategy) it will print out the available parameters and what they
+are and then close the code. This is useful when writing new callbacks.
+
+
 #### Callback Context
 
 Since each solver may provide different context to the callback function, and it may be frustrating to refer to the
@@ -789,7 +807,7 @@ All GridFire C++ types have been bound and can be passed around as one would exp
 This example implements the same logic as the above C++ example
 ```python
 from gridfire.engine import GraphEngine, MultiscalePartitioningEngineView, AdaptiveEngineView
-from gridfire.solver import DirectNetworkSolver
+from gridfire.solver import CVODESolverStrategey
 from gridfire.type import NetIn
 
 from fourdst.composition import Composition
@@ -819,7 +837,7 @@ qseEngine = MultiscalePartitioningEngineView(baseEngine)
 
 adaptiveEngine = AdaptiveEngineView(qseEngine)
 
-solver = DirectNetworkSolver(adaptiveEngine)
+solver = CVODESolverStrategey(adaptiveEngine)
 
 results = solver.evaluate(netIn)
 
@@ -869,12 +887,14 @@ adaptiveEngine = AdaptiveEngineView(qseEngine)
 solver = DirectNetworkSolver(adaptiveEngine)
 
 
+data: List[Tuple[float, Dict[str, Tuple[float, float]]]] = []
 def callback(context):
-    H1Index = context.engine.getSpeciesIndex(species["H-1"])
+    engine = context.engine
-    He4Index = context.engine.getSpeciesIndex(species["He-4"])
+    abundances: Dict[str, Tuple[float, float]] = {}
-    C12ndex = context.engine.getSpeciesIndex(species["C-12"])
+    for species in engine.getNetworkSpecies():
-    Mgh24ndex = context.engine.getSpeciesIndex(species["Mg-24"])
+        sid = engine.getSpeciesIndex(species)
-    print(f"Time: {context.t}, H-1: {context.state[H1Index]}, He-4: {context.state[He4Index]}, C-12: {context.state[C12ndex]}, Mg-24: {context.state[Mgh24ndex]}")
+        abundances[species.name()] = (species.mass(), context.state[sid])
+    data.append((context.t,abundances))
 
 solver.set_callback(callback)
 results = solver.evaluate(netIn)

meson.build

@@ -18,7 +18,7 @@
 #   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 #
 # *********************************************************************** #
-project('GridFire', 'cpp', version: '0.6.0', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
+project('GridFire', 'cpp', version: 'v0.7.0-alpha', 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')