From e87206a4a3812f31855df0dd6ca45ec3873aad75 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 4 Nov 2025 14:03:46 -0500 Subject: [PATCH] docs(readme): updated readme to v0.7.0-alpha --- Doxyfile | 2 +- README.md | 148 +++++++++++++++++++++++++++++----------------------- meson.build | 2 +- 3 files changed, 86 insertions(+), 66 deletions(-) diff --git a/Doxyfile b/Doxyfile index d0b66cfe..94b0a04c 100644 --- a/Doxyfile +++ b/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 diff --git a/README.md b/README.md index 542f5b0e..2fd68fdd 100644 --- a/README.md +++ b/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 | 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 - `rosenbrock4`, optimized for stiff reaction networks with adaptive step - 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. +Further, we use a trigger system to periodically repartition the network as the state of the network changes. This +keeps the stiffness of the network tractable. The algorithm we use for that is -### Algorithmic Workflow in DirectNetworkSolver -1. **Initialization:** Convert input temperature to T9 units, retrieve - tolerances, and initialize state vector `Y` from equilibrated composition. -2. **Integrator Setup:** Construct the controlled Rosenbrock4 stepper and bind - `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. +1. Trigger every 1000th time that the simulation time exceeds the simulationTimeInterval +2. OR if any off-diagonal Jacobian entry exceeds the offDiagonalThreshold +3. OR every 10th time that the timestep growth exceeds the timestepGrowthThreshold (relative or absolute) +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. -### Future Solver Implementations -- **Operator Splitting Solvers:** Strategies to decouple thermodynamics, - screening, and reaction substeps for performance on stiff, multiscale - networks. -- **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. +Moreover, callback functions can be registered in either python or C++ which will take a `const CVODESolverStrategy::TimestepContext&` struct +as argument. This allows for more complex logging logic. Note that callbacks **do not** let you reach inside the +solver and adjust the state of the network. They are only intended for investigation not extension of physics. If you +wish to extend the physics this must be implemented at the engine or engine view level. ## 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::::TimestepContext` @@ -712,7 +713,7 @@ the callback function. #include -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"]) - He4Index = context.engine.getSpeciesIndex(species["He-4"]) - C12ndex = context.engine.getSpeciesIndex(species["C-12"]) - Mgh24ndex = context.engine.getSpeciesIndex(species["Mg-24"]) - 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]}") + 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) results = solver.evaluate(netIn) diff --git a/meson.build b/meson.build index ba0ed33d..80ca4e02 100644 --- a/meson.build +++ b/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')