Merge pull request #10 from tboudreaux/perf/openMP

Performance + openMP support

.gitignore

@@ -82,6 +82,7 @@ subprojects/cvode-*/
 subprojects/kinsol-*/
 subprojects/CLI11-*/
 subprojects/openssl-*/
+subprojects/tomlplusplus-*/
 *.fbundle
 *.wraplock
 
@@ -98,6 +99,8 @@ liblogging.wrap
 libplugin.wrap
 minizip-ng.wrap
 openssl.wrap
+glaze.wrap
+tomlplusplus.wrap
 
 .vscode/
 

build-check/CPPC/meson.build

@@ -0,0 +1,37 @@
+cppc = meson.get_compiler('cpp')
+
+if cppc.get_id() == 'clang'
+    message('disabling bitwise-instead-of-logical warnings for clang')
+    add_project_arguments('-Wno-bitwise-instead-of-logical', language: 'cpp')
+endif
+
+if cppc.get_id() == 'gcc'
+    message('disabling psabi warnings for gcc')
+    add_project_arguments('-Wno-psabi', language: 'cpp')
+    
+    if (cppc.version().version_compare('<14.0'))
+      error('g++ version must be at least 14.0, found ' + cppc.version())
+    endif
+endif
+
+if not cppc.has_header('print')
+    error('C++ standard library header <print> not found. Please ensure your compiler and standard library supports C++23. We have already validated your compiler version so this is likely an issue with your standard library installation.')
+endif
+if not cppc.has_header('format')
+    error('C++ standard library header <format> not found. Please ensure your compiler and standard library supports C++23. We have already validated your compiler version so this is likely an issue with your standard library installation.')
+endif
+
+# For Eigen
+add_project_arguments('-Wno-deprecated-declarations', language: 'cpp')
+
+if get_option('build_python')
+    message('enabling hidden visibility for C++ symbols when building Python extension. This reduces the size of the resulting shared library.')
+    add_project_arguments('-fvisibility=hidden', language: 'cpp')
+else
+    message('enabling default visibility for C++ symbols')
+    add_project_arguments('-fvisibility=default', language: 'cpp')
+endif
+
+if get_option('openmp_support')
+    add_project_arguments('-DGRIDFIRE_USE_OPENMP', language: 'cpp')
+endif

build-check/FC/meson.build

@@ -0,0 +1,15 @@
+if get_option('build_fortran')
+  add_languages('fortran', native: true)
+  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)')
+  fc = meson.get_compiler('fortran')
+  if not get_option('unsafe_fortran')
+    if fc.get_id() != 'gcc'
+      error('The only supported fortran compiler for GridFire is gfortran (version >= 14.0), found ' + fc  + '. GridFire has not been tested with any other compilers. You can disable this check with the -Dunsafe-fortran=true flag to try other compilers')
+    endif
+  endif
+  if (fc.version().version_compare('<14.0'))
+      error('gfortran version must be at least 14.0, found ' + fc.version())
+  endif
+endif

build-check/meson.build

@@ -0,0 +1,15 @@
+message('Found CXX compiler: ' + meson.get_compiler('cpp').get_id())
+message('C++ standard set to: ' + get_option('cpp_std'))
+
+cc = meson.get_compiler('c')
+
+ignore_unused_args = '-Wno-unused-command-line-argument'
+
+add_global_arguments(ignore_unused_args, language: 'cpp')
+add_global_arguments(ignore_unused_args, language: 'c')
+
+
+subdir('CPPC')
+subdir('FC')
+
+

build-config/fourdst/meson.build

@@ -1,17 +1,21 @@
 # bring in all of the fourdst utility repositories
 
 fourdst_build_lib_all = true
-if get_option('unity-safe')
+if not get_option('plugin_support')
     fourdst_build_lib_all=false
+    message('Disabling fourdst plugin support as per user request.')
 endif
 
 fourdst_sp = subproject('fourdst',
                         default_options:
-                            ['build-tests=' + get_option('build-tests').to_string(),
-                             'build-python=' + get_option('build-python').to_string(),
-                             'build-lib-all=' + fourdst_build_lib_all.to_string(),
-                             'pkg-config=' + get_option('pkg-config').to_string(),
-                             'build-lib-log=true'
+                            ['build_tests=' + get_option('build_tests').to_string(),
+                             'build_python=' + get_option('build_python').to_string(),
+                             'build_lib_all=' + fourdst_build_lib_all.to_string(),
+                             'build_lib_comp=true',
+                             'build_lib_config=true',
+                             'build_lib_log=true',
+                             'build_lib_const=true',
+                             'pkg_config=' + get_option('pkg_config').to_string(),
                             ]
 )
 
@@ -19,15 +23,16 @@ composition_dep = fourdst_sp.get_variable('composition_dep')
 log_dep = fourdst_sp.get_variable('log_dep')
 const_dep = fourdst_sp.get_variable('const_dep')
 config_dep = fourdst_sp.get_variable('config_dep')
-if not get_option('unity-safe')
+if get_option('plugin_support')
+    warning('Including plugin library from fourdst. Note this will bring in minizip-ng and openssl, which can cause build issues with cross compilation due to their complexity.')
     plugin_dep = fourdst_sp.get_variable('plugin_dep')
 endif
 
 libcomposition = fourdst_sp.get_variable('libcomposition')
 libconst = fourdst_sp.get_variable('libconst')
-libconfig = fourdst_sp.get_variable('libconfig')
 liblogging = fourdst_sp.get_variable('liblogging')
 
-if not get_option('unity-safe')
+if get_option('plugin_support')
+    warning('Including plugin library from fourdst. Note this will bring in minizip-ng and openssl, which can cause build issues with cross compilation due to their complexity.')
     libplugin = fourdst_sp.get_variable('libplugin')
 endif

build-config/meson.build

@@ -1,7 +1,9 @@
 cmake = import('cmake')
 
-
-subdir('python')
+if get_option('build_python')
+  subdir('python')
+  subdir('pybind')
+endif
 subdir('fourdst')
 subdir('sundials')
 
@@ -11,6 +13,5 @@ subdir('eigen')
 
 subdir('json')
 
-subdir('pybind')
 
 subdir('CLI11')

build-config/sundials/cvode/meson.build

@@ -6,7 +6,7 @@ cvode_cmake_options.add_cmake_defines({
     'CMAKE_C_FLAGS'     : '-Wno-deprecated-declarations',
     'BUILD_SHARED_LIBS' : 'OFF',
     'BUILD_STATIC_LIBS' : 'ON',
-    'EXAMPLES_ENABLE_C': 'OFF',
+    'EXAMPLES_ENABLE_C' : 'OFF',
     'CMAKE_POSITION_INDEPENDENT_CODE': true
 
 })
@@ -16,6 +16,15 @@ cvode_cmake_options.add_cmake_defines({
     'CMAKE_INSTALL_INCLUDEDIR': get_option('includedir')
 })
 
+if meson.is_cross_build() and host_machine.system() == 'emscripten'
+  cvode_cmake_options.add_cmake_defines({
+      'CMAKE_C_FLAGS': '-s MEMORY64=1 -s ALLOW_MEMORY_GROWTH=1',
+      'CMAKE_CXX_FLAGS': '-s MEMORY64=1 -s ALLOW_MEMORY_GROWTH=1',
+      'CMAKE_SHARED_LINKER_FLAGS': '-s MEMORY64=1 -s ALLOW_MEMORY_GROWTH=1',
+      'CMAKE_EXE_LINKER_FLAGS': '-s MEMORY64=1 -s ALLOW_MEMORY_GROWTH=1'
+  })
+endif
+
 cvode_sp = cmake.subproject(
     'cvode',
     options: cvode_cmake_options,

build-extra/log-level/meson.build

@@ -0,0 +1,32 @@
+llevel = get_option('log_level')
+
+logbase='QUILL_COMPILE_ACTIVE_LOG_LEVEL_'
+
+if (llevel == 'traceL3')
+    message('Setting log level to TRACE_L3')
+    log_argument = logbase + 'TRACE_L3'
+elif (llevel == 'traceL2')
+    message('Setting log level to TRACE_L2')
+    log_argument = logbase + 'TRACE_L2'
+elif (llevel == 'traceL1')
+    message('Setting log level to TRACE_L1')
+    log_argument = logbase + 'TRACE_L1'
+elif (llevel == 'debug')
+    message('Setting log level to DEBUG')
+    log_argument = logbase + 'DEBUG'
+elif (llevel == 'info')
+    message('Setting log level to INFO')
+    log_argument = logbase + 'INFO'
+elif (llevel == 'warning')
+    message('Setting log level to WARNING')
+    log_argument = logbase + 'WARNING'
+elif (llevel == 'error')
+    message('Setting log level to ERROR')
+    log_argument = logbase + 'ERROR'
+elif (llevel == 'critical')
+    message('Setting log level to CRITICAL')
+    log_argument = logbase + 'CRITICAL'
+endif
+
+log_argument = '-DQUILL_COMPILE_ACTIVE_LOG_LEVEL=' + log_argument
+add_project_arguments(log_argument, language: 'cpp')

build-extra/pkg-config/meson.build

@@ -0,0 +1,18 @@
+if get_option('pkg_config')
+    message('Generating pkg-config file for GridFire...')
+    pkg = import('pkgconfig')
+    pkg.generate(
+        name: 'gridfire',
+        description: 'GridFire nuclear reaction network solver',
+        version: meson.project_version(),
+        libraries: [
+            libgridfire,
+            libcomposition,
+            libconst,
+            liblogging
+        ],
+        subdirs: ['gridfire'],
+        filebase: 'gridfire',
+        install_dir: join_paths(get_option('libdir'), 'pkgconfig')
+    )
+endif

build-python/meson.build

@@ -1,90 +1,95 @@
+if get_option('build_python')
+  message('Building Python bindings...')
 
-gridfire_py_deps = [
-    py_dep,
-    pybind11_dep,
-    const_dep,
-    config_dep,
-    composition_dep,
-    gridfire_dep
-]
+  gridfire_py_deps = [
+      py_dep,
+      pybind11_dep,
+      const_dep,
+      config_dep,
+      composition_dep,
+      gridfire_dep
+  ]
 
-py_sources =  [
-        meson.project_source_root() + '/src/python/bindings.cpp',
-        meson.project_source_root() + '/src/python/types/bindings.cpp',
-        meson.project_source_root() + '/src/python/partition/bindings.cpp',
-        meson.project_source_root() + '/src/python/partition/trampoline/py_partition.cpp',
-        meson.project_source_root() + '/src/python/reaction/bindings.cpp',
-        meson.project_source_root() + '/src/python/screening/bindings.cpp',
-        meson.project_source_root() + '/src/python/screening/trampoline/py_screening.cpp',
-        meson.project_source_root() + '/src/python/io/bindings.cpp',
-        meson.project_source_root() + '/src/python/io/trampoline/py_io.cpp',
-        meson.project_source_root() + '/src/python/exceptions/bindings.cpp',
-        meson.project_source_root() + '/src/python/engine/bindings.cpp',
-        meson.project_source_root() + '/src/python/engine/trampoline/py_engine.cpp',
-        meson.project_source_root() + '/src/python/solver/bindings.cpp',
-        meson.project_source_root() + '/src/python/solver/trampoline/py_solver.cpp',
-        meson.project_source_root() + '/src/python/policy/bindings.cpp',
-        meson.project_source_root() + '/src/python/policy/trampoline/py_policy.cpp',
-        meson.project_source_root() + '/src/python/utils/bindings.cpp',
-    ]
+  py_sources =  [
+          meson.project_source_root() + '/src/python/bindings.cpp',
+          meson.project_source_root() + '/src/python/types/bindings.cpp',
+          meson.project_source_root() + '/src/python/partition/bindings.cpp',
+          meson.project_source_root() + '/src/python/partition/trampoline/py_partition.cpp',
+          meson.project_source_root() + '/src/python/reaction/bindings.cpp',
+          meson.project_source_root() + '/src/python/screening/bindings.cpp',
+          meson.project_source_root() + '/src/python/screening/trampoline/py_screening.cpp',
+          meson.project_source_root() + '/src/python/io/bindings.cpp',
+          meson.project_source_root() + '/src/python/io/trampoline/py_io.cpp',
+          meson.project_source_root() + '/src/python/exceptions/bindings.cpp',
+          meson.project_source_root() + '/src/python/engine/bindings.cpp',
+          meson.project_source_root() + '/src/python/engine/trampoline/py_engine.cpp',
+          meson.project_source_root() + '/src/python/solver/bindings.cpp',
+          meson.project_source_root() + '/src/python/solver/trampoline/py_solver.cpp',
+          meson.project_source_root() + '/src/python/policy/bindings.cpp',
+          meson.project_source_root() + '/src/python/policy/trampoline/py_policy.cpp',
+          meson.project_source_root() + '/src/python/utils/bindings.cpp',
+      ]
 
 
-if meson.is_cross_build() and host_machine.system() == 'darwin'
-  py_mod = shared_module(
-    '_gridfire',
-    sources: py_sources,
-    dependencies: gridfire_py_deps,
-    name_prefix: '',
-    name_suffix: 'so',
-    install: true,
-    install_dir: py_installation.get_install_dir() + '/gridfire'
-    )
-else
-  py_mod = py_installation.extension_module(
-      '_gridfire', # Name of the generated .so/.pyd file (without extension)
+  if meson.is_cross_build() and host_machine.system() == 'darwin'
+    py_mod = shared_module(
+      '_gridfire',
       sources: py_sources,
-      dependencies : gridfire_py_deps,
-      install : true,
+      dependencies: gridfire_py_deps,
+      name_prefix: '',
+      name_suffix: 'so',
+      install: true,
+      install_dir: py_installation.get_install_dir() + '/gridfire'
+      )
+  else
+    py_mod = py_installation.extension_module(
+        '_gridfire', # Name of the generated .so/.pyd file (without extension)
+        sources: py_sources,
+        dependencies : gridfire_py_deps,
+        install : true,
+        subdir: 'gridfire',
+    )
+  endif
+
+
+  py_installation.install_sources(
+      files(
+          meson.project_source_root() + '/src/python/gridfire/__init__.py',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/__init__.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/exceptions.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/partition.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/reaction.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/screening.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/io.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/solver.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/policy.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/type.pyi'
+      ),
       subdir: 'gridfire',
   )
+
+  py_installation.install_sources(
+      files(
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/engine/__init__.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/engine/diagnostics.pyi',
+      ),
+      subdir: 'gridfire/engine',
+  )
+
+  py_installation.install_sources(
+      files(
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/__init__.pyi',
+      ),
+      subdir: 'gridfire/utils',
+  )
+
+  py_installation.install_sources(
+      files(
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/hashing/__init__.pyi',
+          meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/hashing/reaction.pyi',
+      ),
+      subdir: 'gridfire/utils/hashing',
+  )
+else
+  message('Python bindings disabled')
 endif
-
-
-py_installation.install_sources(
-    files(
-        meson.project_source_root() + '/src/python/gridfire/__init__.py',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/__init__.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/exceptions.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/partition.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/reaction.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/screening.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/io.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/solver.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/policy.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/type.pyi'
-    ),
-    subdir: 'gridfire',
-)
-
-py_installation.install_sources(
-    files(
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/engine/__init__.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/engine/diagnostics.pyi',
-    ),
-    subdir: 'gridfire/engine',
-)
-
-py_installation.install_sources(
-    files(
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/__init__.pyi',
-    ),
-    subdir: 'gridfire/utils',
-)
-
-py_installation.install_sources(
-    files(
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/hashing/__init__.pyi',
-        meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/hashing/reaction.pyi',
-    ),
-    subdir: 'gridfire/utils/hashing',
-)

cross/wasm.ini

@@ -0,0 +1,23 @@
+[binaries]
+c = 'emcc'
+cpp = 'em++'
+ar = 'emar'
+strip = 'emstrip'
+
+exec_wrapper = 'node'
+
+[built-in options]
+c_args = ['-Dpkg_config=false', '-Dbuild_tests=false', '-Dbuild_examples=true', '-Dbuild_fortran=falase', '-Dplugin_support=false', '-s', 'MEMORY64=1', '-pthread', '-DQUILL_NO_THREAD_NAME_SUPPORT', '-DQUILL_IMMEDIATE_FLUSH']
+cpp_args = ['-Dpkg_config=false', '-Dbuild_tests=false', '-Dbuild_examples=true', '-Dbuild_fortran=falase', '-Dplugin_support=false', '-s', 'MEMORY64=1', '-pthread', '-DQUILL_NO_THREAD_NAME_SUPPORT', '-DQUILL_IMMEDIATE_FLUSH']
+c_link_args = ['-s', 'WASM=1', '-s', 'ALLOW_MEMORY_GROWTH=1', '-s', 'MEMORY64=1', '-fwasm-exceptions', '-pthread', '-s', 'EXPORTED_RUNTIME_METHODS=["FS", "callMain"]', '-s', 'STACK_SIZE=10485760']
+cpp_link_args = ['-s', 'WASM=1', '-s', 'ALLOW_MEMORY_GROWTH=1', '-s', 'MEMORY64=1', '-fwasm-exceptions', '-pthread', '-s', 'EXPORTED_RUNTIME_METHODS=["FS", "callMain"]', '-s', 'STACK_SIZE=10485760']
+
+[host_machine]
+system = 'emscripten'
+cpu_family = 'wasm64'
+cpu = 'wasm64'
+endian = 'little'
+
+[properties]
+cmake_toolchain_file = '/home/tboudreaux/Programming/emsdk/upstream/emscripten/cmake/Modules/Platform/Emscripten.cmake'
+

meson.build

@@ -20,141 +20,27 @@
 # *********************************************************************** #
 project('GridFire', ['c', 'cpp'], version: 'v0.7.4_rc2', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
 
-if get_option('build-python')
-    add_project_arguments('-fvisibility=hidden', language: 'cpp')
-else
-    add_project_arguments('-fvisibility=default', language: 'cpp')
-endif
+# Start by running the code which validates the build environment
+subdir('build-check')
 
-message('Found CXX compiler: ' + meson.get_compiler('cpp').get_id())
-message('C++ standard set to: ' + get_option('cpp_std'))
-
-cppc = meson.get_compiler('cpp')
-cc = meson.get_compiler('c')
-
-
-if cppc.get_id() == 'clang'
-    message('disabling bitwise-instead-of-logical warnings for clang')
-    add_project_arguments('-Wno-bitwise-instead-of-logical', language: 'cpp')
-endif
-
-if cppc.get_id() == 'gcc'
-    message('disabling psabi warnings for gcc')
-    add_project_arguments('-Wno-psabi', language: 'cpp')
-    
-    if (cppc.version().version_compare('<14.0'))
-      error('g++ version must be at least 14.0, found ' + cppc.version())
-    endif
-endif
-
-build_fortran = get_option('build-fortran')
-if (build_fortran)
-  add_languages('fortran', native: true)
-  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)')
-  fc = meson.get_compiler('fortran')
-  if not get_option('unsafe-fortran')
-    if fc.get_id() != 'gcc'
-      error('The only supported fortran compiler for GridFire is gfortran (version >= 14.0), found ' + fc  + '. GridFire has not been tested with any other compilers. You can disable this check with the -Dunsafe-fortran=true flag to try other compilers')
-    endif
-  endif
-  if (fc.version().version_compare('<14.0'))
-      error('gfortran version must be at least 14.0, found ' + fc.version())
-  endif
-endif
-
-
-if not cppc.has_header('print')
-    error('C++ standard library header <print> not found. Please ensure your compiler and standard library supports C++23. We have already validated your compiler version so this is likely an issue with your standard library installation.')
-endif
-if not cppc.has_header('format')
-    error('C++ standard library header <format> not found. Please ensure your compiler and standard library supports C++23. We have already validated your compiler version so this is likely an issue with your standard library installation.')
-endif
-
-
-ignore_unused_args = '-Wno-unused-command-line-argument'
-
-add_global_arguments(ignore_unused_args, language: 'cpp')
-add_global_arguments(ignore_unused_args, language: 'c')
-
-
-
-# For Eigen
-add_project_arguments('-Wno-deprecated-declarations', language: 'cpp')
-
-llevel = get_option('log-level')
-
-logbase='QUILL_COMPILE_ACTIVE_LOG_LEVEL_'
-
-if (llevel == 'traceL3')
-    message('Setting log level to TRACE_L3')
-    log_argument = logbase + 'TRACE_L3'
-elif (llevel == 'traceL2')
-    message('Setting log level to TRACE_L2')
-    log_argument = logbase + 'TRACE_L2'
-elif (llevel == 'traceL1')
-    message('Setting log level to TRACE_L1')
-    log_argument = logbase + 'TRACE_L1'
-elif (llevel == 'debug')
-    message('Setting log level to DEBUG')
-    log_argument = logbase + 'DEBUG'
-elif (llevel == 'info')
-    message('Setting log level to INFO')
-    log_argument = logbase + 'INFO'
-elif (llevel == 'warning')
-    message('Setting log level to WARNING')
-    log_argument = logbase + 'WARNING'
-elif (llevel == 'error')
-    message('Setting log level to ERROR')
-    log_argument = logbase + 'ERROR'
-elif (llevel == 'critical')
-    message('Setting log level to CRITICAL')
-    log_argument = logbase + 'CRITICAL'
-endif
-
-log_argument = '-DQUILL_COMPILE_ACTIVE_LOG_LEVEL=' + log_argument
-add_project_arguments(log_argument, language: 'cpp')
-
-cpp = meson.get_compiler('cpp')
+# Configure the logging level
+subdir('build-extra/log-level')
 
+# Then build the external dependencies
 subdir('build-config')
 
+# Build the main source code
 subdir('src')
 
-if get_option('build-python')
-    message('Configuring Python bindings...')
-    subdir('build-python')
-else
-    message('Skipping Python bindings...')
-endif
+# Build the Python bindings
+subdir('build-python')
 
-if get_option('build-tests')
-    message('Setting up tests for GridFire...')
-    subdir('tests')
-else
-    message('Skipping tests for GridFire...')
-endif
-
-
-if get_option('pkg-config')
-    message('Generating pkg-config file for GridFire...')
-    pkg = import('pkgconfig')
-    pkg.generate(
-        name: 'gridfire',
-        description: 'GridFire nuclear reaction network solver',
-        version: meson.project_version(),
-        libraries: [
-            libgridfire,
-            libcomposition,
-            libconfig,
-            libconst,
-            liblogging
-        ],
-        subdirs: ['gridfire'],
-        filebase: 'gridfire',
-        install_dir: join_paths(get_option('libdir'), 'pkgconfig')
-    )
-endif
+# Buil the test suite
+subdir('tests')
+
+subdir('tools')
+
+# Build the pkg-config file
+subdir('build-extra/pkg-config')
 
 

meson_options.txt

@@ -1,8 +1,12 @@
-option('log-level', type: 'combo', choices: ['traceL3', 'traceL2', 'traceL1', 'debug', 'info', 'warning', 'error', 'critial'], value: 'info', description: 'Set the log level for the GridFire library')
-option('pkg-config', type: 'boolean', value: true, description: 'generate pkg-config file for GridFire (gridfire.pc)')
-option('build-python', type: 'boolean', value: false, description: 'build the python bindings so you can use GridFire from python')
-option('build-tests', type: 'boolean', value: true, description: 'build the test suite')
-option('build-fortran', type: 'boolean', value: false, description: 'build fortran module support')
-option('unsafe-fortran', type: 'boolean', value: false, description: 'Allow untested fortran compilers (compilers other than gfortran)')
-option('unity-safe', type: 'boolean', value: false, description: 'Enable safe unity builds for better compatibility across different compilers and platforms')
-option('python-target-version', type: 'string', value: '3.13', description: 'Target version for python compilation, only used for cross compilation')
+option('log_level', type: 'combo', choices: ['traceL3', 'traceL2', 'traceL1', 'debug', 'info', 'warning', 'error', 'critial'], value: 'info', description: 'Set the log level for the GridFire library')
+option('pkg_config', type: 'boolean', value: true, description: 'generate pkg-config file for GridFire (gridfire.pc)')
+option('build_python', type: 'boolean', value: false, description: 'build the python bindings so you can use GridFire from python')
+option('build_tests', type: 'boolean', value: true, description: 'build the test suite')
+option('build_examples', type: 'boolean', value: true, description: 'build example code')
+option('build_fortran', type: 'boolean', value: false, description: 'build fortran module support')
+option('unsafe_fortran', type: 'boolean', value: false, description: 'Allow untested fortran compilers (compilers other than gfortran)')
+option('plugin_support', type: 'boolean', value: false, description: 'Enable support for libplugin plugins')
+option('python_target_version', type: 'string', value: '3.13', description: 'Target version for python compilation, only used for cross compilation')
+option('build_c_api', type: 'boolean', value: true, description: 'compile the C API')
+option('build_tools', type: 'boolean', value: true, description: 'build the GridFire command line tools')
+option('openmp_support', type: 'boolean', value: false, description: 'Enable OpenMP support for parallelization')

src/extern/meson.build

@@ -23,7 +23,7 @@ gridfire_extern_dep = declare_dependency(
 
 install_subdir('include/gridfire', install_dir: get_option('includedir'))
 
-if get_option('build-fortran')
+if get_option('build_fortran')
     message('Configuring Fortran bindings...')
     subdir('fortran')
 endif

src/include/gridfire/config/config.h

@@ -0,0 +1,35 @@
+#pragma once
+
+#include "fourdst/config/config.h"
+
+namespace gridfire::config {
+    struct CVODESolverConfig {
+        double absTol = 1.0e-8;
+        double relTol = 1.0e-5;
+    };
+
+    struct SolverConfig {
+        CVODESolverConfig cvode;
+    };
+
+    struct AdaptiveEngineViewConfig {
+        double relativeCullingThreshold = 1.0e-75;
+    };
+
+    struct EngineViewConfig {
+        AdaptiveEngineViewConfig adaptiveEngineView;
+    };
+
+    struct EngineConfig {
+        EngineViewConfig views;
+
+    };
+
+    struct GridFireConfig {
+        SolverConfig solver;
+        EngineConfig engine;
+    };
+
+
+
+}

src/include/gridfire/engine/engine_abstract.h

@@ -53,7 +53,7 @@ namespace gridfire::engine {
     struct StepDerivatives {
         std::map<fourdst::atomic::Species, T> dydt{}; ///< Derivatives of abundances (dY/dt for each species).
         T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate (e.g., erg/g/s).
-        std::map<fourdst::atomic::Species, std::unordered_map<std::string, T>> reactionContributions{};
+        std::optional<std::map<fourdst::atomic::Species, std::unordered_map<std::string, T>>> reactionContributions = std::nullopt;
         T neutrinoEnergyLossRate = T(0.0); // (erg/g/s)
         T totalNeutrinoFlux = T(0.0);      // (neutrinos/g/s)
 

src/include/gridfire/engine/engine_graph.h

@@ -12,6 +12,7 @@
 #include "gridfire/screening/screening_types.h"
 #include "gridfire/partition/partition_abstract.h"
 #include "gridfire/engine/procedures/construction.h"
+#include "gridfire/config/config.h"
 
 #include <string>
 #include <unordered_map>
@@ -96,7 +97,7 @@ namespace gridfire::engine {
      *
      * @see engine_abstract.h
      */
-    class GraphEngine final : public DynamicEngine{
+    class GraphEngine final : public DynamicEngine {
     public:
         /**
          * @brief Constructs a GraphEngine from a composition.
@@ -753,6 +754,14 @@ namespace gridfire::engine {
         [[nodiscard]]
         SpeciesStatus getSpeciesStatus(const fourdst::atomic::Species &species) const override;
 
+        [[nodiscard]] bool get_store_intermediate_reaction_contributions() const {
+            return m_store_intermediate_reaction_contributions;
+        }
+
+        void set_store_intermediate_reaction_contributions(const bool value) {
+            m_store_intermediate_reaction_contributions = value;
+        }
+
 
     private:
         struct PrecomputedReaction {
@@ -846,8 +855,14 @@ namespace gridfire::engine {
             const reaction::Reaction& m_reaction;
             const GraphEngine& m_engine;
         };
+
+        struct PrecomputationKernelResults {
+            std::vector<double> dydt_vector;
+            double total_neutrino_energy_loss_rate{0.0};
+            double total_neutrino_flux{0.0};
+        };
     private:
-        Config& m_config = Config::getInstance();
+        Config<config::GridFireConfig> m_config;
         quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
 
         constants m_constants;
@@ -866,6 +881,7 @@ namespace gridfire::engine {
         mutable CppAD::ADFun<double> m_rhsADFun; ///< CppAD function for the right-hand side of the ODE.
         mutable CppAD::ADFun<double> m_epsADFun; ///< CppAD function for the energy generation rate.
         mutable CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations.
+        mutable std::vector<double> m_local_abundance_cache;
 
         bool m_has_been_primed = false; ///< Flag indicating if the engine has been primed.
 
@@ -879,6 +895,7 @@ namespace gridfire::engine {
 
         bool m_usePrecomputation = true; ///< Flag to enable or disable using precomputed reactions for efficiency. Mathematically, this should not change the results. Generally end users should not need to change this.
         bool m_useReverseReactions = true; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered.
+        bool m_store_intermediate_reaction_contributions = false; ///< Flag to enable or disable storing intermediate reaction contributions for debugging.
 
         BuildDepthType m_depth;
 
@@ -948,6 +965,42 @@ namespace gridfire::engine {
          */
         [[nodiscard]] bool validateConservation() const;
 
+        double compute_reaction_flow(
+            const std::vector<double> &local_abundances,
+            const std::vector<double> &screening_factors,
+            const std::vector<double> &bare_rates,
+            const std::vector<double> &bare_reverse_rates,
+            double rho,
+            size_t reactionCounter,
+            const reaction::Reaction &reaction,
+            size_t reactionIndex,
+            const PrecomputedReaction &precomputedReaction
+        ) const;
+
+        std::pair<double, double> compute_neutrino_fluxes(
+            double netFlow,
+            const reaction::Reaction &reaction) const;
+
+        PrecomputationKernelResults accumulate_flows_serial(
+            const std::vector<double>& local_abundances,
+            const std::vector<double>& screening_factors,
+            const std::vector<double>& bare_rates,
+            const std::vector<double>& bare_reverse_rates,
+            double rho,
+            const reaction::ReactionSet& activeReactions
+        ) const;
+
+#ifdef GRIDFIRE_USE_OPENMP
+        PrecomputationKernelResults accumulate_flows_parallel(
+            const std::vector<double>& local_abundances,
+            const std::vector<double>& screening_factors,
+            const std::vector<double>& bare_rates,
+            const std::vector<double>& bare_reverse_rates,
+            double rho,
+            const reaction::ReactionSet& activeReactions
+        ) const;
+#endif
+
 
         [[nodiscard]] StepDerivatives<double> calculateAllDerivativesUsingPrecomputation(
             const fourdst::composition::CompositionAbstract &comp,
@@ -1207,7 +1260,10 @@ namespace gridfire::engine {
                 const T nu_ij = static_cast<T>(reaction.stoichiometry(species));
                 const T dydt_increment = threshold_flag * molarReactionFlow * nu_ij;
                 dydt_vec[speciesIdx] += dydt_increment;
-                result.reactionContributions[species][std::string(reaction.id())] = dydt_increment;
+
+                if (m_store_intermediate_reaction_contributions) {
+                    result.reactionContributions.value()[species][std::string(reaction.id())] = dydt_increment;
+                }
             }
 
         }

src/include/gridfire/engine/views/engine_adaptive.h

@@ -4,6 +4,7 @@
 #include "gridfire/screening/screening_abstract.h"
 #include "gridfire/screening/screening_types.h"
 #include "gridfire/types/types.h"
+#include "gridfire/config/config.h"
 
 #include "fourdst/atomic/atomicSpecies.h"
 #include "fourdst/config/config.h"
@@ -386,10 +387,10 @@ namespace gridfire::engine {
          */
         [[nodiscard]] SpeciesStatus getSpeciesStatus(const fourdst::atomic::Species &species) const override;
     private:
-        using Config = fourdst::config::Config;
         using LogManager = fourdst::logging::LogManager;
-        /** @brief A reference to the singleton Config instance, used for retrieving configuration parameters. */
-        Config& m_config = Config::getInstance();
+
+        fourdst::config::Config<config::GridFireConfig> m_config;
+
         /** @brief A pointer to the logger instance, used for logging messages. */
         quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
 

src/include/gridfire/engine/views/engine_defined.h

@@ -6,6 +6,8 @@
 #include "gridfire/io/network_file.h"
 #include "gridfire/types/types.h"
 
+#include "gridfire/config/config.h"
+
 #include "fourdst/config/config.h"
 #include "fourdst/logging/logging.h"
 
@@ -365,9 +367,9 @@ namespace gridfire::engine {
         [[nodiscard]] std::string getNetworkFile() const { return m_fileName; }
         [[nodiscard]] const io::NetworkFileParser& getParser() const { return m_parser; }
     private:
-        using Config = fourdst::config::Config;
-        using LogManager = fourdst::logging::LogManager;
-        Config& m_config = Config::getInstance();
+        using LogManager = LogManager;
+        Config<config::GridFireConfig> m_config;
+
         quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
         std::string m_fileName;
         ///< Parser for the network file.

src/include/gridfire/io/network_file.h

@@ -2,6 +2,7 @@
 
 #include "fourdst/config/config.h"
 #include "fourdst/logging/logging.h"
+#include "gridfire/config/config.h"
 
 #include "quill/Logger.h"
 
@@ -101,9 +102,8 @@ namespace gridfire::io {
          */
         [[nodiscard]] ParsedNetworkData parse(const std::string& filename) const override;
     private:
-        using Config = fourdst::config::Config;
         using LogManager = fourdst::logging::LogManager;
-        Config& m_config = Config::getInstance();
+        fourdst::config::Config<config::GridFireConfig> m_config;
         quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
     };
 
@@ -141,9 +141,8 @@ namespace gridfire::io {
          */
         [[nodiscard]] ParsedNetworkData parse(const std::string& filename) const override;
     private:
-        using Config = fourdst::config::Config;
         using LogManager = fourdst::logging::LogManager;
-        Config& m_config = Config::getInstance();
+        fourdst::config::Config<config::GridFireConfig> m_config;
         quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
 
         std::string m_filename;

src/include/gridfire/reaction/reaction.h

@@ -809,6 +809,8 @@ namespace gridfire::reaction {
         std::vector<RateCoefficientSet> m_rates; ///< List of rate coefficient sets from each source.
         bool m_weak = false;
 
+        mutable std::unordered_map<double, double> m_cached_rates;
+
     private:
         /**
          * @brief Template implementation for calculating the total reaction rate.
@@ -876,6 +878,8 @@ namespace gridfire::reaction {
 
         [[nodiscard]] std::optional<std::unique_ptr<Reaction>> get(const std::string_view& id) const;
 
+        [[nodiscard]] std::unique_ptr<Reaction> get(size_t index) const;
+
         /**
          * @brief Removes a reaction from the set.
          * @param reaction The Reaction to remove.

src/include/gridfire/solver/strategies/CVODE_solver_strategy.h

@@ -4,6 +4,7 @@
 #include "gridfire/engine/engine_abstract.h"
 #include "gridfire/types/types.h"
 #include "gridfire/exceptions/exceptions.h"
+#include "gridfire/config/config.h"
 
 #include "fourdst/atomic/atomicSpecies.h"
 #include "fourdst/config/config.h"
@@ -237,13 +238,13 @@ namespace gridfire::solver {
         };
 
         struct CVODERHSOutputData {
-            std::map<fourdst::atomic::Species, std::unordered_map<std::string, double>> reaction_contribution_map;
+            std::optional<std::map<fourdst::atomic::Species, std::unordered_map<std::string, double>>> reaction_contribution_map;
             double neutrino_energy_loss_rate;
             double total_neutrino_flux;
         };
 
     private:
-        fourdst::config::Config& m_config = fourdst::config::Config::getInstance();
+        fourdst::config::Config<config::GridFireConfig> m_config;
         quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
         /**
          * @brief CVODE RHS C-wrapper that delegates to calculate_rhs and captures exceptions.

src/lib/engine/engine_graph.cpp

@@ -28,6 +28,10 @@
 #include "cppad/utility/sparse_rc.hpp"
 #include "cppad/utility/sparse_rcv.hpp"
 
+#ifdef GRIDFIRE_USE_OPENMP
+    #include <omp.h>
+#endif
+
 
 namespace {
     enum class REACLIB_WEAK_TYPES {
@@ -403,6 +407,167 @@ namespace gridfire::engine {
         return true; // All reactions passed the conservation check
     }
 
+    double GraphEngine::compute_reaction_flow(
+        const std::vector<double> &local_abundances,
+        const std::vector<double> &screening_factors,
+        const std::vector<double> &bare_rates,
+        const std::vector<double> &bare_reverse_rates,
+        const double rho,
+        const size_t reactionCounter,
+        const reaction::Reaction &reaction,
+        const size_t reactionIndex,
+        const PrecomputedReaction &precomputedReaction
+    ) const {
+        double forwardAbundanceProduct = 1.0;
+        for (size_t i = 0; i < precomputedReaction.unique_reactant_indices.size(); ++i) {
+            const size_t reactantIndex = precomputedReaction.unique_reactant_indices[i];
+            const int power = precomputedReaction.reactant_powers[i];
+
+            const double abundance = local_abundances[reactantIndex];
+
+            double factor;
+            if (power == 1) { factor = abundance; }
+            else if (power == 2)  { factor = abundance * abundance; }
+            else { factor = std::pow(abundance, power); }
+
+            if (!std::isfinite(factor)) {
+                LOG_CRITICAL(m_logger, "Non-finite factor encountered in forward abundance product for reaction '{}'. Check input abundances for validity.", reaction.id());
+                throw exceptions::BadRHSEngineError("Non-finite factor encountered in forward abundance product.");
+            }
+
+            forwardAbundanceProduct *= factor;
+        }
+
+        const double bare_rate = bare_rates.at(reactionCounter);
+
+        const double screeningFactor = screening_factors[reactionCounter];
+        const size_t numReactants = m_reactions[reactionIndex].reactants().size();
+        const size_t numProducts = m_reactions[reactionIndex].products().size();
+
+        const double forwardMolarReactionFlow = screeningFactor *
+                                   bare_rate *
+                                   precomputedReaction.symmetry_factor *
+                                   forwardAbundanceProduct *
+                                   std::pow(rho, numReactants >  1 ? static_cast<double>(numReactants) - 1 : 0.0);
+        if (!std::isfinite(forwardMolarReactionFlow)) {
+            LOG_CRITICAL(m_logger, "Non-finite forward molar reaction flow computed for reaction '{}'. Check input abundances and rates for validity.", reaction.id());
+            throw exceptions::BadRHSEngineError("Non-finite forward molar reaction flow computed.");
+        }
+
+
+        double reverseMolarReactionFlow = 0.0;
+        if (precomputedReaction.reverse_symmetry_factor != 0.0 and m_useReverseReactions) {
+            const double bare_reverse_rate = bare_reverse_rates.at(reactionCounter);
+
+            double reverseAbundanceProduct = 1.0;
+            for (size_t i = 0; i < precomputedReaction.unique_product_indices.size(); ++i) {
+                const size_t productIndex = precomputedReaction.unique_product_indices[i];
+                reverseAbundanceProduct *= std::pow(local_abundances[productIndex], precomputedReaction.product_powers[i]);
+            }
+
+            reverseMolarReactionFlow = screeningFactor *
+                                       bare_reverse_rate *
+                                       precomputedReaction.reverse_symmetry_factor *
+                                       reverseAbundanceProduct *
+                                       std::pow(rho, numProducts > 1 ? static_cast<double>(numProducts) - 1 : 0.0);
+        }
+
+        return forwardMolarReactionFlow - reverseMolarReactionFlow;
+    }
+
+    std::pair<double, double> GraphEngine::compute_neutrino_fluxes(
+        const double netFlow,
+        const reaction::Reaction &reaction
+    ) const {
+        if (reaction.type() == reaction::ReactionType::REACLIB_WEAK) {
+            const double q_abs = std::abs(reaction.qValue());
+            const REACLIB_WEAK_TYPES weakType = get_weak_reaclib_reaction_type(reaction);
+            double neutrino_loss_fraction = 0.0;
+            switch (weakType) {
+                case REACLIB_WEAK_TYPES::BETA_PLUS_DECAY:
+                    [[fallthrough]];
+                case REACLIB_WEAK_TYPES::BETA_MINUS_DECAY:
+                    neutrino_loss_fraction = 0.5; // Approximate 50% energy loss to neutrinos for beta decays
+                    break;
+                case REACLIB_WEAK_TYPES::ELECTRON_CAPTURE:
+                    [[fallthrough]];
+                case REACLIB_WEAK_TYPES::POSITRON_CAPTURE:
+                    neutrino_loss_fraction = 1.0;
+                    break;
+                default: ;
+            }
+
+            const double local_neutrino_loss =  netFlow * q_abs * neutrino_loss_fraction * m_constants.Na * m_constants.MeV_to_erg;
+            const double local_neutrino_flux = netFlow * m_constants.Na;
+
+            return {local_neutrino_loss, local_neutrino_flux};
+        }
+        return {0.0, 0.0};
+    }
+
+    GraphEngine::PrecomputationKernelResults GraphEngine::accumulate_flows_serial(
+        const std::vector<double> &local_abundances,
+        const std::vector<double> &screening_factors,
+        const std::vector<double> &bare_rates,
+        const std::vector<double> &bare_reverse_rates,
+        const double rho,
+        const reaction::ReactionSet &activeReactions
+    ) const {
+        PrecomputationKernelResults results;
+        results.dydt_vector.resize(m_networkSpecies.size(), 0.0);
+
+        std::vector<double> molarReactionFlows;
+        molarReactionFlows.reserve(m_precomputedReactions.size());
+
+        size_t reactionCounter = 0;
+
+        for (const auto& reaction : activeReactions) {
+            uint64_t reactionHash = utils::hash_reaction(*reaction);
+            const size_t reactionIndex = m_precomputedReactionIndexMap.at(reactionHash);
+            const PrecomputedReaction& precomputedReaction = m_precomputedReactions[reactionIndex];
+
+            double netFlow = compute_reaction_flow(
+                local_abundances,
+                screening_factors,
+                bare_rates,
+                bare_reverse_rates,
+                rho,
+                reactionCounter,
+                *reaction,
+                reactionIndex,
+                precomputedReaction);
+
+            molarReactionFlows.push_back(netFlow);
+
+            auto [local_neutrino_loss, local_neutrino_flux] = compute_neutrino_fluxes(netFlow, *reaction);
+            results.total_neutrino_energy_loss_rate += local_neutrino_loss;
+            results.total_neutrino_flux += local_neutrino_flux;
+
+            reactionCounter++;
+        }
+
+        LOG_TRACE_L3(m_logger, "Computed {} molar reaction flows for active reactions. Assembling these into RHS", molarReactionFlows.size());
+
+        reactionCounter = 0;
+        for (const auto& reaction: activeReactions) {
+            const size_t j = m_precomputedReactionIndexMap.at(utils::hash_reaction(*reaction));
+            const auto& precomp = m_precomputedReactions[j];
+            const double R_j = molarReactionFlows[reactionCounter];
+
+            for (size_t i = 0; i < precomp.affected_species_indices.size(); ++i) {
+                const size_t speciesIndex = precomp.affected_species_indices[i];
+
+                const int stoichiometricCoefficient = precomp.stoichiometric_coefficients[i];
+
+                const double dydt_increment = static_cast<double>(stoichiometricCoefficient) * R_j;
+                results.dydt_vector[speciesIndex] += dydt_increment;
+            }
+            reactionCounter++;
+        }
+
+        return results;
+    }
+
     double GraphEngine::calculateReverseRate(
         const reaction::Reaction &reaction,
         const double T9,
@@ -655,6 +820,7 @@ namespace gridfire::engine {
 
     }
 
+
     StepDerivatives<double> GraphEngine::calculateAllDerivativesUsingPrecomputation(
         const fourdst::composition::CompositionAbstract &comp,
         const std::vector<double> &bare_rates,
@@ -672,132 +838,43 @@ namespace gridfire::engine {
             T9,
             rho
         );
+        m_local_abundance_cache.clear();
+        for (const auto& species: m_networkSpecies) {
+            m_local_abundance_cache.push_back(comp.contains(species) ? comp.getMolarAbundance(species) : 0.0);
+        }
 
         StepDerivatives<double> result;
+        std::vector<double> dydt_scratch(m_networkSpecies.size(), 0.0);
 
-        // --- Optimized loop ---
-        std::vector<double> molarReactionFlows;
-        molarReactionFlows.reserve(m_precomputedReactions.size());
+#ifndef GRIDFIRE_USE_OPENMP
+        const auto [dydt_vector, total_neutrino_energy_loss_rate, total_neutrino_flux] = accumulate_flows_serial(
+            m_local_abundance_cache,
+            screeningFactors,
+            bare_rates,
+            bare_reverse_rates,
+            rho,
+            activeReactions
+        );
+        dydt_scratch = dydt_vector;
+        result.neutrinoEnergyLossRate = total_neutrino_energy_loss_rate;
+        result.totalNeutrinoFlux = total_neutrino_flux;
+#else
+        const auto [dydt_vector, total_neutrino_energy_loss_rate, total_neutrino_flux] = accumulate_flows_parallel(
+            m_local_abundance_cache,
+            screeningFactors,
+            bare_rates,
+            bare_reverse_rates,
+            rho,
+            activeReactions
+        );
+        dydt_scratch = dydt_vector;
+        result.neutrinoEnergyLossRate = total_neutrino_energy_loss_rate;
+        result.totalNeutrinoFlux = total_neutrino_flux;
+#endif
 
-        size_t reactionCounter = 0;
-        for (const auto& reaction : activeReactions) {
-            // --- Efficient lookup of only the active reactions ---
-            uint64_t reactionHash = utils::hash_reaction(*reaction);
-            const size_t reactionIndex = m_precomputedReactionIndexMap.at(reactionHash);
-            PrecomputedReaction precomputedReaction = m_precomputedReactions[reactionIndex];
-
-            // --- Forward abundance product ---
-            double forwardAbundanceProduct = 1.0;
-            for (size_t i = 0; i < precomputedReaction.unique_reactant_indices.size(); ++i) {
-                const size_t reactantIndex = precomputedReaction.unique_reactant_indices[i];
-                const fourdst::atomic::Species& reactant = m_networkSpecies[reactantIndex];
-                const int power = precomputedReaction.reactant_powers[i];
-
-                if (!comp.contains(reactant)) {
-                    forwardAbundanceProduct = 0.0;
-                    break; // No need to continue if one of the reactants has zero abundance
-                }
-                double factor = std::pow(comp.getMolarAbundance(reactant), power);
-                if (!std::isfinite(factor)) {
-                    LOG_CRITICAL(m_logger, "Non-finite factor encountered in forward abundance product for reaction '{}'. Check input abundances for validity.", reaction->id());
-                    throw exceptions::BadRHSEngineError("Non-finite factor encountered in forward abundance product.");
-                }
-                forwardAbundanceProduct *= std::pow(comp.getMolarAbundance(reactant), power);
-            }
-
-            const double bare_rate = bare_rates.at(reactionCounter);
-
-            const double screeningFactor = screeningFactors[reactionCounter];
-            const size_t numReactants = m_reactions[reactionIndex].reactants().size();
-            const size_t numProducts = m_reactions[reactionIndex].products().size();
-
-            // --- Forward reaction flow ---
-            const double forwardMolarReactionFlow =
-                    screeningFactor *
-                    bare_rate *
-                    precomputedReaction.symmetry_factor *
-                    forwardAbundanceProduct *
-                    std::pow(rho, numReactants >  1 ? static_cast<double>(numReactants) - 1 : 0.0);
-            if (!std::isfinite(forwardMolarReactionFlow)) {
-                LOG_CRITICAL(m_logger, "Non-finite forward molar reaction flow computed for reaction '{}'. Check input abundances and rates for validity.", reaction->id());
-                throw exceptions::BadRHSEngineError("Non-finite forward molar reaction flow computed.");
-            }
-
-
-            // --- Reverse reaction flow ---
-            // Only do this is the reaction has a non-zero reverse symmetry factor (i.e. is reversible)
-            double reverseMolarReactionFlow = 0.0;
-            if (precomputedReaction.reverse_symmetry_factor != 0.0 and m_useReverseReactions) {
-                const double bare_reverse_rate = bare_reverse_rates.at(reactionCounter);
-
-                double reverseAbundanceProduct = 1.0;
-                for (size_t i = 0; i < precomputedReaction.unique_product_indices.size(); ++i) {
-                    const size_t productIndex = precomputedReaction.unique_product_indices[i];
-                    const fourdst::atomic::Species& product = m_networkSpecies[productIndex];
-                    reverseAbundanceProduct *= std::pow(comp.getMolarAbundance(product), precomputedReaction.product_powers[i]);
-                }
-
-                reverseMolarReactionFlow = screeningFactor *
-                    bare_reverse_rate *
-                    precomputedReaction.reverse_symmetry_factor *
-                    reverseAbundanceProduct *
-                    std::pow(rho, numProducts > 1 ? static_cast<double>(numProducts) - 1 : 0.0);
-            }
-
-            molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow);
-
-            if (reaction->type() == reaction::ReactionType::REACLIB_WEAK) {
-                double q_abs = std::abs(reaction->qValue());
-                REACLIB_WEAK_TYPES weakType = get_weak_reaclib_reaction_type(*reaction);
-                double neutrino_loss_fraction = 0.0;
-                switch (weakType) {
-                    case REACLIB_WEAK_TYPES::BETA_PLUS_DECAY:
-                        [[fallthrough]];
-                    case REACLIB_WEAK_TYPES::BETA_MINUS_DECAY:
-                        neutrino_loss_fraction = 0.5; // Approximate 50% energy loss to neutrinos for beta decays
-                        break;
-                    case REACLIB_WEAK_TYPES::ELECTRON_CAPTURE:
-                        [[fallthrough]];
-                    case REACLIB_WEAK_TYPES::POSITRON_CAPTURE:
-                        neutrino_loss_fraction = 1.0;
-                        break;
-                    default: ;
-                }
-
-                double local_neutrino_loss =  molarReactionFlows.back() * q_abs * neutrino_loss_fraction * m_constants.Na * m_constants.MeV_to_erg;
-                double local_neutrino_flux = molarReactionFlows.back() * m_constants.Na;
-
-                result.totalNeutrinoFlux += local_neutrino_flux;
-                result.neutrinoEnergyLossRate += local_neutrino_loss;
-            }
-            reactionCounter++;
-        }
-
-        LOG_TRACE_L3(m_logger, "Computed {} molar reaction flows for active reactions. Assembling these into RHS", molarReactionFlows.size());
-
-        // --- Assemble molar abundance derivatives ---
-        for (const auto& species: m_networkSpecies) {
-            result.dydt[species] = 0.0; // Initialize the change in abundance for each network species to 0
-        }
-
-        reactionCounter = 0;
-        for (const auto& reaction: activeReactions) {
-            size_t j = m_precomputedReactionIndexMap.at(utils::hash_reaction(*reaction));
-            const auto& precomp = m_precomputedReactions[j];
-            const double R_j = molarReactionFlows[reactionCounter];
-
-            for (size_t i = 0; i < precomp.affected_species_indices.size(); ++i) {
-                const size_t speciesIndex = precomp.affected_species_indices[i];
-                const fourdst::atomic::Species& species = m_networkSpecies[speciesIndex];
-
-                const int stoichiometricCoefficient = precomp.stoichiometric_coefficients[i];
-
-                // Update the derivative for this species
-                double dydt_increment = static_cast<double>(stoichiometricCoefficient) * R_j;
-                result.dydt.at(species) += dydt_increment;
-                result.reactionContributions[species][std::string(reaction->id())] = dydt_increment;
-            }
-            reactionCounter++;
+        // load scratch into result.dydt
+        for (size_t i = 0; i < m_networkSpecies.size(); ++i) {
+            result.dydt[m_networkSpecies[i]] = dydt_scratch[i];
         }
 
         // --- Calculate the nuclear energy generation rate ---
@@ -1502,4 +1579,69 @@ namespace gridfire::engine {
 
         return true;
     }
+
+#ifdef GRIDFIRE_USE_OPENMP
+    GraphEngine::PrecomputationKernelResults GraphEngine::accumulate_flows_parallel(
+        const std::vector<double> &local_abundances,
+        const std::vector<double> &screening_factors,
+        const std::vector<double> &bare_rates,
+        const std::vector<double> &bare_reverse_rates,
+        const double rho,
+        const reaction::ReactionSet &activeReactions
+    ) const {
+        int n_threads = omp_get_max_threads();
+        std::vector<std::vector<double>> thread_local_dydt(n_threads, std::vector<double>(m_networkSpecies.size(), 0.0));
+
+        double total_neutrino_energy_loss_rate = 0.0;
+        double total_neutrino_flux = 0.0;
+
+        #pragma omp parallel for schedule(static) reduction(+:total_neutrino_energy_loss_rate, total_neutrino_flux)
+        for (size_t k = 0; k < activeReactions.size(); ++k) {
+            int t_id = omp_get_thread_num();
+            const auto& reaction = activeReactions[k];
+            const size_t reactionIndex = m_precomputedReactionIndexMap.at(utils::hash_reaction(reaction));
+            const PrecomputedReaction& precomputedReaction = m_precomputedReactions[reactionIndex];
+
+            double netFlow = compute_reaction_flow(
+                local_abundances,
+                screening_factors,
+                bare_rates,
+                bare_reverse_rates,
+                rho,
+                reactionIndex,
+                reaction,
+                reactionIndex,
+                precomputedReaction
+            );
+
+            auto [neutrinoEnergyLossRate, neutrinoFlux] = compute_neutrino_fluxes(
+                netFlow,
+                reaction
+            );
+
+            total_neutrino_energy_loss_rate += neutrinoEnergyLossRate;
+            total_neutrino_flux += neutrinoFlux;
+
+            for (size_t i = 0; i < precomputedReaction.affected_species_indices.size(); ++i) {
+                thread_local_dydt[t_id][precomputedReaction.affected_species_indices[i]] +=
+                    netFlow * precomputedReaction.stoichiometric_coefficients[i];
+            }
+
+        }
+        PrecomputationKernelResults results;
+        results.total_neutrino_energy_loss_rate = total_neutrino_energy_loss_rate;
+        results.total_neutrino_flux = total_neutrino_flux;
+
+        results.dydt_vector.resize(m_networkSpecies.size(), 0.0);
+        #pragma omp parallel for schedule(static)
+        for (size_t i = 0; i < m_networkSpecies.size(); ++i) {
+            double sum = 0.0;
+            for (int t = 0; t < n_threads; ++t) sum += thread_local_dydt[t][i];
+            results.dydt_vector[i] = sum;
+        }
+
+        return results;
+    }
+#endif
+
 }

src/lib/engine/views/engine_adaptive.cpp

@@ -394,7 +394,9 @@ namespace gridfire::engine {
         const double maxFlow
     ) const {
         LOG_TRACE_L1(m_logger, "Culling reactions based on flow rates...");
-        const auto relative_culling_threshold = m_config.get<double>("gridfire:AdaptiveEngineView:RelativeCullingThreshold", 1e-75);
+
+        const auto relative_culling_threshold = m_config->engine.views.adaptiveEngineView.relativeCullingThreshold;
+
         double absoluteCullingThreshold = relative_culling_threshold * maxFlow;
         LOG_DEBUG(m_logger, "Relative culling threshold: {:7.3E} ({:7.3E})", relative_culling_threshold, absoluteCullingThreshold);
         std::vector<const reaction::Reaction*> culledReactions;

src/lib/reaction/reaction.cpp

@@ -278,7 +278,12 @@ namespace gridfire::reaction {
         double Ye,
         double mue, const std::vector<double> &Y, const std::unordered_map<size_t, Species>& index_to_species_map
     ) const {
-        return calculate_rate<double>(T9);
+        if (m_cached_rates.contains(T9)) {
+            return m_cached_rates.at(T9);
+        }
+        const double rate = calculate_rate<double>(T9);
+        m_cached_rates[T9] = rate;
+        return rate;
     }
 
     double LogicalReaclibReaction::calculate_log_rate_partial_deriv_wrt_T9(
@@ -455,6 +460,10 @@ namespace gridfire::reaction {
         return std::make_optional(m_reactions[m_reactionNameMap.at(std::string(id))]->clone());
     }
 
+    std::unique_ptr<Reaction> ReactionSet::get(size_t index) const {
+        return m_reactions.at(index)->clone();
+    }
+
     void ReactionSet::remove_reaction(const Reaction& reaction) {
         const size_t rh = reaction.hash(0);
         if (!m_reactionHashes.contains(rh)) {

src/lib/solver/strategies/CVODE_solver_strategy.cpp

@@ -112,8 +112,8 @@ namespace gridfire::solver {
         //  2. If the user has set tolerances in code, those override the config
         //  3. If the user has not set tolerances in code and the config does not have them, use hardcoded defaults
 
-        auto absTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:absTol", 1.0e-8);
-        auto relTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:relTol", 1.0e-5);
+        auto absTol = m_config->solver.cvode.absTol;
+        auto relTol = m_config->solver.cvode.relTol;
 
         if (m_absTol) {
             absTol = *m_absTol;
@@ -935,8 +935,8 @@ namespace gridfire::solver {
         sunrealtype *y_data = N_VGetArrayPointer(m_Y);
         sunrealtype *y_err_data = N_VGetArrayPointer(m_YErr);
 
-        const auto absTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:absTol", 1.0e-8);
-        const auto relTol = m_config.get<double>("gridfire:solver:CVODESolverStrategy:relTol", 1.0e-8);
+        const auto absTol = m_config->solver.cvode.absTol;
+        const auto relTol = m_config->solver.cvode.relTol;
 
         std::vector<double> err_ratios;
         const size_t num_components = N_VGetLength(m_Y);

src/meson.build

@@ -42,10 +42,15 @@ gridfire_build_dependencies = [
     json_dep,
 ]
 
-if not get_option('unity-safe')
+if get_option('plugin_support')
     gridfire_build_dependencies += [plugin_dep]
 endif
 
+if get_option('openmp_support')
+    openmp_dep = dependency('openmp', required: true)
+    gridfire_build_dependencies += [openmp_dep]
+endif
+
 # Define the libnetwork library so it can be linked against by other parts of the build system
 libgridfire = library('gridfire',
     gridfire_sources,
@@ -63,12 +68,11 @@ gridfire_dep = declare_dependency(
 
 install_subdir('include/gridfire', install_dir: get_option('includedir'))
 
-message('Configuring C API...')
-subdir('extern')
-#
-#if get_option('build-python')
-#    message('Configuring Python bindings...')
-#    subdir('python')
-#else
-#    message('Skipping Python bindings...')
-#endif
+
+if not get_option('build_c_api') and get_option('build_fortran')
+  error('Cannot build fortran without C API. Set -Dbuild-c-api=true and -Dbuild-fortran=true')
+endif
+if get_option('build_c_api')
+  message('Configuring C API...')
+  subdir('extern')
+endif

subprojects/fourdst.wrap

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

tests/extern/meson.build

@@ -1,5 +1,7 @@
-subdir('C')
+if get_option('build_c_api')
+  subdir('C')
+endif
 
-if get_option('build-fortran')
+if get_option('build_fortran')
   subdir('fortran')
 endif

tests/graphnet_sandbox/meson.build

@@ -1,5 +1,5 @@
 executable(
     'graphnet_sandbox',
     'main.cpp',
-    dependencies: [gridfire_dep, composition_dep,  cli11_dep],
+    dependencies: [gridfire_dep, cli11_dep],
 )

tests/meson.build

@@ -1,7 +1,7 @@
 # Google Test dependency
-gtest_dep = dependency('gtest', main: true, required : true)
-gtest_main = dependency('gtest_main', required: true)
-gtest_nomain_dep = dependency('gtest', main: false, required : true)
+#gtest_dep = dependency('gtest', main: true, required : true)
+#gtest_main = dependency('gtest_main', required: true)
+#gtest_nomain_dep = dependency('gtest', main: false, required : true)
 
 # Subdirectories for unit and integration tests
 subdir('graphnet_sandbox')

tools/config/generate_config_files.cpp

@@ -0,0 +1,48 @@
+#include "fourdst/config/config.h"
+#include "gridfire/config/config.h"
+#include <source_location>
+#include <filesystem>
+
+#include  "CLI/CLI.hpp"
+
+consteval std::string_view strip_namespaces(const std::string_view fullName) {
+    const size_t pos = fullName.rfind("::");
+    if (pos == std::string_view::npos) {
+        return fullName;
+    }
+    return fullName.substr(pos + 2);
+}
+
+template <typename T>
+consteval std::string_view get_type_name() {
+    constexpr std::string_view name = std::source_location::current().function_name();
+    const auto pos = name.find("T = ");
+    if (pos == std::string_view::npos) return name;
+    const auto start = pos + 4;
+    const auto end = name.rfind(']');
+    return name.substr(start, end - start);
+}
+
+int main(int argc, char** argv) {
+    CLI::App app{"GridFire Sandbox Application."};
+
+    std::string outputPath = ".";
+
+    app.add_option("-p,--path", outputPath, "path to save generated config files (default: current directory)");
+
+    CLI11_PARSE(app, argc, argv);
+
+    const std::filesystem::path outPath(outputPath);
+    if (!std::filesystem::exists(outPath)) {
+        std::cerr << "Error: The specified path does not exist: " << outputPath << std::endl;
+        return 1;
+    }
+
+    fourdst::config::Config<gridfire::config::GridFireConfig> configConfig;
+    const std::string_view name = strip_namespaces(get_type_name<gridfire::config::GridFireConfig>());
+
+    const std::string defaultConfigFilePath = (outPath / (std::string(name) + ".toml")).string();
+    const std::string schemaFilePath = (outPath / (std::string(name) + ".schema.json")).string();
+    configConfig.save(defaultConfigFilePath);
+    configConfig.save_schema(schemaFilePath);
+}

tools/config/meson.build

@@ -0,0 +1 @@
+executable('gf_generate_config_file', 'generate_config_files.cpp', dependencies: [gridfire_dep, cli11_dep], install: true)

tools/meson.build

@@ -0,0 +1,3 @@
+if get_option('build_tools')
+    subdir('config')
+endif