burnman-cpp 0.1.0

C++ translation of burnman (https://github.com/geodynamics/burnman)

Requires: C++17, Eigen >= 3.4, GSL Tests/Benchmarks require: Catch2 >= 3.4

Install: make Build tests/benchmarks: make test

Run benchmarks with: ./bin/run_tests [!benchmark] Save output to xlm also with: ./bin/run_tests [!benchmark] --reporter XML::out=./benchmark-report.xml --reporter console::out=-::colour-mode=ansi Use parse_benchmarks.py to process xml output

Python wrapper

The burnman_cpp package exposes the C++ implementation through pybind11. It includes minerals, assemblages, ideal and symmetric/asymmetric regular solutions, and constrained equilibration. This package has its own import name and can be installed alongside the original Python burnman package.

Install with Python >= 3.10, a C++17 compiler, CMake >= 3.18, Eigen >= 3.4 and GSL. The existing Makefile remains available for C++ builds.

# macOS native dependencies
brew install eigen gsl cmake
# Debian/Ubuntu native dependencies
# sudo apt install libeigen3-dev libgsl-dev cmake g++

python3 -m venv .venv
source .venv/bin/activate
python -m pip install .
python examples/solution_equilibrium.py

For development and tests:

python -m pip install -e '.[test]'
python -m pytest

Editable installs require reinstalling after changing C++ bindings or sources. Nonstandard Eigen/GSL installations can be supplied through CMAKE_PREFIX_PATH, or -Ccmake.define.Eigen3_DIR=/path/to/eigen/cmake and -Ccmake.define.GSL_ROOT_DIR=/path/to/gsl. Python build dependencies are installed automatically by pip. The build follows pybind11's documented CMake and scikit-build-core packaging approach.

Solutions

Create endmember minerals from parameter dictionaries or MineralParams. The dictionaries accept the original BurnMan names n, Debye_0 and EOS strings such as "slb3"; the C++ names napfu, debye_0 and EOSType enums are also accepted. Cp accepts a sequence of four coefficients. Mineral parameters are validated and copied at construction. The params property returns a snapshot; construct a new mineral to change parameters.

Given periclase and wuestite minerals, as constructed in the complete example:

import burnman_cpp as bm

model = bm.SymmetricRegularSolution(
    [(periclase, "[Mg]O"), (wuestite, "[Fe]O")],
    energy_interaction=[[13.0e3]],
)
oxide = bm.Solution(model, [0.8, 0.2], name="Ferropericlase")
oxide.set_state(25.0e9, 2000.0)
print(oxide.density, oxide.molar_gibbs, oxide.activities)
oxide.set_composition([0.6, 0.4])

IdealSolution(endmembers) and AsymmetricRegularSolution(endmembers, alphas, energy_interaction, ...) are also available. Interaction arguments are upper triangular lists with row lengths n-1, n-2, ..., 1. Optional volume_interaction and entropy_interaction have the same format. Each Solution owns a copy of its model's endmembers, so sharing a model between solutions preserves independent pressure and temperature states.

Material properties include volume, density, Gibbs/Helmholtz/internal energy, entropy, enthalpy, heat capacities, elastic moduli, thermal expansivity and seismic velocities. Solutions also expose activities, activity coefficients, excess properties, partial molar properties and Gibbs/entropy/volume Hessians. Properties are available as Python attributes and C++-style get_*() methods. Call set_state() before reading state properties. set_composition() clears the cached solution properties. Eigen arrays and matrices return NumPy copies.

All quantities use SI units: Pa, K, kg/mol, m³/mol and J/mol.

Equilibration

assemblage = bm.Assemblage([oxide], [1.0])
result = bm.equilibrate(
    {"Mg": 0.6, "Fe": 0.4, "O": 1.0},
    assemblage,
    [bm.PressureConstraint(25.0e9), bm.TemperatureConstraint(2000.0)],
    store_iterates=True,
)
solve = result.sol_array[0, 0]
print(solve.success, solve.message, solve.F_norm)
print(result.prm.parameter_names, solve.x)

Supply exactly two equality constraints plus one per free compositional vector. Equilibration uses flat assemblages of mineral and solution phases. Bulk composition dictionaries must include every assemblage element, including elements with zero amount. PressureConstraint uses Pa and TemperatureConstraint uses K. EntropyConstraint uses total entropy in J/K and VolumeConstraint uses total volume in m³ for the specified bulk amount. PTEllipseConstraint, LinearXConstraint, PhaseFractionConstraint and PhaseCompositionConstraint are also exposed. Build the latter two with the parameters from get_equilibration_parameters(assemblage, composition, free_vectors).

A list of constraints in an argument generates a Cartesian grid:

result = bm.equilibrate(
    {"Mg": 0.6, "Fe": 0.4, "O": 1.0}, assemblage,
    [[bm.PressureConstraint(p) for p in (10.0e9, 25.0e9)],
     [bm.TemperatureConstraint(t) for t in (1500.0, 2000.0)]],
)
assert result.sol_array.shape == (2, 2)

sol_array is a NumPy object array of DampedNewtonResult objects containing success, code, message, x, F, J, F_norm, n_iterations and optional iteration_history. Check success for each solve. tol, max_iterations and verbose control the native solver.

Equilibration mutates the assemblage and its phase objects. The C++ solver does not yet store separate assemblage snapshots. Apply a selected grid point using set_composition_and_state_from_parameters(assemblage, solve.x). get_parameter_vector() and get_endmember_amounts() expose the corresponding vectors. Assemblages default to one mole and Voigt-Reuss-Hill averaging; use n_moles, set_fractions() and set_averaging_scheme() to adjust them. After changing a phase directly, call assemblage.reset_cache() before querying cached aggregate properties. Python EOS subclasses, property modifiers and the original Python mineral database are outside this initial API.

For C++ library-only builds through CMake:

cmake -S . -B build/cmake -DBURNMAN_BUILD_PYTHON=OFF
cmake --build build/cmake