diff --git a/.github/workflows/julia-interface.yml b/.github/workflows/julia-interface.yml new file mode 100644 index 0000000000..c7613892c4 --- /dev/null +++ b/.github/workflows/julia-interface.yml @@ -0,0 +1,89 @@ +# CI for the experimental Julia interface (interfaces/julia). +# +# Builds libcantera *from this checkout* (so the generated CLib matches the +# bindings, including functions this branch adds/renames), then runs the Julia +# test suite against the freshly built library on Linux. +name: Julia interface + +on: + push: + paths: + - 'interfaces/julia/**' + - 'interfaces/sourcegen/**' + - 'interfaces/clib/**' + - '.github/workflows/julia-interface.yml' + pull_request: + paths: + - 'interfaces/julia/**' + - 'interfaces/sourcegen/**' + - 'interfaces/clib/**' + - '.github/workflows/julia-interface.yml' + workflow_dispatch: + +permissions: + contents: read + +jobs: + test: + name: Julia ${{ matrix.julia-version }} on ubuntu-latest + runs-on: ubuntu-latest + defaults: + run: + shell: bash -el {0} + strategy: + fail-fast: false + matrix: + julia-version: ['1'] + steps: + - uses: actions/checkout@v4 + with: + persist-credentials: false + + - name: Set up build dependencies (conda) + uses: conda-incubator/setup-miniconda@fc2d68f6413eb2d87b895e92f8584b5b94a10167 # v3.3.0 + with: + miniforge-version: latest + activate-environment: build + python-version: '3.12' + channels: conda-forge + + - name: Install Cantera build dependencies + run: | + # Build toolchain + C++ deps, plus sourcegen's Python deps (jinja2, + # ruamel.yaml) which are needed to regenerate the CLib during the build. + conda install -y scons cython "eigen" sundials fmt yaml-cpp hdf5 \ + "boost-cpp" doxygen openblas jinja2 "ruamel.yaml" + + - name: Build libcantera (shared) + run: | + cat > cantera.conf <"] +version = "0.1.0" + +[compat] +julia = "1.9" + +[extras] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["LinearAlgebra", "Test"] diff --git a/interfaces/julia/README.md b/interfaces/julia/README.md new file mode 100644 index 0000000000..e1f3aaa4e7 --- /dev/null +++ b/interfaces/julia/README.md @@ -0,0 +1,47 @@ +# Cantera.jl + +A Julia interface to [Cantera](https://cantera.org) for chemical +kinetics, thermodynamics, and transport. It links Cantera's native `libcantera` +directly through the generated CLib API — it is **not** a reimplementation of +Cantera and has **no** Python dependency. + +## Install + +Requires Julia 1.9+ and a compiled `libcantera`. Point the package at the +library (a directory or a full path) and the mechanism data: + +```bash +export CANTERA_LIBRARY_PATH=/path/to/cantera/lib # or use a conda env +export CANTERA_DATA=/path/to/cantera/data # for gri30.yaml, etc. +``` + +The low-level CLib bindings are scaffolded by `scons build` and are kept up to +date by the build process. + +```julia +using Pkg +Pkg.activate("interfaces/julia") +Pkg.instantiate() +``` + +## Quickstart + +```julia +using Cantera +gas = Solution("gri30.yaml") +set_TPX!(gas, 1000.0, one_atm, "H2:2, O2:1, N2:4") +net = ReactorNet(IdealGasReactor(gas)) +advance!(net, 1e-3) +temperature(gas) # ignited temperature [K] +``` + +## Documentation + +Exported functions carry docstrings, accessible from the Julia REPL help mode +(type `?` followed by the function name). See the [`examples/`](examples) +directory for runnable usage. + +## Status + +Experimental. This interface targets Cantera's experimental CLib backend, and +APIs may change. diff --git a/interfaces/julia/SConscript b/interfaces/julia/SConscript new file mode 100644 index 0000000000..29448f8e33 --- /dev/null +++ b/interfaces/julia/SConscript @@ -0,0 +1,26 @@ +"""SConscript file for the generated Julia CLib bindings.""" +from pathlib import Path + +Import("env", "build") + +# Generated file names can be anticipated from the CLib specifications +auto_path = Path(Dir("#interfaces/sourcegen/src/sourcegen/headers").abspath) +yaml_files = sorted(auto_path.glob("ct*.yaml")) +generated_files = ["#interfaces/julia/src/generated/_manifest.jl"] +generated_files.extend( + f"#interfaces/julia/src/generated/lib{yaml_file.stem}.jl" + for yaml_file in yaml_files) + +sourcegen = env.Command( + generated_files, + "#build/doc/Cantera.tag", + ("$python_cmd -m interfaces.sourcegen.src.sourcegen " + f"--api=julia --output={Path(Dir('#interfaces/julia').abspath)}") +) +env.Depends(sourcegen, [File(ff) for ff in yaml_files]) +env.Depends(sourcegen, env.Glob("#interfaces/sourcegen/src/sourcegen/*.py")) +env.Depends(sourcegen, env.Glob("#interfaces/sourcegen/src/sourcegen/julia/*.*")) + +build(sourcegen) + +Return("sourcegen") diff --git a/interfaces/julia/docs/Project.toml b/interfaces/julia/docs/Project.toml new file mode 100644 index 0000000000..66ce5ee495 --- /dev/null +++ b/interfaces/julia/docs/Project.toml @@ -0,0 +1,6 @@ +[deps] +Cantera = "c8cccd67-5bbd-4b73-b352-8d151246ed36" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "1" diff --git a/interfaces/julia/docs/make.jl b/interfaces/julia/docs/make.jl new file mode 100644 index 0000000000..7467b4b8f6 --- /dev/null +++ b/interfaces/julia/docs/make.jl @@ -0,0 +1,32 @@ +# Build the Cantera.jl documentation with Documenter. +# +# julia --project=interfaces/julia/docs interfaces/julia/docs/make.jl +# +# The generated HTML is written to `docs/build`, which the Cantera doc build +# (doc/SConscript) copies into `build/doc/html/julia` so it is deployed +# alongside the rest of the Cantera documentation. There is intentionally no +# `deploydocs` call: deployment is handled by the top-level docs workflow. +# +# Requires a working libcantera and the generated CLib bindings (see the +# interface README); `using Cantera` loads the library to introspect docstrings. + +using Documenter +using Cantera + +makedocs( + sitename = "Cantera.jl", + modules = [Cantera], + authors = "Cantera Developers", + pages = [ + "Home" => "index.md", + "API Reference" => "reference.md", + ], + format = Documenter.HTML( + prettyurls = get(ENV, "CI", "false") == "true", + canonical = "https://cantera.org/dev/julia", + # The API reference is a single @autodocs page; allow it to exceed the + # default 200 KiB cap rather than splitting the generated listing. + size_threshold = 512 * 1024, + ), + warnonly = true, # the CLib backend is experimental; don't fail on missing refs +) diff --git a/interfaces/julia/docs/src/index.md b/interfaces/julia/docs/src/index.md new file mode 100644 index 0000000000..44c6de6ac0 --- /dev/null +++ b/interfaces/julia/docs/src/index.md @@ -0,0 +1,43 @@ +# Cantera.jl + +A Julia interface to [Cantera](https://cantera.org) for chemical kinetics, +thermodynamics, and transport. It links Cantera's native `libcantera` directly +through the generated CLib API — it is **not** a reimplementation of Cantera and +has **no** Python dependency. + +!!! warning "Experimental" + This interface targets Cantera's experimental CLib backend. The API may + change between releases. + +## Installation + +Requires Julia 1.9+ and a compiled `libcantera`. Point the package at the +library and the mechanism data: + +```bash +export CANTERA_LIBRARY_PATH=/path/to/cantera/lib # or use a conda env +export CANTERA_DATA=/path/to/cantera/data # for gri30.yaml, etc. +``` + +The low-level CLib bindings are scaffolded by `scons build`, so only the +environment needs to be instantiated: + +```bash +julia --project=interfaces/julia -e 'using Pkg; Pkg.instantiate()' +``` + +## Quickstart + +```julia +using Cantera + +gas = Solution("gri30.yaml") +set_TPX!(gas, 1000.0, one_atm, "H2:2, O2:1, N2:4") +net = ReactorNet(IdealGasReactor(gas)) +advance!(net, 1e-3) +temperature(gas) # ignited temperature [K] +``` + +See the [`examples/`](https://github.com/Cantera/cantera/tree/main/interfaces/julia/examples) +directory for more runnable usage, and the [API Reference](reference.md) for the +full list of documented functions. diff --git a/interfaces/julia/docs/src/reference.md b/interfaces/julia/docs/src/reference.md new file mode 100644 index 0000000000..77c319c3b8 --- /dev/null +++ b/interfaces/julia/docs/src/reference.md @@ -0,0 +1,7 @@ +# API Reference + +Documentation for the functions and types exported by `Cantera.jl`. + +```@autodocs +Modules = [Cantera] +``` diff --git a/interfaces/julia/examples/basic.jl b/interfaces/julia/examples/basic.jl new file mode 100644 index 0000000000..992e4bddb8 --- /dev/null +++ b/interfaces/julia/examples/basic.jl @@ -0,0 +1,21 @@ +# Basic usage of the Cantera Julia interface. +# +# Run from the repository root with: +# CANTERA_LIBRARY_PATH=/path/to/lib julia --project=interfaces/julia \ +# interfaces/julia/examples/basic.jl + +using Cantera + +gas = Solution("gri30.yaml") +set_TPY!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") + +println("T = ", temperature(gas), " K") +println("p = ", pressure(gas), " Pa") +println("rho = ", density(gas), " kg/m^3") +println("cp = ", cp_mass(gas), " J/kg/K") +println("lambda = ", thermal_conductivity(gas), " W/m/K") +println("mu = ", viscosity(gas), " Pa*s") + +wdot = net_production_rates(gas) +ich4 = species_index(gas, "CH4") +println("wdot(CH4) = ", wdot[ich4], " kmol/m^3/s") diff --git a/interfaces/julia/examples/flame.jl b/interfaces/julia/examples/flame.jl new file mode 100644 index 0000000000..95ef32c974 --- /dev/null +++ b/interfaces/julia/examples/flame.jl @@ -0,0 +1,19 @@ +# Freely-propagating premixed laminar flame, mirroring Cantera's `adiabatic_flame` +# example. Computes the laminar flame speed of a stoichiometric methane/air +# mixture on the GRI-Mech 3.0 mechanism. + +using Cantera + +gas = Solution("gri30.yaml") +set_TPX!(gas, 300.0, one_atm, "CH4:1, O2:2, N2:7.52") + +flame = FreeFlame(gas; width=0.03) +solve!(flame; auto=true) # staged energy-off -> energy-on -> refine + +println("laminar flame speed Su = ", flame_speed(flame), " m/s") +println("adiabatic flame T = ", maximum(flame_T(flame)), " K") +println("grid points = ", length(grid(flame))) + +# Peak OH mole fraction through the flame front +oh = flame_X(flame, "OH") +println("peak X(OH) = ", maximum(oh)) diff --git a/interfaces/julia/examples/reactor.jl b/interfaces/julia/examples/reactor.jl new file mode 100644 index 0000000000..0369923548 --- /dev/null +++ b/interfaces/julia/examples/reactor.jl @@ -0,0 +1,22 @@ +# Constant-volume adiabatic ignition of a hydrogen/air mixture, integrated with +# a ReactorNet. Mirrors Cantera's canonical `reactor1.py` example. + +using Cantera + +gas = Solution("gri30.yaml") +set_TPX!(gas, 1000.0, one_atm, "H2:2, O2:1, N2:4") + +reactor = IdealGasReactor(gas) +net = ReactorNet(reactor) + +println("# t [s] T [K] p [Pa]") +t = 0.0 +tend = 1e-3 +while t < tend + global t = step!(net) + println(rpad(round(t, sigdigits=4), 10), " ", + rpad(round(temperature(reactor), digits=2), 10), " ", + round(pressure(reactor), digits=1)) +end + +println("\nFinal temperature: ", temperature(reactor), " K") diff --git a/interfaces/julia/examples/reactor_network.jl b/interfaces/julia/examples/reactor_network.jl new file mode 100644 index 0000000000..4cf4f2a917 --- /dev/null +++ b/interfaces/julia/examples/reactor_network.jl @@ -0,0 +1,37 @@ +# Reactor network with mass flow: a well-stirred-reactor-like setup where a hot +# reservoir feeds an IdealGasReactor through a MassFlowController and the reactor +# exhausts to a downstream reservoir through a Valve. Demonstrates reservoirs, +# flow devices, and time integration. + +using Cantera + +# Upstream reservoir: fresh reactants. +upstream_gas = Solution("gri30.yaml") +set_TPX!(upstream_gas, 300.0, one_atm, "CH4:1, O2:2, N2:7.52") +upstream = Reservoir(upstream_gas) + +# The reactor, initially filled with hot products to ignite the incoming flow. +reactor_gas = Solution("gri30.yaml") +set_TPX!(reactor_gas, 1800.0, one_atm, "CH4:1, O2:2, N2:7.52") +equilibrate!(reactor_gas, "TP") +reactor = IdealGasReactor(reactor_gas) + +# Downstream exhaust reservoir. +downstream_gas = Solution("gri30.yaml") +set_TPX!(downstream_gas, 300.0, one_atm, "N2:1") +downstream = Reservoir(downstream_gas) + +# Connect: fixed inlet mass flow, pressure-controlled outlet. +mfc = MassFlowController(upstream, reactor; mdot=0.05) +valve = Valve(reactor, downstream; K=1.0) + +net = ReactorNet(reactor) +set_tolerances!(net; rtol=1e-8, atol=1e-14) + +println("# t [s] T [K] inlet mdot [kg/s]") +for t in range(0, 2e-2; length=6) + advance!(net, t) + println(rpad(round(t, sigdigits=3), 10), " ", + rpad(round(temperature(reactor), digits=1), 10), " ", + round(mass_flow_rate(mfc), sigdigits=4)) +end diff --git a/interfaces/julia/src/Cantera.jl b/interfaces/julia/src/Cantera.jl new file mode 100644 index 0000000000..d0a694d765 --- /dev/null +++ b/interfaces/julia/src/Cantera.jl @@ -0,0 +1,51 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +""" + Cantera + +Julia interface to [Cantera](https://cantera.org), built on top of +Cantera's generated CLib API. + +```julia +using Cantera +gas = Solution("gri30.yaml") +set_TPY!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") +temperature(gas), pressure(gas), net_production_rates(gas) +``` +""" +module Cantera + +include("LibCantera.jl") +using .LibCantera: LibCantera, libcantera + +include("errors.jl") +include("handles.jl") +include("solution.jl") +include("thermo.jl") +include("kinetics.jl") +include("transport.jl") +include("reaction.jl") +include("func1.jl") +include("reactor.jl") +include("reactornet.jl") +include("connectors.jl") +include("onedim.jl") +include("multiphase.jl") +include("rdiag.jl") +include("solutionarray.jl") +include("utils.jl") + +# Register Cantera's bundled data directory so mechanisms can be loaded by name. +function __init__() + d = LibCantera.default_data_directory() + if d !== nothing + try + LibCantera.ct_addDataDirectory(d) + catch + end + end + return nothing +end + +end # module Cantera diff --git a/interfaces/julia/src/LibCantera.jl b/interfaces/julia/src/LibCantera.jl new file mode 100644 index 0000000000..3e6d176a93 --- /dev/null +++ b/interfaces/julia/src/LibCantera.jl @@ -0,0 +1,109 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +""" + LibCantera + +Internal module holding the low-level `ccall` bindings to Cantera's generated +CLib API. Everything here is an implementation detail: user code should use the +API exported from the top-level [`Cantera`](@ref) module instead of the +raw handle-based functions defined here. + +The wrapper functions in `generated/` are scaffolded by `sourcegen` from the +CLib specifications as part of the Cantera build process and must not be +hand-edited. This module is responsible only for: + + * locating the Cantera shared library (`libcantera`); + * locating Cantera's data directory (for `gri30.yaml` and friends); + * `include`-ing the generated binding files. +""" +module LibCantera + +export libcantera + +""" + find_library() -> String + +Locate the Cantera shared library. Search order: + + 1. `\$CANTERA_LIBRARY_PATH` (a directory or a full path to the library); + 2. `\$CONDA_PREFIX/lib` for a conda-installed Cantera; + 3. a sibling Cantera build tree (`../../../build/lib`); + 4. the system loader's default search path (`libcantera`). + +Throws an informative error if nothing is found. +""" +function find_library() + # Installed names first, then the in-tree build-tree name (libcantera_shared). + names = ["libcantera.so", "libcantera.dylib", "libcantera.dll", + "libcantera.so.3", "libcantera_shared.so", "libcantera_shared.dylib", + "libcantera_shared.dll", "cantera"] + + function probe(dir) + isdir(dir) || return nothing + for n in names + p = joinpath(dir, n) + isfile(p) && return p + end + # also accept versioned sonames like libcantera.so.3.2.0 + for f in readdir(dir) + startswith(f, "libcantera") && occursin(".so", f) && return joinpath(dir, f) + startswith(f, "libcantera") && endswith(f, ".dylib") && + return joinpath(dir, f) + end + return nothing + end + + if haskey(ENV, "CANTERA_LIBRARY_PATH") + p = ENV["CANTERA_LIBRARY_PATH"] + isfile(p) && return p + hit = probe(p) + hit === nothing || return hit + end + + dirs = String[] + haskey(ENV, "CONDA_PREFIX") && push!(dirs, joinpath(ENV["CONDA_PREFIX"], "lib")) + push!(dirs, normpath(joinpath(@__DIR__, "..", "..", "..", "build", "lib"))) + append!(dirs, ["/usr/local/lib", "/usr/lib"]) + for d in dirs + hit = probe(d) + hit === nothing || return hit + end + + # Fall back to the loader's search path; dlopen will error if truly missing. + return "libcantera" +end + +# Resolved at runtime in `__init__` (below), NOT baked into the precompile image +# — otherwise the discovered path would be frozen at precompile time and later +# changes to CANTERA_LIBRARY_PATH would be silently ignored. Generated `ccall` +# wrappers dereference this Ref as `libcantera[]`. +const libcantera = Ref{String}("libcantera") + +function __init__() + libcantera[] = find_library() + return nothing +end + +# ---- generated low-level bindings ------------------------------------------- +include("generated/_manifest.jl") +for f in GENERATED_FILES + include(joinpath("generated", f)) +end + +""" + default_data_directory() -> Union{String,Nothing} + +Best-effort discovery of this repository's `data` directory (containing +`gri30.yaml`), for the case of running against a source checkout rather than +an installed Cantera. `CANTERA_DATA` and the conda `share/cantera/data` +directory are already registered automatically by Cantera's +`Application::setDefaultDirectories` and don't need to be added here. +""" +function default_data_directory() + d = normpath(joinpath(@__DIR__, "..", "..", "..", "data")) + isdir(d) && isfile(joinpath(d, "gri30.yaml")) && return d + return nothing +end + +end # module LibCantera diff --git a/interfaces/julia/src/connectors.jl b/interfaces/julia/src/connectors.jl new file mode 100644 index 0000000000..f6ddab10dc --- /dev/null +++ b/interfaces/julia/src/connectors.jl @@ -0,0 +1,246 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Reactor connectors: walls and flow devices. +# +# A connector links two reactors in a network. Each Julia wrapper stores the +# CLib connector handle and keeps references to both reactors alive for its +# lifetime. Connectors are freed with `connector_del` via a finalizer. + +""" + Connector + +Abstract supertype of objects that link two reactors in a network +([`Wall`](@ref), [`MassFlowController`](@ref), [`Valve`](@ref), +[`PressureController`](@ref)). +""" +abstract type Connector <: CanteraObject end + +""" + FlowDevice <: Connector + +Abstract supertype of flow devices that move mass between reactors. +""" +abstract type FlowDevice <: Connector end + +# Resolve a Func1-like object (or raw handle) to its CLib handle. +_func_handle(f::Integer) = Int32(f) +_func_handle(f) = Int32(f.handle) + +# ---- concrete types --------------------------------------------------------- + +for T in (:Wall, :MassFlowController, :Valve, :PressureController) + super = T === :Wall ? :Connector : :FlowDevice + @eval begin + mutable struct $T <: $super + handle::Int32 + r0::Reactor # upstream / left + r1::Reactor # downstream / right + closed::Bool + end + end +end + +function _new_connector(::Type{T}, model, left::Reactor, right::Reactor, + name) where {T<:Connector} + h = check(LibCantera.connector_new(model, left.handle, right.handle, name)) + c = T(h, left, right, false) + finalizer(close!, c) + return c +end + +function close!(c::Connector) + c.closed && return c + c.closed = true + try + LibCantera.connector_del(c.handle) + catch + end + return c +end + +"Connector name." +name(c::Connector) = get_string((n, b) -> LibCantera.connector_name(c.handle, n, b)) + +"Connector type string." +connector_type(c::Connector) = + get_string((n, b) -> LibCantera.connector_type(c.handle, n, b)) + +"Set the connector name." +function set_name!(c::Connector, nm::AbstractString) + check(LibCantera.connector_setName(c.handle, nm)) + return c +end + +# ---- walls ------------------------------------------------------------------ + +""" + Wall(left, right; A=1.0, U=0.0, K=0.0, Q=nothing, velocity=nothing, + expansion_rate_coeff=nothing, emissivity=nothing, name="") + +A wall separating reactors `left` and `right`. `A` is the wall area [m^2], `U` +the heat-transfer coefficient [W/m^2/K], `K` the expansion-rate coefficient +[m/s/Pa]. `Q` (heat flux) and `velocity` may be `Func1` objects; `emissivity` +enables radiative transfer. +""" +function Wall(left::Reactor, right::Reactor; A=1.0, U=0.0, K=0.0, Q=nothing, + velocity=nothing, expansion_rate_coeff=nothing, emissivity=nothing, + name::AbstractString="") + w = _new_connector(Wall, "Wall", left, right, name) + set_area!(w, A) + U != 0 && set_heat_transfer_coeff!(w, U) + K != 0 && set_expansion_rate_coeff!(w, K) + expansion_rate_coeff !== nothing && + set_expansion_rate_coeff!(w, expansion_rate_coeff) + emissivity !== nothing && set_emissivity!(w, emissivity) + Q !== nothing && set_heat_flux!(w, Q) + velocity !== nothing && set_velocity!(w, velocity) + return w +end + +"Rate of volumetric expansion of the left reactor [m^3/s]." +expansion_rate(w::Wall) = checkd(LibCantera.wall_expansionRate(w.handle)) + +"Rate of heat transfer through the wall (left to right) [W]." +heat_rate(w::Wall) = checkd(LibCantera.wall_heatRate(w.handle)) + +"Wall area [m^2]." +area(w::Wall) = checkd(LibCantera.wall_area(w.handle)) + +"Set the wall area [m^2]." +function set_area!(w::Wall, a) + check(LibCantera.wall_setArea(w.handle, Float64(a))) + return w +end + +"Set the wall heat-transfer coefficient [W/m^2/K]." +function set_heat_transfer_coeff!(w::Wall, U) + check(LibCantera.wall_setHeatTransferCoeff(w.handle, Float64(U))) + return w +end + +"Set the wall thermal resistance [m^2*K/W]." +function set_thermal_resistance!(w::Wall, Rth) + check(LibCantera.wall_setThermalResistance(w.handle, Float64(Rth))) + return w +end + +"Set the wall expansion-rate coefficient [m/s/Pa]." +function set_expansion_rate_coeff!(w::Wall, k) + check(LibCantera.wall_setExpansionRateCoeff(w.handle, Float64(k))) + return w +end + +"Set the wall emissivity for radiative heat transfer (0..1)." +function set_emissivity!(w::Wall, epsilon) + check(LibCantera.wall_setEmissivity(w.handle, Float64(epsilon))) + return w +end + +"Set the wall heat flux as a time-dependent `Func1` (or handle) [W/m^2]." +function set_heat_flux!(w::Wall, q) + check(LibCantera.wall_setHeatFlux(w.handle, _func_handle(q))) + return w +end + +"Set the wall velocity as a time-dependent `Func1` (or handle) [m/s]." +function set_velocity!(w::Wall, f) + check(LibCantera.wall_setVelocity(w.handle, _func_handle(f))) + return w +end + +# ---- flow devices ----------------------------------------------------------- + +"Mass flow rate through the flow device [kg/s]." +mass_flow_rate(d::FlowDevice) = checkd(LibCantera.flowdev_massFlowRate(d.handle)) + +"Device coefficient of the flow device (meaning depends on the device type)." +device_coefficient(d::FlowDevice) = + checkd(LibCantera.flowdev_deviceCoefficient(d.handle)) + +"Set the device coefficient of the flow device." +function set_device_coefficient!(d::FlowDevice, c) + check(LibCantera.flowdev_setDeviceCoefficient(d.handle, Float64(c))) + return d +end + +"Set the pressure function of the flow device as a `Func1` (or handle)." +function set_pressure_function!(d::FlowDevice, f) + check(LibCantera.flowdev_setPressureFunction(d.handle, _func_handle(f))) + return d +end + +"Set the time function of the flow device as a `Func1` (or handle)." +function set_time_function!(d::FlowDevice, g) + check(LibCantera.flowdev_setTimeFunction(d.handle, _func_handle(g))) + return d +end + +"Set the primary flow device (used by [`PressureController`](@ref))." +function set_primary!(d::FlowDevice, primary::FlowDevice) + check(LibCantera.flowdev_setPrimary(d.handle, primary.handle)) + return d +end + +""" + MassFlowController(upstream, downstream; mdot=0.0, name="") + +A flow device that maintains a specified mass flow rate `mdot` [kg/s] from +`upstream` to `downstream`, independent of the pressure difference. +""" +function MassFlowController(upstream::Reactor, downstream::Reactor; mdot=0.0, + name::AbstractString="") + d = _new_connector(MassFlowController, "MassFlowController", upstream, downstream, + name) + mdot != 0 && set_mass_flow_rate!(d, mdot) + return d +end + +""" + set_mass_flow_rate!(mfc, mdot) + +Set the (constant) mass flow rate of a [`MassFlowController`](@ref) [kg/s]. +""" +set_mass_flow_rate!(d::MassFlowController, mdot) = set_device_coefficient!(d, mdot) + +""" + Valve(upstream, downstream; K=0.0, name="") + +A flow device whose mass flow rate is proportional to the pressure difference, +`mdot = K * (P_upstream - P_downstream)`, with valve coefficient `K`. +""" +function Valve(upstream::Reactor, downstream::Reactor; K=0.0, name::AbstractString="") + d = _new_connector(Valve, "Valve", upstream, downstream, name) + K != 0 && set_device_coefficient!(d, K) + return d +end + +""" + PressureController(upstream, downstream; primary=nothing, K=0.0, name="") + +A flow device that regulates the pressure difference across it relative to a +`primary` flow device: `mdot = primary.mdot + K * (P_upstream - P_downstream)`. +""" +function PressureController(upstream::Reactor, downstream::Reactor; primary=nothing, + K=0.0, name::AbstractString="") + d = _new_connector(PressureController, "PressureController", upstream, downstream, + name) + K != 0 && set_device_coefficient!(d, K) + primary !== nothing && set_primary!(d, primary) + return d +end + +# ---- display ---------------------------------------------------------------- + +for T in (:Wall, :MassFlowController, :Valve, :PressureController) + @eval Base.show(io::IO, c::$T) = + print(io, c.closed ? string($(QuoteNode(T)), "()") : + string($(QuoteNode(T)), "(\"", name(c), "\")")) +end + +export Connector, FlowDevice, Wall, MassFlowController, Valve, PressureController, + device_coefficient, set_device_coefficient!, set_pressure_function!, + set_time_function!, set_primary!, expansion_rate, heat_rate, + set_heat_transfer_coeff!, set_thermal_resistance!, set_expansion_rate_coeff!, + set_emissivity!, set_heat_flux!, set_velocity!, connector_type, set_name!, + mass_flow_rate, set_mass_flow_rate!, area, set_area! diff --git a/interfaces/julia/src/errors.jl b/interfaces/julia/src/errors.jl new file mode 100644 index 0000000000..9d4aab2904 --- /dev/null +++ b/interfaces/julia/src/errors.jl @@ -0,0 +1,63 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Error handling for the Cantera CLib boundary. +# +# CLib functions signal failure with sentinel return values: +# * `int32_t` functions return -1 (or ERR = -999); +# * `double` functions return DERR = -999.999. +# When such a sentinel is seen, the human-readable message is retrieved with +# `ct_getCanteraError` and raised as a `CanteraError`. + +const _ERR = Int32(-999) # ERR from clib_defs.h +const _DERR = -999.999 # DERR from clib_defs.h + +""" + CanteraError(msg) + +Exception raised when an underlying Cantera CLib call reports an error. The +message is the exception text produced by the C++ Cantera core. +""" +struct CanteraError <: Exception + msg::String +end + +Base.showerror(io::IO, e::CanteraError) = print(io, "CanteraError: ", e.msg) + +export CanteraError + +"Retrieve and clear the last error message stored by the Cantera CLib." +function last_cantera_error() + buflen = LibCantera.ct_getCanteraError(Int32(0), Ptr{UInt8}(C_NULL)) + buflen <= 0 && return "unknown Cantera error (no message available)" + buf = Vector{UInt8}(undef, buflen) + LibCantera.ct_getCanteraError(Int32(buflen), pointer(buf)) + return unsafe_string(pointer(buf)) +end + +""" + check(code::Integer) -> Int32 + +Verify the return code of an `int32_t` CLib call and raise [`CanteraError`](@ref) +on the error sentinel. Returns the code unchanged on success so it can be used +inline, e.g. `n = check(LibCantera.thermo_nSpecies(h))`. +""" +function check(code::Integer) + # All successful `int32_t` CLib returns are non-negative (handles, counts, + # indices, string lengths, or the 0 status). Any negative value signals an + # exception. Callers that must distinguish a "not found" npos result (also + # negative) do so *before* calling `check` (see `species_index`). + code < 0 && throw(CanteraError(last_cantera_error())) + return Int32(code) +end + +""" + checkd(value::Float64) -> Float64 + +Verify the return value of a `double` CLib call and raise [`CanteraError`](@ref) +on the `DERR` sentinel. +""" +function checkd(value::Float64) + value == _DERR && throw(CanteraError(last_cantera_error())) + return value +end diff --git a/interfaces/julia/src/func1.jl b/interfaces/julia/src/func1.jl new file mode 100644 index 0000000000..f5dc954ee7 --- /dev/null +++ b/interfaces/julia/src/func1.jl @@ -0,0 +1,132 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Func1 -- functors of a single variable. +# +# Cantera's `Func1` objects wrap C++ function objects (sin, cos, polynomials, +# compound sum/product functions, ...) that can be evaluated, differentiated +# symbolically and combined arithmetically. Each object is referenced by an +# `Int32` handle into the `func1` cabinet. + +""" + Func1 + +A wrapper around a Cantera `Func1` function object. A `Func1` is +callable: `f(t)` evaluates it at `t`. Instances combine with `+`, `-`, `*` and +`/` to build compound functions and support symbolic [`derivative`](@ref). + +# Examples +```julia +f = Func1("sin", 2.0) # t -> sin(2t) +f(pi/4) # == sin(pi/2) +g = Func1("polynomial3", [1.0, 2.0, 3.0, 4.0]) # t^3 + 2t^2 + 3t + 4 +h = Func1("constant", 2.0) + Func1("constant", 3.0) +``` +""" +mutable struct Func1 <: CanteraObject + handle::Int32 + closed::Bool + + function Func1(handle::Integer) + f = new(Int32(handle), false) + finalizer(close!, f) + return f + end +end + +""" + Func1(type::String, coeff::Real=1.0) + +Construct a basic functor of the given `type` (e.g. `"sin"`, `"cos"`, `"exp"`, +`"log"`, `"pow"`, `"constant"`). For `"sin"` with coefficient `w`, evaluating +at `t` returns `sin(w*t)`. +""" +Func1(type::AbstractString, coeff::Real=1.0) = + Func1(check(LibCantera.func1_newBasic(type, Float64(coeff)))) + +""" + Func1(type::String, coeffs::AbstractVector) + +Construct an advanced functor parametrized by an array of coefficients, e.g. +`"polynomial3"` (coefficients from highest to lowest degree) or `"Fourier"`. +""" +function Func1(type::AbstractString, coeffs::AbstractVector{<:Real}) + c = as_f64(coeffs) + return Func1(check( + LibCantera.func1_newAdvanced(type, Int32(length(c)), pointer(c)))) +end + +""" + constant_function(c) -> Func1 + +Convenience constructor for the constant functor `t -> c`. +""" +constant_function(c::Real) = Func1("constant", c) + +""" + close!(f::Func1) + +Release the Cantera-side object backing `f`. Idempotent; called automatically +by the finalizer. +""" +function close!(f::Func1) + f.closed && return f + f.closed = true + try + LibCantera.func1_del(f.handle) + catch + end + return f +end + +""" + evaluate(f::Func1, t) -> Float64 + +Evaluate the functor at `t`. +""" +evaluate(f::Func1, t::Real) = checkd(LibCantera.func1_eval(f.handle, Float64(t))) + +(f::Func1)(t::Real) = evaluate(f, t) + +""" + derivative(f::Func1) -> Func1 + +Return a new `Func1` representing the symbolic derivative of `f`. +""" +derivative(f::Func1) = Func1(check(LibCantera.func1_derivative(f.handle))) + +""" + func_type(f::Func1) -> String + +The functor's type string (e.g. `"sin"`, `"sum"`). +""" +func_type(f::Func1) = get_string((n, b) -> LibCantera.func1_type(f.handle, n, b)) + +""" + write(f::Func1, arg="t") -> String + +Return a string (LaTeX-style) representation of `f` using `arg` as the +independent-variable symbol. +""" +write(f::Func1, arg::AbstractString="t") = + get_string((n, b) -> LibCantera.func1_write(f.handle, arg, n, b)) + +# Arithmetic combinations build compound functors. +Base.:+(f::Func1, g::Func1) = + Func1(check(LibCantera.func1_newSumFunction(f.handle, g.handle))) +Base.:-(f::Func1, g::Func1) = + Func1(check(LibCantera.func1_newDiffFunction(f.handle, g.handle))) +Base.:*(f::Func1, g::Func1) = + Func1(check(LibCantera.func1_newProdFunction(f.handle, g.handle))) +Base.:/(f::Func1, g::Func1) = + Func1(check(LibCantera.func1_newRatioFunction(f.handle, g.handle))) + +function Base.show(io::IO, f::Func1) + if f.closed + print(io, "Func1()") + else + print(io, "Func1(\"", func_type(f), "\": ", write(f), ")") + end +end + +export Func1, evaluate, derivative, func_type, constant_function diff --git a/interfaces/julia/src/handles.jl b/interfaces/julia/src/handles.jl new file mode 100644 index 0000000000..c551ce36b9 --- /dev/null +++ b/interfaces/julia/src/handles.jl @@ -0,0 +1,88 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Wrapper types and low-level handle helpers. +# +# Every Cantera CLib object is referenced by an `Int32` handle into a C++-side +# "cabinet". The Julia wrapper types below store that handle and are the only +# things user code touches; raw handles are never part of the public API. + +""" + CanteraObject + +Abstract supertype of all Julia wrappers around Cantera CLib objects. +""" +abstract type CanteraObject end + +""" + ThermoPhase + +Thermodynamic-phase view of a [`Solution`](@ref). Wraps a `ThermoPhase` CLib +handle and holds a reference to the owning `Solution`. +""" +mutable struct ThermoPhase{P} <: CanteraObject + handle::Int32 + parent::P +end + +""" + Kinetics + +Kinetics view of a [`Solution`](@ref). Wraps a `Kinetics` CLib handle and holds +a reference to the owning `Solution`. +""" +mutable struct Kinetics{P} <: CanteraObject + handle::Int32 + parent::P +end + +""" + Transport + +Transport view of a [`Solution`](@ref). Wraps a `Transport` CLib handle and +holds a reference to the owning `Solution`. +""" +mutable struct Transport{P} <: CanteraObject + handle::Int32 + parent::P +end + +"`handle(obj)` returns the raw CLib handle backing a wrapper (internal)." +handle(o::CanteraObject) = o.handle + +export CanteraObject, ThermoPhase, Kinetics, Transport + +# ---- generic string / array marshalling ------------------------------------ + +""" + get_string(f) -> String + +Call a CLib string getter through the closure `f(bufLen::Int32, buf::Ptr{UInt8})` +using the size-query protocol: first with a null buffer to obtain the required +length, then again with a correctly sized buffer. +""" +function get_string(f) + n = check(f(Int32(0), Ptr{UInt8}(C_NULL))) + n <= 0 && return "" + buf = Vector{UInt8}(undef, n) + check(f(Int32(n), pointer(buf))) + return unsafe_string(pointer(buf)) +end + +""" + get_array!(out::Vector{Float64}, f) -> out + +Fill `out` in place through the closure `f(len::Int32, ptr::Ptr{Float64})` and +return it. Used by the performance-oriented `foo!(out, gas)` variants. +""" +function get_array!(out::Vector{Float64}, f) + check(f(Int32(length(out)), pointer(out))) + return out +end + +"`get_array(n, f)` allocates a `Vector{Float64}` of length `n` and fills it." +get_array(n::Integer, f) = get_array!(Vector{Float64}(undef, n), f) + +"Convert an arbitrary real vector to a contiguous `Vector{Float64}` at the boundary." +as_f64(x::Vector{Float64}) = x +as_f64(x::AbstractVector{<:Real}) = convert(Vector{Float64}, x) diff --git a/interfaces/julia/src/kinetics.jl b/interfaces/julia/src/kinetics.jl new file mode 100644 index 0000000000..350751587f --- /dev/null +++ b/interfaces/julia/src/kinetics.jl @@ -0,0 +1,266 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Kinetics accessors. All functions accept a `Solution` or a `Kinetics`. + +"Number of reactions in the mechanism." +n_reactions(g::KineticsLike) = + Int(check(LibCantera.kin_nReactions(_kinetics_handle(g)))) + +"Total number of species across all phases in the Kinetics manager." +n_total_species(g::KineticsLike) = + Int(check(LibCantera.kin_nTotalSpecies(_kinetics_handle(g)))) + +"Reaction equation string for reaction `i` (1-based)." +function reaction_equation(g::KineticsLike, i::Integer) + r = check(LibCantera.kin_reaction(_kinetics_handle(g), Int32(i - 1))) + eq = get_string((n, b) -> LibCantera.rxn_equation(r, n, b)) + LibCantera.rxn_del(r) + return eq +end + +"Vector of all reaction equation strings." +reaction_equations(g::KineticsLike) = + [reaction_equation(g, i) for i in 1:n_reactions(g)] + +"Whether reaction `i` (1-based) is reversible." +is_reversible(g::KineticsLike, i::Integer) = + check(LibCantera.kin_isReversible(_kinetics_handle(g), Int32(i - 1))) != 0 + +# Length helpers: rate/ROP arrays are per-reaction; production arrays per-species. +_nrxn(g) = n_reactions(g) +_nsp(g) = n_total_species(g) + +# Per-reaction kinetic arrays. +for (jl, c) in ( + (:forward_rates_of_progress, :kin_getFwdRatesOfProgress), + (:reverse_rates_of_progress, :kin_getRevRatesOfProgress), + (:net_rates_of_progress, :kin_getNetRatesOfProgress), + (:equilibrium_constants, :kin_getEquilibriumConstants), + (:forward_rate_constants, :kin_getFwdRateConstants), + (:delta_enthalpy, :kin_getDeltaEnthalpy), + (:delta_gibbs, :kin_getDeltaGibbs), + (:delta_entropy, :kin_getDeltaEntropy), + ) + bang = Symbol(jl, :!) + @eval begin + $jl(g::KineticsLike) = + get_array(_nrxn(g), (n, b) -> LibCantera.$c(_kinetics_handle(g), n, b)) + $bang(out, g::KineticsLike) = + get_array!(out, (n, b) -> LibCantera.$c(_kinetics_handle(g), n, b)) + end +end + +""" + reverse_rate_constants(gas; do_irreversible=false) + +Reverse rate constants for all reactions. +""" +function reverse_rate_constants(g::KineticsLike; do_irreversible::Bool=false) + h = _kinetics_handle(g) + return get_array(_nrxn(g), + (n, b) -> LibCantera.kin_getRevRateConstants(h, n, b, Int32(do_irreversible))) +end + +# Per-species production arrays. +for (jl, c) in ( + (:net_production_rates, :kin_getNetProductionRates), + (:creation_rates, :kin_getCreationRates), + (:destruction_rates, :kin_getDestructionRates), + ) + bang = Symbol(jl, :!) + @eval begin + $jl(g::KineticsLike) = + get_array(_nsp(g), (n, b) -> LibCantera.$c(_kinetics_handle(g), n, b)) + $bang(out, g::KineticsLike) = + get_array!(out, (n, b) -> LibCantera.$c(_kinetics_handle(g), n, b)) + end +end + +# ---- analytic kinetic derivatives ------------------------------------------- +# These wrap Cantera's native analytic derivative getters, exposed through the +# CLib via recipes added to interfaces/sourcegen/.../headers/ctkin.yaml (the +# `kin_get*_dd{T,P,C}` functions). Each is an array getter with the standard +# `(len, ptr)` array convention, so they reuse `get_array`/`get_array!`. +# +# Derivatives are taken at constant volume/composition for `_ddT` and `_ddP`, +# and with respect to molar concentration for `_ddC`, following Cantera's C++ +# convention. The composition Jacobian (`_ddX`) returns a *sparse* matrix in +# C++; it is wrapped separately below, densified into a flat column-major +# buffer to fit the array-out CLib convention. + +for (jl, c, len) in ( + (:net_production_rates_ddT, :kin_getNetProductionRates_ddT, :_nsp), + (:net_production_rates_ddP, :kin_getNetProductionRates_ddP, :_nsp), + (:net_production_rates_ddC, :kin_getNetProductionRates_ddC, :_nsp), + (:creation_rates_ddT, :kin_getCreationRates_ddT, :_nsp), + (:creation_rates_ddP, :kin_getCreationRates_ddP, :_nsp), + (:creation_rates_ddC, :kin_getCreationRates_ddC, :_nsp), + (:destruction_rates_ddT, :kin_getDestructionRates_ddT, :_nsp), + (:destruction_rates_ddP, :kin_getDestructionRates_ddP, :_nsp), + (:destruction_rates_ddC, :kin_getDestructionRates_ddC, :_nsp), + (:net_rates_of_progress_ddT, :kin_getNetRatesOfProgress_ddT, :_nrxn), + (:net_rates_of_progress_ddP, :kin_getNetRatesOfProgress_ddP, :_nrxn), + (:net_rates_of_progress_ddC, :kin_getNetRatesOfProgress_ddC, :_nrxn), + ) + bang = Symbol(jl, :!) + @eval begin + $jl(g::KineticsLike) = + get_array($len(g), (n, b) -> LibCantera.$c(_kinetics_handle(g), n, b)) + $bang(out, g::KineticsLike) = + get_array!(out, (n, b) -> LibCantera.$c(_kinetics_handle(g), n, b)) + end +end + +# ---- composition Jacobians (_ddX) ------------------------------------------- +# The C++ `_ddX` getters return sparse matrices; the CLib densifies them into a +# flat column-major buffer (see interfaces/sourcegen/.../ctkin.yaml), which we +# reshape into a Julia matrix. + +""" + net_production_rates_ddX(gas) -> Matrix{Float64} + +Dense composition Jacobian d(wdot_k)/dX_j of the net production rates with +respect to species mole fractions, size nSpecies x nSpecies. +""" +function net_production_rates_ddX(g::KineticsLike) + n = _nsp(g) + buf = get_array(n * n, + (len, b) -> LibCantera.kin_getNetProductionRates_ddX( + _kinetics_handle(g), len, b)) + return reshape(buf, n, n) +end + +""" + net_rates_of_progress_ddX(gas) -> Matrix{Float64} + +Dense composition Jacobian d(q_i)/dX_j of the net rates of progress, size +nReactions x nSpecies. +""" +function net_rates_of_progress_ddX(g::KineticsLike) + nr = _nrxn(g); ns = _nsp(g) + buf = get_array(nr * ns, + (len, b) -> LibCantera.kin_getNetRatesOfProgress_ddX( + _kinetics_handle(g), len, b)) + return reshape(buf, nr, ns) +end + +""" + heat_release_rate(gas) -> Float64 + +Volumetric heat release rate [W/m^3], `-Σ_k h_k · ẇ_k`, where `h_k` are the +partial molar enthalpies and `ẇ_k` the net production rates. Matches Python's +`gas.heat_release_rate` (the CLib exposes no direct getter). +""" +heat_release_rate(g::Solution) = + -sum(partial_molar_enthalpies(g) .* net_production_rates(g)) + +# ---- stoichiometry ---------------------------------------------------------- + +""" +Stoichiometric coefficient of reactant species `k` (1-based) in reaction `i` +(1-based). +""" +reactant_stoich_coeff(g::KineticsLike, k::Integer, i::Integer) = + checkd(LibCantera.kin_reactantStoichCoeff( + _kinetics_handle(g), Int32(k - 1), Int32(i - 1))) + +""" +Stoichiometric coefficient of product species `k` (1-based) in reaction `i` +(1-based). +""" +product_stoich_coeff(g::KineticsLike, k::Integer, i::Integer) = + checkd(LibCantera.kin_productStoichCoeff( + _kinetics_handle(g), Int32(k - 1), Int32(i - 1))) + +"`nSpecies × nReactions` matrix of reactant stoichiometric coefficients." +function reactant_stoich_coeffs(g::KineticsLike) + ns = _nsp(g); nr = _nrxn(g) + [reactant_stoich_coeff(g, k, i) for k in 1:ns, i in 1:nr] +end + +"`nSpecies × nReactions` matrix of product stoichiometric coefficients." +function product_stoich_coeffs(g::KineticsLike) + ns = _nsp(g); nr = _nrxn(g) + [product_stoich_coeff(g, k, i) for k in 1:ns, i in 1:nr] +end + +# ---- reaction multipliers --------------------------------------------------- + +"Rate multiplier for reaction `i` (1-based)." +multiplier(g::KineticsLike, i::Integer) = + checkd(LibCantera.kin_multiplier(_kinetics_handle(g), Int32(i - 1))) + +"Set the rate multiplier for reaction `i` (1-based)." +function set_multiplier!(g::KineticsLike, i::Integer, v) + check(LibCantera.kin_setMultiplier(_kinetics_handle(g), Int32(i - 1), Float64(v))) + return g +end + +# ---- standard-state reaction deltas ----------------------------------------- + +"Standard-state enthalpy change for each reaction [J/kmol]." +delta_standard_enthalpy(g::KineticsLike) = + get_array(_nrxn(g), + (n, b) -> LibCantera.kin_getDeltaSSEnthalpy(_kinetics_handle(g), n, b)) + +"Standard-state Gibbs free energy change for each reaction [J/kmol]." +delta_standard_gibbs(g::KineticsLike) = + get_array(_nrxn(g), + (n, b) -> LibCantera.kin_getDeltaSSGibbs(_kinetics_handle(g), n, b)) + +"Standard-state entropy change for each reaction [J/kmol/K]." +delta_standard_entropy(g::KineticsLike) = + get_array(_nrxn(g), + (n, b) -> LibCantera.kin_getDeltaSSEntropy(_kinetics_handle(g), n, b)) + +# ---- forward/reverse rate-of-progress derivatives --------------------------- + +for (jl, c) in ( + (:forward_rates_of_progress_ddT, :kin_getFwdRatesOfProgress_ddT), + (:forward_rates_of_progress_ddP, :kin_getFwdRatesOfProgress_ddP), + (:forward_rates_of_progress_ddC, :kin_getFwdRatesOfProgress_ddC), + (:reverse_rates_of_progress_ddT, :kin_getRevRatesOfProgress_ddT), + (:reverse_rates_of_progress_ddP, :kin_getRevRatesOfProgress_ddP), + (:reverse_rates_of_progress_ddC, :kin_getRevRatesOfProgress_ddC), + ) + @eval $jl(g::KineticsLike) = + get_array(_nrxn(g), (n, b) -> LibCantera.$c(_kinetics_handle(g), n, b)) +end + +""" + heat_production_rates(gas) -> Vector{Float64} + +Per-reaction volumetric heat production rates [W/m^3], `-q_i · Δh_i` (net rate +of progress times reaction enthalpy). Their sum equals [`heat_release_rate`](@ref). +""" +heat_production_rates(g::KineticsLike) = + -net_rates_of_progress(g) .* delta_enthalpy(g) + +export n_reactions, n_total_species, reaction_equation, reaction_equations, + is_reversible, + forward_rate_constants, reverse_rate_constants, equilibrium_constants, + forward_rates_of_progress, reverse_rates_of_progress, + net_rates_of_progress, net_production_rates, net_production_rates!, + creation_rates, creation_rates!, destruction_rates, destruction_rates!, + forward_rates_of_progress!, reverse_rates_of_progress!, + net_rates_of_progress!, forward_rate_constants!, equilibrium_constants!, + delta_enthalpy, delta_gibbs, delta_entropy, + net_production_rates_ddT, net_production_rates_ddP, net_production_rates_ddC, + creation_rates_ddT, creation_rates_ddP, creation_rates_ddC, + destruction_rates_ddT, destruction_rates_ddP, destruction_rates_ddC, + net_rates_of_progress_ddT, net_rates_of_progress_ddP, net_rates_of_progress_ddC, + net_production_rates_ddX, net_rates_of_progress_ddX, + net_production_rates_ddT!, net_production_rates_ddP!, net_production_rates_ddC!, + creation_rates_ddT!, creation_rates_ddP!, creation_rates_ddC!, + destruction_rates_ddT!, destruction_rates_ddP!, destruction_rates_ddC!, + net_rates_of_progress_ddT!, net_rates_of_progress_ddP!, + net_rates_of_progress_ddC!, + heat_release_rate, heat_production_rates, + reactant_stoich_coeff, product_stoich_coeff, + reactant_stoich_coeffs, product_stoich_coeffs, + multiplier, set_multiplier!, + delta_standard_enthalpy, delta_standard_gibbs, delta_standard_entropy, + forward_rates_of_progress_ddT, forward_rates_of_progress_ddP, + forward_rates_of_progress_ddC, reverse_rates_of_progress_ddT, + reverse_rates_of_progress_ddP, reverse_rates_of_progress_ddC diff --git a/interfaces/julia/src/multiphase.jl b/interfaces/julia/src/multiphase.jl new file mode 100644 index 0000000000..223e0a214a --- /dev/null +++ b/interfaces/julia/src/multiphase.jl @@ -0,0 +1,250 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# MultiPhase mixtures. +# +# A `MultiPhase` holds one or more phases (each a `Solution`) together with the +# number of moles of each, and computes multi-phase equilibrium and mixture +# thermodynamic properties. It mirrors Python's `cantera.Mixture`. + +""" + MultiPhase() + MultiPhase(phases::AbstractVector{<:Pair{Solution,<:Real}}) + +A mixture of one or more phases. Construct empty and add phases with +[`add_phase!`](@ref) followed by [`init!`](@ref), or pass a vector of +`solution => moles` pairs to add and initialize in one step. + +```julia +gas = Solution("gri30.yaml") +mix = MultiPhase([gas => 1.0]) +set_temperature!(mix, 1200.0); set_pressure!(mix, one_atm) +equilibrate!(mix, "TP") +``` +""" +mutable struct MultiPhase <: CanteraObject + handle::Int32 + phases::Vector{Solution} # keep added Solutions alive + initialized::Bool + closed::Bool + + function MultiPhase() + h = check(LibCantera.mix_new()) + mp = new(h, Solution[], false, false) + finalizer(close!, mp) + return mp + end +end + +function MultiPhase(phases::AbstractVector{<:Pair{Solution,<:Real}}) + mp = MultiPhase() + for (sol, moles) in phases + add_phase!(mp, sol, moles) + end + init!(mp) + return mp +end + +""" + close!(mp::MultiPhase) + +Release the Cantera-side object backing `mp`. Idempotent. +""" +function close!(mp::MultiPhase) + mp.closed && return mp + mp.closed = true + try + LibCantera.mix_del(mp.handle) + catch + end + return mp +end + +""" + add_phase!(mp::MultiPhase, solution::Solution, moles::Real) + +Add `solution` to the mixture with the given number of `moles`. Call +[`init!`](@ref) after all phases have been added. +""" +function add_phase!(mp::MultiPhase, sol::Solution, moles::Real) + check(LibCantera.mix_addPhase(mp.handle, sol.handle, Float64(moles))) + push!(mp.phases, sol) + return mp +end + +""" + init!(mp::MultiPhase) + +Finalize the mixture after all phases have been added. Required before any +property or equilibrium calculation. +""" +function init!(mp::MultiPhase) + check(LibCantera.mix_init(mp.handle)) + mp.initialized = true + return mp +end + +""" + update_phases!(mp::MultiPhase) + +Synchronize the constituent phase objects with the mixture's current state. +""" +function update_phases!(mp::MultiPhase) + check(LibCantera.mix_updatePhases(mp.handle)) + return mp +end + +# ---- sizes ------------------------------------------------------------------ + +"Number of phases in the mixture." +n_phases(mp::MultiPhase) = Int(check(LibCantera.mix_nPhases(mp.handle))) + +"Total number of species across all phases." +n_species(mp::MultiPhase) = Int(check(LibCantera.mix_nSpecies(mp.handle))) + +"Number of elements in the mixture." +n_elements(mp::MultiPhase) = Int(check(LibCantera.mix_nElements(mp.handle))) + +""" + element_index(mp::MultiPhase, name) -> Int + +1-based global index of element `name`, or `0` if not present. +""" +function element_index(mp::MultiPhase, name::AbstractString) + idx = LibCantera.mix_elementIndex(mp.handle, name) + idx == -1 && return 0 + return Int(idx) + 1 +end + +# ---- scalar state ----------------------------------------------------------- + +"Temperature [K]." +temperature(mp::MultiPhase) = checkd(LibCantera.mix_temperature(mp.handle)) + +"Set the mixture temperature [K]." +function set_temperature!(mp::MultiPhase, T::Real) + check(LibCantera.mix_setTemperature(mp.handle, Float64(T))) + return mp +end + +"Pressure [Pa]." +pressure(mp::MultiPhase) = checkd(LibCantera.mix_pressure(mp.handle)) + +"Set the mixture pressure [Pa]." +function set_pressure!(mp::MultiPhase, P::Real) + check(LibCantera.mix_setPressure(mp.handle, Float64(P))) + return mp +end + +"Minimum temperature [K] of the mixture's valid range." +min_temp(mp::MultiPhase) = checkd(LibCantera.mix_minTemp(mp.handle)) + +"Maximum temperature [K] of the mixture's valid range." +max_temp(mp::MultiPhase) = checkd(LibCantera.mix_maxTemp(mp.handle)) + +"Total electric charge [C/mol] of the mixture." +charge(mp::MultiPhase) = checkd(LibCantera.mix_charge(mp.handle)) + +# ---- mixture thermodynamic properties --------------------------------------- + +"Total enthalpy [J]." +enthalpy(mp::MultiPhase) = checkd(LibCantera.mix_enthalpy(mp.handle)) + +"Total entropy [J/K]." +entropy(mp::MultiPhase) = checkd(LibCantera.mix_entropy(mp.handle)) + +"Total Gibbs free energy [J]." +gibbs(mp::MultiPhase) = checkd(LibCantera.mix_gibbs(mp.handle)) + +"Total heat capacity at constant pressure [J/K]." +cp(mp::MultiPhase) = checkd(LibCantera.mix_cp(mp.handle)) + +"Total volume [m^3]." +volume(mp::MultiPhase) = checkd(LibCantera.mix_volume(mp.handle)) + +# ---- composition ------------------------------------------------------------ + +"Number of moles in phase `n` (1-based)." +phase_moles(mp::MultiPhase, n::Integer) = + checkd(LibCantera.mix_phaseMoles(mp.handle, Int32(n - 1))) + +"Set the number of moles in phase `n` (1-based)." +function set_phase_moles!(mp::MultiPhase, n::Integer, moles::Real) + check(LibCantera.mix_setPhaseMoles(mp.handle, Int32(n - 1), Float64(moles))) + return mp +end + +"Moles of global species `k` (1-based)." +species_moles(mp::MultiPhase, k::Integer) = + checkd(LibCantera.mix_speciesMoles(mp.handle, Int32(k - 1))) + +"Moles of element `m` (1-based)." +element_moles(mp::MultiPhase, m::Integer) = + checkd(LibCantera.mix_elementMoles(mp.handle, Int32(m - 1))) + +"Mole fraction of global species `k` (1-based)." +mole_fraction(mp::MultiPhase, k::Integer) = + checkd(LibCantera.mix_moleFraction(mp.handle, Int32(k - 1))) + +""" + n_atoms(mp::MultiPhase, k, m) -> Float64 + +Number of atoms of element `m` in global species `k` (both 1-based). +""" +n_atoms(mp::MultiPhase, k::Integer, m::Integer) = + checkd(LibCantera.mix_nAtoms(mp.handle, Int32(k - 1), Int32(m - 1))) + +""" + set_moles!(mp::MultiPhase, moles) + +Set the moles of all global species from a numeric vector, or from a +composition string (e.g. `"CH4:1, O2:2"`). +""" +function set_moles!(mp::MultiPhase, moles::AbstractVector{<:Real}) + n = as_f64(moles) + check(LibCantera.mix_setMoles(mp.handle, Int32(length(n)), pointer(n))) + return mp +end +function set_moles!(mp::MultiPhase, moles::AbstractString) + check(LibCantera.mix_setMolesByName(mp.handle, moles)) + return mp +end + +""" + chemical_potentials(mp::MultiPhase) -> Vector{Float64} + +Chemical potentials [J/kmol] of all global species. +""" +chemical_potentials(mp::MultiPhase) = + get_array(n_species(mp), + (n, b) -> LibCantera.mix_getChemPotentials(mp.handle, n, b)) + +""" + equilibrate!(mp::MultiPhase, XY; solver="auto", rtol=1e-9, max_steps=1000, + max_iter=200, estimate_equil=0) + +Bring the mixture to chemical equilibrium holding the property pair `XY` +(e.g. `"TP"`) fixed. +""" +function equilibrate!(mp::MultiPhase, XY::AbstractString; solver="auto", + rtol=1e-9, max_steps=1000, max_iter=200, estimate_equil=0) + check(LibCantera.mix_equilibrate(mp.handle, XY, solver, Float64(rtol), + Int32(max_steps), Int32(max_iter), + Int32(estimate_equil))) + return mp +end + +function Base.show(io::IO, mp::MultiPhase) + if mp.closed + print(io, "MultiPhase()") + else + print(io, "MultiPhase(", n_phases(mp), " phases, ", n_species(mp), + " species)") + end +end + +# note: `cp` is available as `Cantera.cp` to avoid the clash with `Base.cp` +export MultiPhase, add_phase!, init!, update_phases!, n_phases, + set_temperature!, set_pressure!, min_temp, max_temp, charge, + enthalpy, entropy, gibbs, volume, phase_moles, set_phase_moles!, + species_moles, element_moles, mole_fraction, n_atoms, set_moles! diff --git a/interfaces/julia/src/onedim.jl b/interfaces/julia/src/onedim.jl new file mode 100644 index 0000000000..7f4414767d --- /dev/null +++ b/interfaces/julia/src/onedim.jl @@ -0,0 +1,736 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# One-dimensional flames (Domain1D / Sim1D). +# +# Parity wrapper over the `ctdomain` and `ctonedim` CLib libraries, mirroring +# Python/MATLAB `FreeFlame`. A [`FreeFlame`](@ref) owns three [`Domain1D`](@ref) +# objects (inlet / flow / outlet) plus a `Sim1D` solver handle, and keeps a +# reference to the gas [`Solution`](@ref) alive for their lifetime. + +# --------------------------------------------------------------------------- +# Domain1D: a generic wrapper around a single 1-D domain handle. +# --------------------------------------------------------------------------- + +""" + Domain1D + +Wrapper around a single Cantera 1-D domain (a boundary such as an inlet/outlet, +or a flow region). Domains are normally created and managed by a +[`FreeFlame`](@ref); the accessors here let you inspect and set profiles on an +individual domain. Domains borrow the gas [`Solution`](@ref) handle and are +freed together with the owning simulation. +""" +mutable struct Domain1D <: CanteraObject + handle::Int32 + solution::Solution # keep the phase alive + closed::Bool +end + +"Number of solution points (grid nodes) in the domain." +n_points(d::Domain1D) = Int(check(LibCantera.domain_nPoints(d.handle))) + +"Number of solution components in the domain." +n_components(d::Domain1D) = Int(check(LibCantera.domain_nComponents(d.handle))) + +"Name of component `n` (1-based)." +function component_name(d::Domain1D, n::Integer) + return get_string( + (bl, b) -> LibCantera.domain_componentName(d.handle, Int32(n - 1), bl, b)) +end + +"Vector of all component names in the domain." +component_names(d::Domain1D) = [component_name(d, n) for n in 1:n_components(d)] + +"Domain type string (e.g. \"free-flow\", \"inlet\", \"outlet\")." +domain_type(d::Domain1D) = + get_string((bl, b) -> LibCantera.domain_domainType(d.handle, bl, b)) + +"Grid points of the domain [m]." +grid(d::Domain1D) = + get_array(n_points(d), (bl, b) -> LibCantera.domain_grid(d.handle, bl, b)) + +"Scalar value of `component` in a boundary domain (a single-point domain)." +value(d::Domain1D, component::AbstractString) = + checkd(LibCantera.domain_value(d.handle, component)) + +"Per-point profile of `component` across the domain grid." +function solution_profile(d::Domain1D, component::AbstractString) + np = n_points(d) + return get_array(np, + (n, b) -> LibCantera.domain_getValues(d.handle, component, n, b)) +end + +""" + set_profile!(domain, component, positions, values) + +Set the profile of `component` over normalized positions `positions` (in +`[0, 1]`, spanning the domain) with the given `values`. +""" +function set_profile!(d::Domain1D, component::AbstractString, positions, values) + pos = as_f64(collect(positions)) + val = as_f64(collect(values)) + check(LibCantera.domain_setProfile(d.handle, component, + Int32(length(pos)), pointer(pos), + Int32(length(val)), pointer(val))) + return d +end + +"Set a spatially uniform value of `component` across the whole domain." +function set_flat_profile!(d::Domain1D, component::AbstractString, value) + check(LibCantera.domain_setFlatProfile(d.handle, component, Float64(value))) + return d +end + +"Set a uniform grid of `points` nodes spanning `[start, start+length]` [m]." +function setup_uniform_grid!(d::Domain1D, points::Integer, length::Real, + start::Real=0.0) + check(LibCantera.domain_setupUniformGrid(d.handle, Int32(points), Float64(length), + Float64(start))) + return d +end + +"Set an explicit grid from a vector of positions [m]." +function setup_grid!(d::Domain1D, z) + zz = as_f64(collect(z)) + check(LibCantera.domain_setupGrid(d.handle, Int32(length(zz)), pointer(zz))) + return d +end + +function Base.show(io::IO, d::Domain1D) + if d.closed + print(io, "Domain1D()") + else + print(io, "Domain1D(\"", domain_type(d), "\", ", n_points(d), " points)") + end +end + +# --------------------------------------------------------------------------- +# FreeFlame: a freely-propagating premixed flame. +# --------------------------------------------------------------------------- + +""" + FreeFlame(gas; width=0.03, points=8) + +Assemble a freely-propagating premixed flame from the current state of `gas`, +mirroring Python's `cantera.FreeFlame`. The unburned mixture state (temperature, +pressure, composition) is taken from `gas` at construction time. + +The flame consists of three domains in solver order: an inlet boundary, a +`free-flow` region of the given `width` [m] discretized with `points` initial +grid nodes, and an outlet boundary. Call [`solve!`](@ref) to compute the +solution; the laminar flame speed is then available from [`flame_speed`](@ref). +""" +mutable struct FreeFlame <: CanteraObject + gas::Solution + inlet::Domain1D + flow::Domain1D + outlet::Domain1D + sim::Int32 + rho_u::Float64 # unburned density [kg/m^3], for the flame-speed eigenvalue + Tu::Float64 # unburned temperature [K] + P::Float64 # pressure [Pa] + Yu::Vector{Float64} # unburned mass fractions + closed::Bool +end + +# Initial (non-uniform) flame grid used by Python's FreeFlame for a given width. +_initial_flame_grid(width) = as_f64([0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0] .* width) + +function FreeFlame(gas::Solution; width::Real=0.03) + # Capture the unburned state from the incoming gas. + Tu = temperature(gas) + P = pressure(gas) + Xu = mole_fractions(gas) + Yu = mass_fractions(gas) + rho_u = density(gas) + + # Create the three domains (order matters: inlet, flow, outlet). + hin = check(LibCantera.domain_newBoundary1D("inlet", gas.handle, "reactants")) + hflw = check(LibCantera.domain_newFlow1D("free-flow", gas.handle, "flame")) + hout = check(LibCantera.domain_newBoundary1D("outlet", gas.handle, "products")) + + inlet = Domain1D(hin, gas, false) + flow = Domain1D(hflw, gas, false) + outlet = Domain1D(hout, gas, false) + + # Flow domain: operating pressure and the initial non-uniform flame grid. + check(LibCantera.flow_setPressure(hflw, Float64(P))) + z0 = _initial_flame_grid(width) + GC.@preserve z0 check(LibCantera.domain_setupGrid(hflw, Int32(length(z0)), + pointer(z0))) + + # Inlet boundary: unburned temperature and composition. The mass flux is + # left at zero here; `set_initial_guess!` seeds it (mdot = 1*rho_u) exactly + # as Python does. + check(LibCantera.bdry_setTemperature(hin, Float64(Tu))) + xu = as_f64(Xu) + GC.@preserve xu check(LibCantera.bdry_setMoleFractions(hin, Int32(length(xu)), + pointer(xu))) + + # Build the simulation from the three domain handles (order matters). + # GC.@preserve keeps `doms` alive across the ccall.s + doms = Int32[hin, hflw, hout] + hsim = GC.@preserve doms check(LibCantera.sim1D_newSim1D(Int32(length(doms)), + pointer(doms))) + + flame = FreeFlame(gas, inlet, flow, outlet, hsim, rho_u, Tu, P, Yu, false) + finalizer(close!, flame) + + _set_initial_guess!(flame) + return flame +end + +# Replicates Python `FreeFlame.set_initial_guess` (locs = [0, 0.3, 0.5, 1.0]). +function _set_initial_guess!(flame::FreeFlame) + hflw = flame.flow.handle + gas = flame.gas + locs = as_f64([0.0, 0.3, 0.5, 1.0]) + + _setprof(comp, vals) = GC.@preserve locs vals check( + LibCantera.domain_setProfile(hflw, comp, Int32(length(locs)), pointer(locs), + Int32(length(vals)), pointer(vals))) + + # Unburned state at the inlet. + set_TPY!(gas, flame.Tu, flame.P, flame.Yu) + + mdot = checkd(LibCantera.bdry_mdot(flame.inlet.handle)) + if mdot == 0.0 + # A nonzero initial guess increases the likelihood of convergence. + mdot = 1.0 * density(gas) + check(LibCantera.bdry_setMdot(flame.inlet.handle, mdot)) + end + + Y0 = flame.Yu + u0 = mdot / density(gas) + T0 = flame.Tu + + # Burned (adiabatic-flame) state by equilibrating at constant H, P. + equilibrate!(gas, "HP") + Teq = temperature(gas) + Yeq = mass_fractions(gas) + u1 = mdot / density(gas) + + _setprof("velocity", as_f64([u0, u0, u1, u1])) + _setprof("T", as_f64([T0, T0, Teq, Teq])) + + # Pick the fixed-temperature point, reusing an existing grid point when a + # reasonable choice exists (Python's exact criterion). + T = solution_profile(flame.flow, "T") + Tmid = 0.75 * T0 + 0.25 * Teq + i = findlast(<(Tmid), T) # 1-based analogue of numpy flatnonzero[-1] + Tfixed = if Tmid - T[i] < 0.2 * (Tmid - T0) + T[i] + elseif T[i + 1] - Tmid < 0.2 * (Teq - Tmid) + T[i + 1] + else + Tmid + end + check(LibCantera.sim1D_setFixedTemperature(flame.sim, Float64(Tfixed))) + + # Species mass-fraction profiles: unburned -> burned. + for n in 1:n_species(gas) + _setprof(species_name(gas, n), as_f64([Y0[n], Y0[n], Yeq[n], Yeq[n]])) + end + return flame +end + +""" + set_refine_criteria!(flame; ratio=10.0, slope=0.8, curve=0.8, prune=0.0) + +Set the grid-refinement criteria on the flow domain (domain index 1). Defaults +match Python's `FlameBase.set_refine_criteria`. +""" +function set_refine_criteria!(flame::FreeFlame; ratio=10.0, slope=0.8, curve=0.8, + prune=0.0) + check(LibCantera.sim1D_setRefineCriteria(flame.sim, Int32(1), Float64(ratio), + Float64(slope), Float64(curve), + Float64(prune))) + return flame +end + +""" + set_inlet!(flame; T=nothing, X=nothing, mdot=nothing) + +Update the inlet (reactants) boundary temperature, composition and/or mass flux. +""" +function set_inlet!(flame::FreeFlame; T=nothing, X=nothing, mdot=nothing) + h = flame.inlet.handle + T === nothing || check(LibCantera.bdry_setTemperature(h, Float64(T))) + mdot === nothing || check(LibCantera.bdry_setMdot(h, Float64(mdot))) + if X !== nothing + if X isa AbstractString + check(LibCantera.bdry_setMoleFractionsByName(h, X)) + else + x = as_f64(collect(X)) + check(LibCantera.bdry_setMoleFractions(h, Int32(length(x)), pointer(x))) + end + end + return flame +end + +"Fix the flame temperature at `T` [K] (used by the free-flame eigenvalue solve)." +function set_fixed_temperature!(flame::FreeFlame, T::Real) + check(LibCantera.sim1D_setFixedTemperature(flame.sim, Float64(T))) + return flame +end + +# Regrid the flow domain to a uniform `N`-point grid over `[zmin, zmax]`. +function _regrid!(flame::FreeFlame, N::Integer, zmin::Real, zmax::Real) + z = collect(range(zmin, zmax; length=N)) + GC.@preserve z check(LibCantera.domain_setupGrid(flame.flow.handle, + Int32(length(z)), pointer(z))) + return flame +end + +# One `sim1D_solve`; returns `true` on success, `false` on a CanteraError. +function _try_solve(sim, ll, refine) + try + check(LibCantera.sim1D_solve(sim, ll, Int32(refine))) + return true + catch + return false + end +end + +# Signals that the flow domain is too narrow to contain the flame; caught by +# `solve!`, which doubles the domain width (mirrors Python's `DomainTooNarrow`). +struct _DomainTooNarrow <: Exception end + +# Python's `check_width` criterion: the domain is too narrow if the temperature +# gradient at either edge is significant compared to the average gradient. +function _width_ok(flame::FreeFlame) + T = solution_profile(flame.flow, "T") + x = grid(flame.flow) + length(x) < 3 && return true + mRef = (T[end] - T[1]) / (x[end] - x[1]) + mRef == 0 && return true + mLeft = (T[2] - T[1]) / (x[2] - x[1]) / mRef + mRight = (T[end-2] - T[end]) / (x[end-2] - x[end]) / mRef + return !(mLeft > 0.02 || mRight > 0.02) +end + +# Replicates the Cython `Sim1D.solve(auto=True)` staged loop: progressively +# finer uniform grids (8, 12, 24, 48 points), re-seeding the initial guess on +# each, trying energy-on first and falling back to a fixed-temperature solve, +# then a refinement pass. Returns `true` once a solution is found. +function _staged_solve!(flame::FreeFlame, ll::Int32, refine_grid::Bool) + sim = flame.sim + flow = flame.flow.handle + g = grid(flame.flow) + zmin, zmax = g[1], g[end] + + nPoints = [n_points(flame.flow), 12, 24, 48] + solved = false + for N in nPoints + if N != n_points(flame.flow) + _regrid!(flame, N, zmin, zmax) + end + _set_initial_guess!(flame) + + # Try with the energy equation enabled first (usually works). + check(LibCantera.flow_setEnergyEnabled(flow, Int32(1))) + solved = _try_solve(sim, ll, false) + + # Fall back to a fixed-temperature solve, then re-enable energy. + if !solved + check(LibCantera.flow_setEnergyEnabled(flow, Int32(0))) + solved = _try_solve(sim, ll, false) + if solved + check(LibCantera.flow_setEnergyEnabled(flow, Int32(1))) + solved = _try_solve(sim, ll, false) + end + end + + # Emulate Python's steady-state width callback: it fires after this + # (pre-refinement) steady solve on the coarse grid, so check here. + solved && !_width_ok(flame) && throw(_DomainTooNarrow()) + + # Refinement pass; a non-extinct converged solution ends the loop. + # (`extinct()` is always false for a free flame.) + if solved && refine_grid + solved = _try_solve(sim, ll, true) + solved && !_width_ok(flame) && throw(_DomainTooNarrow()) + solved && break + end + + refine_grid || break + end + solved || throw(CanteraError("Could not find a solution for the 1D problem")) + return solved +end + +""" + solve!(flame; loglevel=0, refine_grid=true, auto=true) + +Solve the flame. + +With `auto=true` (the default) this reproduces Python's +`FreeFlame.solve(auto=True)` exactly: an internal staged multi-grid schedule +wrapped in a domain-widening loop. After each staged solve the temperature +gradients at the domain edges are checked; if the flame is too close to a +boundary the grid is doubled (and refined) and the staged solve is repeated, +up to 12 times. + +!!! note + Python additionally installs the width check as a *steady-state callback* + inside the C++ solver, so it can abort and widen mid-solve. The CLib does + not expose steady callbacks, so here the identical width criterion is applied + *after* each completed staged solve instead. This is the one piece of the + reference algorithm that is emulated externally rather than in-solver; the + final converged result is unaffected for the usual case. + +With `auto=false` a single `sim1D_solve(loglevel, refine_grid)` is issued using +the current refine criteria (advanced use). +""" +function solve!(flame::FreeFlame; loglevel::Integer=0, refine_grid::Bool=true, + auto::Bool=true) + flame.closed && throw(CanteraError("FreeFlame is closed")) + ll = Int32(loglevel) + sim = flame.sim + flow = flame.flow.handle + + if !auto + if !_try_solve(sim, ll, refine_grid) + throw(CanteraError("Could not find a solution for the 1D problem")) + end + return flame + end + + # Domain-widening loop (mirrors FreeFlame.solve check_width / grid *= 2). + for _ in 1:12 + try + _staged_solve!(flame, ll, refine_grid) + break + catch err + err isa _DomainTooNarrow || rethrow() + # Too narrow: double the domain, then re-run the *staged* solve from + # scratch. Mirrors Python's `self.flame.grid *= 2; self.refine(...)`: + # `sim1D_refine` only adapts the grid (inserting points per the + # refine criteria), it does not solve. + znew = grid(flame.flow) .* 2 + GC.@preserve znew check(LibCantera.domain_setupGrid(flow, + Int32(length(znew)), pointer(znew))) + refine_grid && check(LibCantera.sim1D_refine(sim, ll)) + end + end + return flame +end + +# ---- post-solve accessors --------------------------------------------------- + +"Grid points of the flow (flame) domain [m]." +grid(flame::FreeFlame) = grid(flame.flow) + +"Per-point profile of `component` in the flow domain." +solution_profile(flame::FreeFlame, component::AbstractString) = + solution_profile(flame.flow, component) + +"Temperature profile across the flame [K]." +flame_T(flame::FreeFlame) = solution_profile(flame.flow, "T") + +"Axial velocity profile across the flame [m/s]." +flame_velocity(flame::FreeFlame) = solution_profile(flame.flow, "velocity") + +"Mole-fraction profile of `species` across the flame." +flame_X(flame::FreeFlame, species::AbstractString) = + solution_profile(flame.flow, species) + +""" + flame_speed(flame) -> Float64 + +Laminar flame speed [m/s]: the inlet (unburned) axial velocity. For a freely +propagating flame the inlet mass flux is an eigenvalue updated by the solver, so +this reads the current inlet `mdot` and divides by the unburned density. +""" +flame_speed(flame::FreeFlame) = + checkd(LibCantera.bdry_mdot(flame.inlet.handle)) / flame.rho_u + +"Number of grid points in the flow domain." +n_points(flame::FreeFlame) = n_points(flame.flow) + +function close!(flame::FreeFlame) + flame.closed && return flame + flame.closed = true + try + LibCantera.sim1D_del(flame.sim) + catch + end + for d in (flame.inlet, flame.flow, flame.outlet) + d.closed = true + try + LibCantera.domain_del(d.handle) + catch + end + end + return flame +end + +function Base.show(io::IO, flame::FreeFlame) + if flame.closed + print(io, "FreeFlame()") + else + print(io, "FreeFlame(", n_points(flame), " grid points)") + end +end + +# --------------------------------------------------------------------------- +# BurnerFlame: a burner-stabilized flat flame. +# --------------------------------------------------------------------------- + +""" + BurnerFlame(gas; width=0.03, mdot=nothing) + +Assemble a burner-stabilized flat flame from the current state of `gas`, +mirroring Python's `cantera.BurnerFlame`. The unburned mixture state +(temperature, pressure, composition) is taken from `gas` at construction time. + +Unlike a [`FreeFlame`](@ref), the burner inlet has a user-prescribed mass flux +`mdot` [kg/m^2/s]; there is no flame-speed eigenvalue and no fixed-temperature +anchor. If `mdot` is not given it defaults to `0.4 * ρ_unburned`. + +The flame consists of three domains in solver order: a burner (inlet) boundary, +an `unstrained-flow` region of the given `width` [m], and an outlet boundary. +Call [`solve!`](@ref) to compute the solution. Use [`set_burner!`](@ref) to +change the burner temperature, composition, or mass flux. +""" +mutable struct BurnerFlame <: CanteraObject + gas::Solution + burner::Domain1D + flow::Domain1D + outlet::Domain1D + sim::Int32 + rho_u::Float64 # unburned density [kg/m^3] + Tu::Float64 # unburned temperature [K] + P::Float64 # pressure [Pa] + Yu::Vector{Float64} # unburned mass fractions + closed::Bool +end + +# Initial (non-uniform) grid used by Python's BurnerFlame for a given width. +_initial_burner_grid(width) = as_f64([0.0, 0.1, 0.2, 0.3, 0.5, 0.7, 1.0] .* width) + +function BurnerFlame(gas::Solution; width::Real=0.03, mdot=nothing) + # Capture the unburned (burner) state from the incoming gas. + Tu = temperature(gas) + P = pressure(gas) + Xu = mole_fractions(gas) + Yu = mass_fractions(gas) + rho_u = density(gas) + + # Create the three domains (order matters: burner, flow, outlet). + hbrn = check(LibCantera.domain_newBoundary1D("inlet", gas.handle, "burner")) + hflw = check(LibCantera.domain_newFlow1D("unstrained-flow", gas.handle, "flame")) + hout = check(LibCantera.domain_newBoundary1D("outlet", gas.handle, "products")) + + burner = Domain1D(hbrn, gas, false) + flow = Domain1D(hflw, gas, false) + outlet = Domain1D(hout, gas, false) + + # Flow domain: operating pressure and the initial non-uniform grid. + check(LibCantera.flow_setPressure(hflw, Float64(P))) + z0 = _initial_burner_grid(width) + GC.@preserve z0 check(LibCantera.domain_setupGrid(hflw, Int32(length(z0)), + pointer(z0))) + + # Burner boundary: unburned temperature, composition, and prescribed mdot. + check(LibCantera.bdry_setTemperature(hbrn, Float64(Tu))) + xu = as_f64(Xu) + GC.@preserve xu check(LibCantera.bdry_setMoleFractions(hbrn, Int32(length(xu)), + pointer(xu))) + md = mdot === nothing ? 0.4 * rho_u : Float64(mdot) + check(LibCantera.bdry_setMdot(hbrn, md)) + + doms = Int32[hbrn, hflw, hout] + hsim = GC.@preserve doms check(LibCantera.sim1D_newSim1D(Int32(length(doms)), + pointer(doms))) + + flame = BurnerFlame(gas, burner, flow, outlet, hsim, rho_u, Tu, P, Yu, false) + finalizer(close!, flame) + + _set_initial_guess!(flame) + return flame +end + +# Replicates Python `BurnerFlame.set_initial_guess` (locs = [0, 0.2, 1.0]). +function _set_initial_guess!(flame::BurnerFlame) + hflw = flame.flow.handle + gas = flame.gas + locs = as_f64([0.0, 0.2, 1.0]) + + _setprof(comp, vals) = GC.@preserve locs vals check( + LibCantera.domain_setProfile(hflw, comp, Int32(length(locs)), pointer(locs), + Int32(length(vals)), pointer(vals))) + + # Unburned state at the burner. + set_TPY!(gas, flame.Tu, flame.P, flame.Yu) + mdot = checkd(LibCantera.bdry_mdot(flame.burner.handle)) + + Y0 = flame.Yu + u0 = mdot / density(gas) + T0 = flame.Tu + + # Burned (adiabatic-flame) state by equilibrating at constant H, P. + equilibrate!(gas, "HP") + Teq = temperature(gas) + Yeq = mass_fractions(gas) + u1 = mdot / density(gas) + + _setprof("velocity", as_f64([u0, u1, u1])) + _setprof("T", as_f64([T0, Teq, Teq])) + for n in 1:n_species(gas) + _setprof(species_name(gas, n), as_f64([Y0[n], Yeq[n], Yeq[n]])) + end + return flame +end + +""" + set_burner!(flame; T=nothing, X=nothing, mdot=nothing) + +Update the burner (reactants) boundary temperature, composition and/or mass flux +[kg/m^2/s]. +""" +function set_burner!(flame::BurnerFlame; T=nothing, X=nothing, mdot=nothing) + h = flame.burner.handle + if T !== nothing + check(LibCantera.bdry_setTemperature(h, Float64(T))) + flame.Tu = Float64(T) + end + mdot === nothing || check(LibCantera.bdry_setMdot(h, Float64(mdot))) + if X !== nothing + if X isa AbstractString + check(LibCantera.bdry_setMoleFractionsByName(h, X)) + else + x = as_f64(collect(X)) + check(LibCantera.bdry_setMoleFractions(h, Int32(length(x)), pointer(x))) + end + # Refresh the cached unburned composition/density from the new state. + set_TPX!(flame.gas, flame.Tu, flame.P, X isa AbstractString ? X : collect(X)) + flame.Yu = mass_fractions(flame.gas) + flame.rho_u = density(flame.gas) + end + return flame +end + +""" + set_refine_criteria!(flame::BurnerFlame; ratio=10.0, slope=0.8, curve=0.8, + prune=0.0) + +Set the grid-refinement criteria on the flow domain (domain index 1). +""" +function set_refine_criteria!(flame::BurnerFlame; ratio=10.0, slope=0.8, curve=0.8, + prune=0.0) + check(LibCantera.sim1D_setRefineCriteria(flame.sim, Int32(1), Float64(ratio), + Float64(slope), Float64(curve), + Float64(prune))) + return flame +end + +_regrid!(flame::BurnerFlame, N::Integer, zmin::Real, zmax::Real) = + (z = collect(range(zmin, zmax; length=N)); + GC.@preserve z check(LibCantera.domain_setupGrid(flame.flow.handle, + Int32(length(z)), pointer(z))); + flame) + +# Staged multi-grid solve for a burner flame: no fixed-temperature eigenvalue. +function _staged_solve!(flame::BurnerFlame, ll::Int32, refine_grid::Bool) + sim = flame.sim + flow = flame.flow.handle + g = grid(flame.flow) + zmin, zmax = g[1], g[end] + + nPoints = [n_points(flame.flow), 12, 24, 48] + solved = false + for N in nPoints + if N != n_points(flame.flow) + _regrid!(flame, N, zmin, zmax) + end + _set_initial_guess!(flame) + + check(LibCantera.flow_setEnergyEnabled(flow, Int32(1))) + solved = _try_solve(sim, ll, false) + if !solved + check(LibCantera.flow_setEnergyEnabled(flow, Int32(0))) + solved = _try_solve(sim, ll, false) + if solved + check(LibCantera.flow_setEnergyEnabled(flow, Int32(1))) + solved = _try_solve(sim, ll, false) + end + end + if solved && refine_grid + solved = _try_solve(sim, ll, true) + solved && break + end + refine_grid || break + end + solved || throw(CanteraError("Could not find a solution for the 1D problem")) + return solved +end + +""" + solve!(flame::BurnerFlame; loglevel=0, refine_grid=true, auto=true) + +Solve the burner-stabilized flame. With `auto=true` a staged multi-grid +schedule is used; with `auto=false` a single `sim1D_solve` is issued using the +current refine criteria. Unlike a free flame there is no flame-speed eigenvalue +or fixed-temperature anchor; the burner mass flux is a fixed input. +""" +function solve!(flame::BurnerFlame; loglevel::Integer=0, refine_grid::Bool=true, + auto::Bool=true) + flame.closed && throw(CanteraError("BurnerFlame is closed")) + ll = Int32(loglevel) + if !auto + if !_try_solve(flame.sim, ll, refine_grid) + throw(CanteraError("Could not find a solution for the 1D problem")) + end + return flame + end + _staged_solve!(flame, ll, refine_grid) + return flame +end + +# ---- post-solve accessors --------------------------------------------------- + +grid(flame::BurnerFlame) = grid(flame.flow) +solution_profile(flame::BurnerFlame, component::AbstractString) = + solution_profile(flame.flow, component) +flame_T(flame::BurnerFlame) = solution_profile(flame.flow, "T") +flame_velocity(flame::BurnerFlame) = solution_profile(flame.flow, "velocity") +flame_X(flame::BurnerFlame, species::AbstractString) = + solution_profile(flame.flow, species) +n_points(flame::BurnerFlame) = n_points(flame.flow) + +"Prescribed burner mass flux [kg/m^2/s]." +burner_mdot(flame::BurnerFlame) = checkd(LibCantera.bdry_mdot(flame.burner.handle)) + +function close!(flame::BurnerFlame) + flame.closed && return flame + flame.closed = true + try + LibCantera.sim1D_del(flame.sim) + catch + end + for d in (flame.burner, flame.flow, flame.outlet) + d.closed = true + try + LibCantera.domain_del(d.handle) + catch + end + end + return flame +end + +function Base.show(io::IO, flame::BurnerFlame) + if flame.closed + print(io, "BurnerFlame()") + else + print(io, "BurnerFlame(", n_points(flame), " grid points)") + end +end + +export Domain1D, FreeFlame, BurnerFlame, set_burner!, burner_mdot, + n_points, n_components, component_name, + component_names, domain_type, grid, value, solution_profile, set_profile!, + set_flat_profile!, setup_uniform_grid!, setup_grid!, set_refine_criteria!, + solve!, flame_speed, flame_T, flame_X, flame_velocity, + set_fixed_temperature!, set_inlet! diff --git a/interfaces/julia/src/rdiag.jl b/interfaces/julia/src/rdiag.jl new file mode 100644 index 0000000000..178df02827 --- /dev/null +++ b/interfaces/julia/src/rdiag.jl @@ -0,0 +1,189 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# ReactionPathDiagram: a parity wrapper over the CLib `rdiag_*` reaction-path +# diagram API (ctrdiag.h). A diagram holds a handle into the CLib rdiag cabinet +# and a reference to the Solution it was built from (to keep the underlying +# Kinetics manager alive for the diagram's lifetime). + +""" + ReactionPathDiagram(gas::Solution, element::AbstractString) + +Reaction-path diagram tracing the flow of a chemical `element` (e.g. `"C"`, +`"H"`, `"O"`, `"N"`) through the reaction network of `gas` at its current state. + +Set options (see [`set_threshold!`](@ref), [`set_flow_type!`](@ref), ...), then +call [`build!`](@ref) and retrieve the result with [`get_dot`](@ref) (Graphviz), +[`get_data`](@ref) or [`get_log`](@ref). + +```julia +gas = Solution("gri30.yaml") +set_TPX!(gas, 2000.0, one_atm, "CH4:1, O2:2, N2:7.52") +d = ReactionPathDiagram(gas, "C") +set_threshold!(d, 0.01) +set_flow_type!(d, "NetFlow") +build!(d) +print(get_dot(d)) +``` +""" +mutable struct ReactionPathDiagram <: CanteraObject + handle::Int32 + solution::Solution # keep the Kinetics manager alive for the diagram + closed::Bool +end + +function ReactionPathDiagram(gas::Solution, element::AbstractString) + kin = _kinetics_handle(gas) + h = check(LibCantera.rdiag_newReactionPathDiagram(kin, element)) + d = ReactionPathDiagram(h, gas, false) + finalizer(close!, d) + return d +end + +""" + close!(d::ReactionPathDiagram) + +Release the Cantera-side object backing `d`. Safe to call multiple times; the +finalizer calls this automatically. +""" +function close!(d::ReactionPathDiagram) + d.closed && return d + d.closed = true + try + LibCantera.rdiag_del(d.handle) + catch + end + return d +end + +function Base.show(io::IO, d::ReactionPathDiagram) + if d.closed + print(io, "ReactionPathDiagram()") + else + print(io, "ReactionPathDiagram(handle=", d.handle, ")") + end +end + +# ---- build / output --------------------------------------------------------- + +""" + build!(d::ReactionPathDiagram) + +Build (or rebuild) the diagram from the current state of the owning Solution. +Call after setting options and before retrieving output. +""" +function build!(d::ReactionPathDiagram) + check(LibCantera.rdiag_build(d.handle)) + return d +end + +"Graphviz `dot` description of the diagram (call [`build!`](@ref) first)." +get_dot(d::ReactionPathDiagram) = + get_string((n, b) -> LibCantera.rdiag_getDot(d.handle, n, b)) + +"Raw flow data underlying the diagram (call [`build!`](@ref) first)." +get_data(d::ReactionPathDiagram) = + get_string((n, b) -> LibCantera.rdiag_getData(d.handle, n, b)) + +"Diagnostic log accumulated while building the diagram." +get_log(d::ReactionPathDiagram) = + get_string((n, b) -> LibCantera.rdiag_getLog(d.handle, n, b)) + +# ---- scalar (double) properties --------------------------------------------- + +for (getter, setter, cget, cset, doc) in ( + (:threshold, :set_threshold!, :rdiag_threshold, :rdiag_setThreshold, + "Minimum relative flux value that will be plotted."), + (:bold_threshold, :set_bold_threshold!, :rdiag_boldThreshold, + :rdiag_setBoldThreshold, "Minimum relative flux for bold lines."), + (:normal_threshold, :set_normal_threshold!, :rdiag_normalThreshold, + :rdiag_setNormalThreshold, "Maximum relative flux for dashed lines."), + (:label_threshold, :set_label_threshold!, :rdiag_labelThreshold, + :rdiag_setLabelThreshold, "Minimum relative flux for labels."), + (:scale, :set_scale!, :rdiag_scale, :rdiag_setScale, + "Scaling factor for the fluxes; -1 normalizes by the maximum net flux."), + (:arrow_width, :set_arrow_width!, :rdiag_arrowWidth, :rdiag_setArrowWidth, + "Arrow width; if negative, arrows scale with the flux value instead."), + ) + setter_doc = "Set [`$getter`](@ref)." + @eval begin + $getter(d::ReactionPathDiagram) = checkd(LibCantera.$cget(d.handle)) + @doc $doc $getter + function $setter(d::ReactionPathDiagram, v::Real) + check(LibCantera.$cset(d.handle, Float64(v))) + return d + end + @doc $setter_doc $setter + end +end + +# ---- string properties ------------------------------------------------------ + +for (getter, setter, cget, cset, doc) in ( + (:flow_type, :set_flow_type!, :rdiag_flowType, :rdiag_setFlowType, + "Way flows are drawn: `\"NetFlow\"` or `\"OneWayFlow\"`."), + (:title, :set_title!, :rdiag_title, :rdiag_setTitle, + "Reaction path diagram title."), + (:font, :set_font!, :rdiag_font, :rdiag_setFont, + "Font used for the diagram."), + (:bold_color, :set_bold_color!, :rdiag_boldColor, :rdiag_setBoldColor, + "Color for bold lines."), + (:normal_color, :set_normal_color!, :rdiag_normalColor, + :rdiag_setNormalColor, "Color for normal-weight lines."), + (:dashed_color, :set_dashed_color!, :rdiag_dashedColor, + :rdiag_setDashedColor, "Color for dashed lines."), + (:dot_options, :set_dot_options!, :rdiag_dotOptions, :rdiag_setDotOptions, + "Options passed through to the `dot` program."), + ) + setter_doc = "Set [`$getter`](@ref)." + @eval begin + $getter(d::ReactionPathDiagram) = + get_string((n, b) -> LibCantera.$cget(d.handle, n, b)) + @doc $doc $getter + function $setter(d::ReactionPathDiagram, v::AbstractString) + check(LibCantera.$cset(d.handle, v)) + return d + end + @doc $setter_doc $setter + end +end + +# ---- boolean property ------------------------------------------------------- + +"Whether reaction details are shown as edge labels." +show_details(d::ReactionPathDiagram) = + check(LibCantera.rdiag_showDetails(d.handle)) != 0 + +"Enable or disable showing reaction details as edge labels." +function set_show_details!(d::ReactionPathDiagram, v::Bool) + check(LibCantera.rdiag_setShowDetails(d.handle, Int32(v))) + return d +end + +# ---- other operations ------------------------------------------------------- + +""" + display_only!(d::ReactionPathDiagram, species_index) + +Restrict the diagram to flows into/out of a single species (1-based index), or +pass `:all` (or a negative index) to show all species again. +""" +function display_only!(d::ReactionPathDiagram, k::Integer) + # CLib expects a 0-based index; -1 means "all". + idx = k < 0 ? Int32(-1) : Int32(k - 1) + check(LibCantera.rdiag_displayOnly(d.handle, idx)) + return d +end +function display_only!(d::ReactionPathDiagram, s::Symbol) + s === :all || throw(ArgumentError("expected a species index or :all, got :$s")) + return display_only!(d, -1) +end + +export ReactionPathDiagram, build!, get_dot, get_data, get_log, + threshold, set_threshold!, bold_threshold, set_bold_threshold!, + normal_threshold, set_normal_threshold!, label_threshold, + set_label_threshold!, scale, set_scale!, arrow_width, set_arrow_width!, + flow_type, set_flow_type!, title, set_title!, font, set_font!, + show_details, set_show_details!, bold_color, set_bold_color!, + normal_color, set_normal_color!, dashed_color, set_dashed_color!, + dot_options, set_dot_options!, display_only! diff --git a/interfaces/julia/src/reaction.jl b/interfaces/julia/src/reaction.jl new file mode 100644 index 0000000000..e4d94694a7 --- /dev/null +++ b/interfaces/julia/src/reaction.jl @@ -0,0 +1,56 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Reaction object wrapper. Reactions are usually accessed through a Kinetics +# manager (see `reaction_equations`), but a standalone wrapper is provided for +# completeness and parity with the CLib `ctrxn` library. + +""" + Reaction + +Wrapper around a Cantera `Reaction` CLib handle. +""" +mutable struct Reaction <: CanteraObject + handle::Int32 + closed::Bool + + function Reaction(handle::Integer) + r = new(Int32(handle), false) + finalizer(close!, r) + return r + end +end + +""" + reaction(gas, i) -> Reaction + +Return reaction `i` (1-based) from a Kinetics manager as a [`Reaction`](@ref). +""" +function reaction(g::KineticsLike, i::Integer) + h = check(LibCantera.kin_reaction(_kinetics_handle(g), Int32(i - 1))) + return Reaction(h) +end + +function close!(r::Reaction) + r.closed && return r + r.closed = true + try + LibCantera.rxn_del(r.handle) + catch + end + return r +end + +"Reaction equation string." +equation(r::Reaction) = get_string((n, b) -> LibCantera.rxn_equation(r.handle, n, b)) + +"Reaction type string (e.g. `\"reaction\"`, `\"three-body\"`)." +reaction_type(r::Reaction) = get_string((n, b) -> LibCantera.rxn_type(r.handle, n, b)) + +"Whether the reaction involves a third body." +uses_third_body(r::Reaction) = check(LibCantera.rxn_usesThirdBody(r.handle)) != 0 + +Base.show(io::IO, r::Reaction) = + print(io, "Reaction(\"", r.closed ? "" : equation(r), "\")") + +export Reaction, reaction, equation, reaction_type, uses_third_body diff --git a/interfaces/julia/src/reactor.jl b/interfaces/julia/src/reactor.jl new file mode 100644 index 0000000000..408c0fbc63 --- /dev/null +++ b/interfaces/julia/src/reactor.jl @@ -0,0 +1,199 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Reactor wrappers. A Reactor holds a handle into the CLib reactor cabinet and +# a reference to the Solution it was built from (to keep it alive). + +""" + Reactor + +Wrapper around a Cantera reactor CLib handle. Construct one of the concrete +helpers ([`IdealGasReactor`](@ref), [`Reactor`](@ref), [`ConstPressureReactor`](@ref), +[`IdealGasConstPressureReactor`](@ref)) rather than this type directly. +""" +mutable struct Reactor <: CanteraObject + handle::Int32 + solution::Solution # keep the phase alive for the reactor's lifetime + closed::Bool +end + +function _new_reactor(model::AbstractString, gas::Solution, name::AbstractString) + h = check(LibCantera.reactor_new(model, gas.handle, Int32(1), name)) + r = Reactor(h, gas, false) + finalizer(close!, r) + return r +end + +""" + IdealGasReactor(gas; name="") + +Zero-dimensional constant-volume reactor with an ideal-gas energy equation. +""" +IdealGasReactor(gas::Solution; name::AbstractString="") = + _new_reactor("IdealGasReactor", gas, name) + +""" + Reactor(gas; name="") + +General constant-volume reactor. +""" +Reactor(gas::Solution; name::AbstractString="") = + _new_reactor("Reactor", gas, name) + +""" + ConstPressureReactor(gas; name="") + +Constant-pressure reactor. +""" +ConstPressureReactor(gas::Solution; name::AbstractString="") = + _new_reactor("ConstPressureReactor", gas, name) + +""" + IdealGasConstPressureReactor(gas; name="") + +Constant-pressure reactor with an ideal-gas energy equation. +""" +IdealGasConstPressureReactor(gas::Solution; name::AbstractString="") = + _new_reactor("IdealGasConstPressureReactor", gas, name) + +function close!(r::Reactor) + r.closed && return r + r.closed = true + try + LibCantera.reactor_del(r.handle) + catch + end + return r +end + +# ---- scalar reactor state --------------------------------------------------- + +for (jl, c) in ( + (:mass, :reactor_mass), + (:volume, :reactor_volume), + (:density, :reactor_density), + (:temperature, :reactor_temperature), + (:pressure, :reactor_pressure), + (:enthalpy_mass, :reactor_enthalpy_mass), + (:internal_energy_mass,:reactor_intEnergy_mass), + ) + @eval $jl(r::Reactor) = checkd(LibCantera.$c(r.handle)) +end + +"Reactor name." +name(r::Reactor) = get_string((n, b) -> LibCantera.reactor_name(r.handle, n, b)) + +"Reactor type string." +reactor_type(r::Reactor) = get_string((n, b) -> LibCantera.reactor_type(r.handle, n, b)) + +""" + reactor_phase(r::Reactor) -> Solution + +The `Solution` actually integrated by the reactor. This is the reactor's own +clone of the phase it was constructed from (see [`IdealGasReactor`](@ref) and +friends), and reflects its current state, unlike the original `Solution` +passed to the constructor. +""" +reactor_phase(r::Reactor) = Solution(check(LibCantera.reactor_phase(r.handle))) + +"Mass fractions in the reactor." +function mass_fractions(r::Reactor) + nsp = n_species(r.solution) + return get_array(nsp, (n, b) -> LibCantera.reactor_massFractions(r.handle, n, b)) +end + +"Enable/disable chemistry in the reactor." +function set_energy_enabled!(r::Reactor, flag::Bool) + check(LibCantera.reactor_setEnergyEnabled(r.handle, Int32(flag))) + return r +end + +"Enable/disable the energy equation in the reactor." +function set_chemistry_enabled!(r::Reactor, flag::Bool) + check(LibCantera.reactor_setChemistryEnabled(r.handle, Int32(flag))) + return r +end + +"Set the initial reactor volume [m^3]." +function set_initial_volume!(r::Reactor, vol) + check(LibCantera.reactor_setInitialVolume(r.handle, Float64(vol))) + return r +end + +Base.show(io::IO, r::Reactor) = + print(io, r.closed ? "Reactor()" : "$(reactor_type(r))(\"$(name(r))\")") + +# ---- reservoirs ------------------------------------------------------------- + +""" + Reservoir(gas; name="") + +A reactor with a fixed thermodynamic state, used as a boundary (source or sink) +for [`FlowDevice`](@ref)s and [`Wall`](@ref)s. Its state never changes during +integration. +""" +Reservoir(gas::Solution; name::AbstractString="") = + _new_reactor("Reservoir", gas, name) + +# ---- surfaces --------------------------------------------------------------- + +""" + reactor_area(r) -> Float64 + +Wall/surface area associated with the reactor [m^2]. +""" +area(r::Reactor) = checkd(LibCantera.reactor_area(r.handle)) + +"Set the wall/surface area of the reactor [m^2]." +function set_area!(r::Reactor, a) + check(LibCantera.reactor_setArea(r.handle, Float64(a))) + return r +end + +""" + ReactorSurface(surf, reactor; name="") + +A surface on which heterogeneous reactions take place, coupling the surface +phase `surf` (a [`Solution`](@ref) for an interface phase) to a bulk `reactor`. +""" +function ReactorSurface(surf::Solution, reactor::Reactor; name::AbstractString="") + reactors = Int32[reactor.handle] + h = check(LibCantera.reactor_newSurface(surf.handle, Int32(length(reactors)), + pointer(reactors), Int32(0), name)) + rs = Reactor(h, surf, false) + finalizer(close!, rs) + return rs +end + +# ---- mass flow rate (for flow reactors / surfaces) -------------------------- + +"Mass flow rate through the reactor [kg/s]." +mass_flow_rate(r::Reactor) = checkd(LibCantera.reactor_massFlowRate(r.handle)) + +"Set the mass flow rate through the reactor [kg/s]." +function set_mass_flow_rate!(r::Reactor, mdot) + check(LibCantera.reactor_setMassFlowRate(r.handle, Float64(mdot))) + return r +end + +# ---- sensitivities ---------------------------------------------------------- + +""" + add_sensitivity_reaction!(r, i) + +Mark reaction `i` (1-based) as a sensitivity parameter for reactor `r`. The +reactor must already be part of a [`ReactorNet`](@ref). +""" +function add_sensitivity_reaction!(r::Reactor, i::Integer) + check(LibCantera.reactor_addSensitivityReaction(r.handle, Int32(i - 1))) + return r +end + +"Number of sensitivity parameters associated with the reactor." +n_sens_params(r::Reactor) = Int(check(LibCantera.reactor_nSensParams(r.handle))) + +export Reactor, IdealGasReactor, ConstPressureReactor, IdealGasConstPressureReactor, + set_energy_enabled!, set_chemistry_enabled!, set_initial_volume!, + mass, volume, reactor_phase +export Reservoir, ReactorSurface, mass_flow_rate, set_mass_flow_rate!, + area, set_area!, add_sensitivity_reaction!, n_sens_params diff --git a/interfaces/julia/src/reactornet.jl b/interfaces/julia/src/reactornet.jl new file mode 100644 index 0000000000..67f7a7bccb --- /dev/null +++ b/interfaces/julia/src/reactornet.jl @@ -0,0 +1,156 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# ReactorNet: time integration of a network of reactors. + +""" + ReactorNet(reactors::Vector{Reactor}) + ReactorNet(reactor::Reactor) + +Create a reactor network for time integration. The network keeps references to +its reactors so they are not finalized while integration is in progress. +""" +mutable struct ReactorNet <: CanteraObject + handle::Int32 + reactors::Vector{Reactor} + closed::Bool + + function ReactorNet(reactors::Vector{Reactor}) + isempty(reactors) && + throw(ArgumentError("ReactorNet needs at least one reactor")) + handles = Int32[r.handle for r in reactors] + h = check(LibCantera.reactornet_new(Int32(length(handles)), pointer(handles))) + net = new(h, reactors, false) + finalizer(close!, net) + return net + end +end + +ReactorNet(r::Reactor) = ReactorNet([r]) + +function close!(net::ReactorNet) + net.closed && return net + net.closed = true + try + LibCantera.reactornet_del(net.handle) + catch + end + return net +end + +""" + advance!(net, t) + +Advance the network state to absolute time `t` [s]. +""" +function advance!(net::ReactorNet, t) + check(LibCantera.reactornet_advance(net.handle, Float64(t))) + return net +end + +""" + step!(net) -> Float64 + +Take one internal timestep and return the new time [s]. +""" +step!(net::ReactorNet) = checkd(LibCantera.reactornet_step(net.handle)) + +"Current network time [s]. Extends `Base.time` for `ReactorNet`." +Base.time(net::ReactorNet) = checkd(LibCantera.reactornet_time(net.handle)) + +"Set the initial integration time [s]." +function set_initial_time!(net::ReactorNet, t) + check(LibCantera.reactornet_setInitialTime(net.handle, Float64(t))) + return net +end + +"Set the maximum internal timestep [s]." +function set_max_time_step!(net::ReactorNet, dt) + check(LibCantera.reactornet_setMaxTimeStep(net.handle, Float64(dt))) + return net +end + +""" + set_tolerances!(net; rtol=1e-9, atol=1e-15) + +Set the relative and absolute integrator tolerances. +""" +function set_tolerances!(net::ReactorNet; rtol=1e-9, atol=1e-15) + check(LibCantera.reactornet_setTolerances(net.handle, Float64(rtol), Float64(atol))) + return net +end + +"Relative error tolerance of the network integrator." +rtol(net::ReactorNet) = checkd(LibCantera.reactornet_rtol(net.handle)) + +"Absolute error tolerance of the network integrator." +atol(net::ReactorNet) = checkd(LibCantera.reactornet_atol(net.handle)) + +""" + set_sensitivity_tolerances!(net; rtol=1e-6, atol=1e-6) + +Set the relative and absolute tolerances used for sensitivity analysis. +""" +function set_sensitivity_tolerances!(net::ReactorNet; rtol=1e-6, atol=1e-6) + check(LibCantera.reactornet_setSensitivityTolerances(net.handle, + Float64(rtol), Float64(atol))) + return net +end + +""" + sensitivity(net, component, p, reactor) -> Float64 + +Normalized sensitivity of `component` (e.g. `"temperature"` or a species name) +in `reactor` with respect to sensitivity parameter `p` (1-based). `reactor` may +be a [`Reactor`](@ref) belonging to the network or its 1-based position. +""" +function sensitivity(net::ReactorNet, component::AbstractString, p::Integer, + reactor::Reactor) + idx = findfirst(r -> r === reactor, net.reactors) + idx === nothing && throw(ArgumentError("reactor is not part of this network")) + return sensitivity(net, component, p, idx) +end + +function sensitivity(net::ReactorNet, component::AbstractString, p::Integer, + reactor::Integer) + return checkd(LibCantera.reactornet_sensitivity(net.handle, component, + Int32(p - 1), Int32(reactor - 1))) +end + +""" + n_components(net) -> Int + +Number of state variables (equations) integrated by the network. +""" +n_components(net::ReactorNet) = Int(check(LibCantera.reactornet_neq(net.handle))) + +""" + state(net) -> Vector{Float64} + +Current network state vector, of length [`n_components`](@ref). +""" +state(net::ReactorNet) = + get_array(n_components(net), + (n, b) -> LibCantera.reactornet_getState(net.handle, n, b)) + +"Name of state-vector component `i` (1-based)." +function component_name(net::ReactorNet, i::Integer) + return get_string( + (n, b) -> LibCantera.reactornet_componentName(net.handle, Int32(i - 1), n, b)) +end + +""" + component_names(net) -> Vector{String} + +Names of all state-vector components, aligned with [`state`](@ref). +""" +component_names(net::ReactorNet) = [component_name(net, i) for i in 1:n_components(net)] + +Base.show(io::IO, net::ReactorNet) = + print(io, net.closed ? "ReactorNet()" : + "ReactorNet($(length(net.reactors)) reactor(s), t=$(time(net)) s)") + +export ReactorNet, + advance!, step!, set_initial_time!, set_max_time_step!, set_tolerances!, + rtol, atol, set_sensitivity_tolerances!, sensitivity, state, + n_components, component_name, component_names diff --git a/interfaces/julia/src/solution.jl b/interfaces/julia/src/solution.jl new file mode 100644 index 0000000000..757fcd0c2d --- /dev/null +++ b/interfaces/julia/src/solution.jl @@ -0,0 +1,152 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# The Solution type: the primary user-facing object, mirroring Python's +# `cantera.Solution`. A Solution bundles a ThermoPhase, a Kinetics manager and +# a Transport manager over a single phase, and owns their CLib handles. + +""" + Solution(infile; name="", transport="") + Solution(infile, name; transport="") + +Construct a `Solution` from a Cantera input (YAML) file. + +`name` selects the phase within the file (the first phase is used when empty); +`transport` selects a transport model (the file's default when empty). + +```julia +gas = Solution("gri30.yaml") +gas = Solution("gri30.yaml", "gri30") +``` +""" +mutable struct Solution <: CanteraObject + handle::Int32 + thermo::Int32 + kinetics::Int32 + transport::Int32 + closed::Bool + + function Solution(infile::AbstractString, name::AbstractString=""; + transport::AbstractString="default") + h = check(LibCantera.sol_newSolution(infile, name, transport)) + return Solution(h) + end + + # Internal: wrap an already-registered CLib Solution handle (e.g. one + # obtained from `reactor_phase`) without creating a new one. + function Solution(h::Integer) + th = check(LibCantera.sol_thermo(h)) + # Kinetics/Transport may be absent for some phases; tolerate failures. + kin = _try_handle(() -> LibCantera.sol_kinetics(h)) + tr = _try_handle(() -> LibCantera.sol_transport(h)) + s = new(Int32(h), th, kin, tr, false) + finalizer(close!, s) + return s + end +end + +"Run a handle-returning CLib call, returning -1 instead of throwing on failure." +function _try_handle(f) + h = try + f() + catch + return Int32(-1) + end + return (h == -1 || h == _ERR) ? Int32(-1) : Int32(h) +end + +""" + close!(gas::Solution) + +Release the Cantera-side objects backing `gas`. Safe to call multiple times; +the finalizer calls this automatically, so explicit use is only needed for +deterministic cleanup. +""" +function close!(s::Solution) + s.closed && return s + s.closed = true + # Best-effort deletion; never throw from a finalizer path. + for (h, del) in ((s.transport, LibCantera.trans_del), + (s.kinetics, LibCantera.kin_del), + (s.thermo, LibCantera.thermo_del), + (s.handle, LibCantera.sol_del)) + h == -1 && continue + try + del(h) + catch + end + end + return s +end + +function Base.show(io::IO, s::Solution) + if s.closed + print(io, "Solution()") + else + print(io, "Solution(", n_species(s), " species, ", + n_reactions(s), " reactions)") + end +end + +# ---- handle accessors ------------------------------------------------------- +# Thermo/kinetics/transport functions accept either a Solution or the specific +# sub-phase wrapper. These helpers extract the right handle for dispatch. + +_thermo_handle(s::Solution) = s.thermo +_thermo_handle(t::ThermoPhase) = t.handle + +function _kinetics_handle(s::Solution) + s.kinetics == -1 && + throw(CanteraError("this Solution has no Kinetics manager")) + return s.kinetics +end +_kinetics_handle(k::Kinetics) = k.handle + +function _transport_handle(s::Solution) + s.transport == -1 && + throw(CanteraError("this Solution has no Transport manager")) + return s.transport +end +_transport_handle(t::Transport) = t.handle + +const ThermoLike = Union{Solution,ThermoPhase} +const KineticsLike = Union{Solution,Kinetics} +const TransportLike = Union{Solution,Transport} + +""" + thermo(gas::Solution) -> ThermoPhase + +Return the [`ThermoPhase`](@ref) view of `gas`. +""" +thermo(s::Solution) = ThermoPhase(s.thermo, s) + +""" + kinetics(gas::Solution) -> Kinetics + +Return the [`Kinetics`](@ref) view of `gas`. +""" +kinetics(s::Solution) = Kinetics(_kinetics_handle(s), s) + +""" + transport(gas::Solution) -> Transport + +Return the [`Transport`](@ref) view of `gas`. +""" +transport(s::Solution) = Transport(_transport_handle(s), s) + +""" + name(gas::Solution) -> String + +Name of the Solution object. +""" +name(s::Solution) = get_string((n, b) -> LibCantera.sol_name(s.handle, n, b)) + +""" + transport_model(gas::Solution) -> String + +Name of the active transport model. +""" +transport_model(s::Solution) = + get_string((n, b) -> LibCantera.sol_transportModel(s.handle, n, b)) + +export Solution, close!, thermo, kinetics, transport, name, transport_model diff --git a/interfaces/julia/src/solutionarray.jl b/interfaces/julia/src/solutionarray.jl new file mode 100644 index 0000000000..0ad4e0862e --- /dev/null +++ b/interfaces/julia/src/solutionarray.jl @@ -0,0 +1,150 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# SolutionArray: a lightweight container for a batch of thermodynamic states of +# a single phase, mirroring the common uses of Python's `cantera.SolutionArray` +# (storing a reactor-integration history, a parameter sweep, or a flame profile +# and reading properties or writing CSV). +# +# Unlike the Python version this is a pure-Julia convenience: states are stored +# as (T, P, mass-fractions) and properties are evaluated by restoring each state +# into the backing `Solution` on demand. It needs no CLib support. + +""" + SolutionArray(gas::Solution, n::Integer) + SolutionArray(gas::Solution) + +A collection of `n` thermodynamic states of `gas` (or an empty, growable array). +States are stored as temperature, pressure, and mass fractions, and are restored +into `gas` on demand to evaluate properties. + +```julia +states = SolutionArray(gas) +for t in times + advance!(net, t) + append!(states, gas) # snapshot the current gas state +end +temperature(states) # Vector of T for every stored state +write_csv(states, "history.csv") +``` +""" +mutable struct SolutionArray + gas::Solution + T::Vector{Float64} + P::Vector{Float64} + Y::Matrix{Float64} # nSpecies x nStates +end + +function SolutionArray(gas::Solution, n::Integer) + nsp = n_species(gas) + T = fill(temperature(gas), n) + P = fill(pressure(gas), n) + Y = repeat(mass_fractions(gas), 1, n) + return SolutionArray(gas, T, P, Y) +end + +SolutionArray(gas::Solution) = + SolutionArray(gas, Float64[], Float64[], Matrix{Float64}(undef, n_species(gas), 0)) + +Base.length(sa::SolutionArray) = length(sa.T) +Base.size(sa::SolutionArray) = (length(sa),) +Base.isempty(sa::SolutionArray) = length(sa) == 0 + +"Restore stored state `i` (1-based) into the backing `Solution` and return it." +function restore!(sa::SolutionArray, i::Integer) + set_TPY!(sa.gas, sa.T[i], sa.P[i], @view sa.Y[:, i]) + return sa.gas +end + +""" + append!(sa::SolutionArray, gas::Solution) + +Append a snapshot of the current state of `gas` to `sa`. +""" +function Base.append!(sa::SolutionArray, gas::Solution) + push!(sa.T, temperature(gas)) + push!(sa.P, pressure(gas)) + sa.Y = hcat(sa.Y, mass_fractions(gas)) + return sa +end + +""" + set_state!(sa, i; T, P, X=nothing, Y=nothing) + +Set stored state `i` (1-based). Provide composition as mole fractions `X` or +mass fractions `Y` (string or vector); if neither is given the composition is +left unchanged. +""" +function set_state!(sa::SolutionArray, i::Integer; T, P, X=nothing, Y=nothing) + if X !== nothing + set_TPX!(sa.gas, T, P, X) + elseif Y !== nothing + set_TPY!(sa.gas, T, P, Y) + else + set_TP!(sa.gas, T, P) + end + sa.T[i] = temperature(sa.gas) + sa.P[i] = pressure(sa.gas) + sa.Y[:, i] = mass_fractions(sa.gas) + return sa +end + +# ---- vectorized property access -------------------------------------------- + +"Temperatures [K] of all stored states." +temperature(sa::SolutionArray) = copy(sa.T) + +"Pressures [Pa] of all stored states." +pressure(sa::SolutionArray) = copy(sa.P) + +"Mass fractions of all stored states, `nSpecies x nStates`." +mass_fractions(sa::SolutionArray) = copy(sa.Y) + +""" + extract(sa, f) -> Vector + +Evaluate the scalar property function `f` (e.g. `density`, `enthalpy_mass`) for +every stored state and return the results as a vector. +""" +function extract(sa::SolutionArray, f) + return [f(restore!(sa, i)) for i in 1:length(sa)] +end + +"Mass densities [kg/m^3] of all stored states." +density(sa::SolutionArray) = extract(sa, density) + +"Mole fractions of all stored states, `nSpecies x nStates`." +function mole_fractions(sa::SolutionArray) + nsp = n_species(sa.gas) + X = Matrix{Float64}(undef, nsp, length(sa)) + for i in 1:length(sa) + X[:, i] = mole_fractions(restore!(sa, i)) + end + return X +end + +""" + write_csv(sa::SolutionArray, path) + +Write the stored states to a CSV file with columns `T`, `P`, and one column per +species mass fraction. +""" +function write_csv(sa::SolutionArray, path::AbstractString) + names = species_names(sa.gas) + open(path, "w") do io + println(io, "T,P,", join(names, ",")) + for i in 1:length(sa) + print(io, sa.T[i], ",", sa.P[i]) + for k in 1:size(sa.Y, 1) + print(io, ",", sa.Y[k, i]) + end + println(io) + end + end + return path +end + +Base.show(io::IO, sa::SolutionArray) = + print(io, "SolutionArray(", length(sa), " states, ", n_species(sa.gas), " species)") + +export SolutionArray, restore!, set_state!, extract, write_csv diff --git a/interfaces/julia/src/thermo.jl b/interfaces/julia/src/thermo.jl new file mode 100644 index 0000000000..7d2371d76c --- /dev/null +++ b/interfaces/julia/src/thermo.jl @@ -0,0 +1,477 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Thermodynamic state, composition and property accessors. +# +# Every function accepts either a `Solution` or a `ThermoPhase` (collectively +# `ThermoLike`) and operates on the underlying ThermoPhase handle. Scalar +# getters return `Float64`; array getters return freshly allocated +# `Vector{Float64}` and have in-place `!` variants for hot loops. + +# ---- sizes and names -------------------------------------------------------- + +"Number of chemical species in the phase." +n_species(g::ThermoLike) = Int(check(LibCantera.thermo_nSpecies(_thermo_handle(g)))) + +"Number of elements in the phase." +n_elements(g::ThermoLike) = Int(check(LibCantera.thermo_nElements(_thermo_handle(g)))) + +"Name of species `k` (1-based)." +function species_name(g::ThermoLike, k::Integer) + h = _thermo_handle(g) + return get_string((n, b) -> LibCantera.thermo_speciesName(h, Int32(k - 1), n, b)) +end + +"Vector of all species names." +species_names(g::ThermoLike) = [species_name(g, k) for k in 1:n_species(g)] + +""" + species_index(gas, name) -> Int + +1-based index of species `name`, or `0` if it is not present. +""" +function species_index(g::ThermoLike, nm::AbstractString) + idx = LibCantera.thermo_speciesIndex(_thermo_handle(g), nm) + idx == -1 && return 0 # not found (CLib returns npos-as-(-1)) + return Int(idx) + 1 +end + +"Name of element `m` (1-based)." +function element_name(g::ThermoLike, m::Integer) + h = _thermo_handle(g) + return get_string((n, b) -> LibCantera.thermo_elementName(h, Int32(m - 1), n, b)) +end + +"Vector of all element names." +element_names(g::ThermoLike) = [element_name(g, m) for m in 1:n_elements(g)] + +# ---- scalar state ----------------------------------------------------------- + +"Temperature [K]." +temperature(g::ThermoLike) = checkd(LibCantera.thermo_temperature(_thermo_handle(g))) + +"Pressure [Pa]." +pressure(g::ThermoLike) = checkd(LibCantera.thermo_pressure(_thermo_handle(g))) + +"Mass density [kg/m^3]." +density(g::ThermoLike) = checkd(LibCantera.thermo_density(_thermo_handle(g))) + +"Molar density [kmol/m^3]." +molar_density(g::ThermoLike) = checkd(LibCantera.thermo_molarDensity(_thermo_handle(g))) + +"Mean molecular weight [kg/kmol]." +mean_molecular_weight(g::ThermoLike) = + checkd(LibCantera.thermo_meanMolecularWeight(_thermo_handle(g))) + +# ---- molar / specific properties ------------------------------------------- + +for (jl, c, doc) in ( + (:enthalpy_mass, :thermo_enthalpy_mass, "Specific enthalpy [J/kg]."), + (:enthalpy_mole, :thermo_enthalpy_mole, "Molar enthalpy [J/kmol]."), + (:internal_energy_mass, :thermo_intEnergy_mass, + "Specific internal energy [J/kg]."), + (:internal_energy_mole, :thermo_intEnergy_mole, + "Molar internal energy [J/kmol]."), + (:entropy_mass, :thermo_entropy_mass, "Specific entropy [J/kg/K]."), + (:entropy_mole, :thermo_entropy_mole, "Molar entropy [J/kmol/K]."), + (:gibbs_mass, :thermo_gibbs_mass, "Specific Gibbs free energy [J/kg]."), + (:gibbs_mole, :thermo_gibbs_mole, "Molar Gibbs free energy [J/kmol]."), + (:cp_mass, :thermo_cp_mass, + "Specific heat capacity at constant pressure [J/kg/K]."), + (:cp_mole, :thermo_cp_mole, + "Molar heat capacity at constant pressure [J/kmol/K]."), + (:cv_mass, :thermo_cv_mass, + "Specific heat capacity at constant volume [J/kg/K]."), + (:cv_mole, :thermo_cv_mole, + "Molar heat capacity at constant volume [J/kmol/K]."), + ) + @eval begin + $jl(g::ThermoLike) = checkd(LibCantera.$c(_thermo_handle(g))) + @doc $doc $jl + end +end + +"Isothermal compressibility [1/Pa]." +isothermal_compressibility(g::ThermoLike) = + checkd(LibCantera.thermo_isothermalCompressibility(_thermo_handle(g))) + +"Thermal (volumetric) expansion coefficient [1/K]." +thermal_expansion_coeff(g::ThermoLike) = + checkd(LibCantera.thermo_thermalExpansionCoeff(_thermo_handle(g))) + +"Maximum temperature [K] for which the phase's thermo data are valid." +max_temp(g::ThermoLike) = + checkd(LibCantera.thermo_maxTemp(_thermo_handle(g), Int32(-1))) + +"Minimum temperature [K] for which the phase's thermo data are valid." +min_temp(g::ThermoLike) = + checkd(LibCantera.thermo_minTemp(_thermo_handle(g), Int32(-1))) + +"Reference pressure [Pa] used for standard-state thermo data." +reference_pressure(g::ThermoLike) = + checkd(LibCantera.thermo_refPressure(_thermo_handle(g))) + +"Electric potential [V] of the phase." +electric_potential(g::ThermoLike) = + checkd(LibCantera.thermo_electricPotential(_thermo_handle(g))) + +"Set the electric potential [V] of the phase." +function set_electric_potential!(g::ThermoLike, v) + check(LibCantera.thermo_setElectricPotential(_thermo_handle(g), Float64(v))) + return g +end + +"Speed of sound [m/s] (`sqrt(cp/cv · p/ρ)`)." +sound_speed(g::ThermoLike) = sqrt(cp_mass(g) / cv_mass(g) * pressure(g) / density(g)) + +"Specific volume [m^3/kg]." +volume_mass(g::ThermoLike) = 1.0 / density(g) + +"Molar volume [m^3/kmol]." +volume_mole(g::ThermoLike) = mean_molecular_weight(g) / density(g) + +# ---- critical / saturation properties (raise for ideal-gas phases) ---------- + +"Critical temperature [K]." +critical_temperature(g::ThermoLike) = + checkd(LibCantera.thermo_critTemperature(_thermo_handle(g))) +"Critical pressure [Pa]." +critical_pressure(g::ThermoLike) = + checkd(LibCantera.thermo_critPressure(_thermo_handle(g))) +"Critical density [kg/m^3]." +critical_density(g::ThermoLike) = + checkd(LibCantera.thermo_critDensity(_thermo_handle(g))) +"Vapor fraction (quality) of a two-phase state." +vapor_fraction(g::ThermoLike) = + checkd(LibCantera.thermo_vaporFraction(_thermo_handle(g))) +"Saturation temperature [K] at pressure `p` [Pa]." +sat_temperature(g::ThermoLike, p) = + checkd(LibCantera.thermo_satTemperature(_thermo_handle(g), Float64(p))) +"Saturation pressure [Pa] at temperature `T` [K]." +sat_pressure(g::ThermoLike, T) = + checkd(LibCantera.thermo_satPressure(_thermo_handle(g), Float64(T))) + +# ---- elements / atoms ------------------------------------------------------- + +""" + element_index(gas, name) -> Int + +1-based index of element `name`, or `0` if it is not present. +""" +function element_index(g::ThermoLike, nm::AbstractString) + idx = LibCantera.thermo_elementIndex(_thermo_handle(g), nm) + idx == -1 && return 0 + return Int(idx) + 1 +end + +"Number of atoms of element `m` (1-based) in species `k` (1-based)." +n_atoms(g::ThermoLike, k::Integer, m::Integer) = + checkd(LibCantera.thermo_nAtoms(_thermo_handle(g), Int32(k - 1), Int32(m - 1))) + +"Atomic weights of all elements [kg/kmol]." +atomic_weights(g::ThermoLike) = + get_array(n_elements(g), + (n, b) -> LibCantera.thermo_atomicWeights(_thermo_handle(g), n, b)) + +"Atomic weight of element `m` (1-based) [kg/kmol]." +atomic_weight(g::ThermoLike, m::Integer) = atomic_weights(g)[m] + +"Species electric charges (per elementary charge)." +charges(g::ThermoLike) = + get_array(n_species(g), + (n, b) -> LibCantera.thermo_getCharges(_thermo_handle(g), n, b)) + +""" + elemental_mole_fraction(gas, m) -> Float64 + +Mole fraction of element `m` (1-based) — moles of that element per mole of all +elements in the mixture. +""" +function elemental_mole_fraction(g::ThermoLike, m::Integer) + X = mole_fractions(g); ne = n_elements(g); nsp = n_species(g) + num = sum(n_atoms(g, k, m) * X[k] for k in 1:nsp) + den = sum(sum(n_atoms(g, k, j) * X[k] for k in 1:nsp) for j in 1:ne) + return num / den +end + +""" + elemental_mass_fraction(gas, m) -> Float64 + +Mass fraction of element `m` (1-based) in the mixture. +""" +function elemental_mass_fraction(g::ThermoLike, m::Integer) + X = mole_fractions(g); nsp = n_species(g) + num = atomic_weight(g, m) * sum(n_atoms(g, k, m) * X[k] for k in 1:nsp) + return num / mean_molecular_weight(g) +end + +""" + equivalence_ratio(gas) -> Float64 + +Fuel/air equivalence ratio from the elemental composition (C, H, O, S), using +complete-combustion available-oxygen accounting: `φ = (2·Z_C + Z_H/2 + 2·Z_S) / Z_O`. +""" +function equivalence_ratio(g::ThermoLike) + Z(sym) = (i = element_index(g, sym); i == 0 ? 0.0 : elemental_mole_fraction(g, i)) + ZO = Z("O") + ZO == 0 && return Inf + return (2Z("C") + Z("H") / 2 + 2Z("S")) / ZO +end + +# ---- composition setters ---------------------------------------------------- + +""" + set_equivalence_ratio!(gas, phi, fuel, oxidizer) + +Set the mixture composition to the given equivalence ratio `phi` for the given +`fuel` and `oxidizer` compositions (name-value strings), holding T and p. +""" +function set_equivalence_ratio!(g::ThermoLike, phi, fuel::AbstractString, + oxidizer::AbstractString) + check(LibCantera.thermo_setEquivalenceRatio(_thermo_handle(g), Float64(phi), + fuel, oxidizer)) + return g +end + +# ---- standard-state (per species) ------------------------------------------- + +for (jl, c, doc) in ( + (:standard_enthalpies_RT, :thermo_getEnthalpy_RT, "enthalpies h°_k/RT"), + (:standard_entropies_R, :thermo_getEntropy_R, "entropies s°_k/R"), + (:standard_gibbs_RT, :thermo_getGibbs_RT, "Gibbs energies g°_k/RT"), + (:standard_int_energies_RT, :thermo_getIntEnergy_RT, + "internal energies u°_k/RT"), + (:standard_cp_R, :thermo_getCp_R, "heat capacities cp°_k/R"), + ) + @eval begin + @doc $("Nondimensional standard-state species " * doc * ".") + $jl(g::ThermoLike) = + get_array(n_species(g), (n, b) -> LibCantera.$c(_thermo_handle(g), n, b)) + end +end + +# ---- composition ------------------------------------------------------------ + +"Molecular weights of all species [kg/kmol]." +molecular_weights(g::ThermoLike) = + get_array(n_species(g), + (n, b) -> LibCantera.thermo_getMolecularWeights(_thermo_handle(g), n, b)) + +"In-place [`molecular_weights`](@ref)." +molecular_weights!(out, g::ThermoLike) = + get_array!(out, + (n, b) -> LibCantera.thermo_getMolecularWeights(_thermo_handle(g), n, b)) + +"Mole fractions of all species." +mole_fractions(g::ThermoLike) = + get_array(n_species(g), + (n, b) -> LibCantera.thermo_getMoleFractions(_thermo_handle(g), n, b)) + +"In-place [`mole_fractions`](@ref)." +mole_fractions!(out, g::ThermoLike) = + get_array!(out, + (n, b) -> LibCantera.thermo_getMoleFractions(_thermo_handle(g), n, b)) + +"Mass fractions of all species." +mass_fractions(g::ThermoLike) = + get_array(n_species(g), + (n, b) -> LibCantera.thermo_getMassFractions(_thermo_handle(g), n, b)) + +"In-place [`mass_fractions`](@ref)." +mass_fractions!(out, g::ThermoLike) = + get_array!(out, + (n, b) -> LibCantera.thermo_getMassFractions(_thermo_handle(g), n, b)) + +"Species concentrations [kmol/m^3]." +concentrations(g::ThermoLike) = + get_array(n_species(g), + (n, b) -> LibCantera.thermo_getConcentrations(_thermo_handle(g), n, b)) + +# ---- partial molar / potentials -------------------------------------------- + +for (jl, c, doc) in ( + (:partial_molar_enthalpies, :thermo_getPartialMolarEnthalpies, + "Partial molar enthalpies of the species [J/kmol]."), + (:partial_molar_entropies, :thermo_getPartialMolarEntropies, + "Partial molar entropies of the species [J/kmol/K]."), + (:partial_molar_int_energies, :thermo_getPartialMolarIntEnergies, + "Partial molar internal energies of the species [J/kmol]."), + (:partial_molar_cp, :thermo_getPartialMolarCp, + "Partial molar heat capacities at constant pressure [J/kmol/K]."), + (:partial_molar_volumes, :thermo_getPartialMolarVolumes, + "Partial molar volumes of the species [m^3/kmol]."), + (:chemical_potentials, :thermo_getChemPotentials, + "Chemical potentials of the species [J/kmol]."), + (:electrochemical_potentials, :thermo_getElectrochemPotentials, + "Electrochemical potentials of the species [J/kmol]."), + ) + bang = Symbol(jl, :!) + bang_doc = "In-place [`$jl`](@ref)." + @eval begin + $jl(g::ThermoLike) = + get_array(n_species(g), (n, b) -> LibCantera.$c(_thermo_handle(g), n, b)) + @doc $doc $jl + $bang(out, g::ThermoLike) = + get_array!(out, (n, b) -> LibCantera.$c(_thermo_handle(g), n, b)) + @doc $bang_doc $bang + end +end + +# ---- composition setters ---------------------------------------------------- + +""" + set_mole_fractions!(gas, X) + +Set mole fractions from a numeric vector or a composition string +(e.g. `"CH4:1, O2:2"`). The values are normalized by Cantera. +""" +function set_mole_fractions!(g::ThermoLike, X::AbstractVector{<:Real}) + x = as_f64(X) + check(LibCantera.thermo_setMoleFractions(_thermo_handle(g), Int32(length(x)), + pointer(x))) + return g +end +function set_mole_fractions!(g::ThermoLike, X::AbstractString) + check(LibCantera.thermo_setMoleFractionsByName(_thermo_handle(g), X)) + return g +end + +""" + set_mass_fractions!(gas, Y) + +Set mass fractions from a numeric vector or a composition string. +""" +function set_mass_fractions!(g::ThermoLike, Y::AbstractVector{<:Real}) + y = as_f64(Y) + check(LibCantera.thermo_setMassFractions(_thermo_handle(g), Int32(length(y)), + pointer(y))) + return g +end +function set_mass_fractions!(g::ThermoLike, Y::AbstractString) + check(LibCantera.thermo_setMassFractionsByName(_thermo_handle(g), Y)) + return g +end + +# ---- state setters ---------------------------------------------------------- + +"Set temperature [K] and pressure [Pa]." +function set_TP!(g::ThermoLike, T, P) + check(LibCantera.thermo_setState_TP(_thermo_handle(g), Float64(T), Float64(P))) + return g +end + +""" + set_TPX!(gas, T, p, X) + +Set temperature, pressure and mole fractions in one call. `X` may be a numeric +vector or a composition string. +""" +function set_TPX!(g::ThermoLike, T, P, X::AbstractVector{<:Real}) + x = as_f64(X) + check(LibCantera.thermo_setState_TPX(_thermo_handle(g), Float64(T), Float64(P), + Int32(length(x)), pointer(x))) + return g +end +function set_TPX!(g::ThermoLike, T, P, X::AbstractString) + check(LibCantera.thermo_setState_TPX_byName(_thermo_handle(g), Float64(T), + Float64(P), X)) + return g +end + +""" + set_TPY!(gas, T, p, Y) + +Set temperature, pressure and mass fractions in one call. `Y` may be a numeric +vector or a composition string. +""" +function set_TPY!(g::ThermoLike, T, P, Y::AbstractVector{<:Real}) + y = as_f64(Y) + check(LibCantera.thermo_setState_TPY(_thermo_handle(g), Float64(T), Float64(P), + Int32(length(y)), pointer(y))) + return g +end +function set_TPY!(g::ThermoLike, T, P, Y::AbstractString) + check(LibCantera.thermo_setState_TPY_byName(_thermo_handle(g), Float64(T), + Float64(P), Y)) + return g +end + +# Two-property setters that map directly onto CLib `setState_XY`. +for (jl, c, a, b, doc) in ( + (:set_HP!, :thermo_setState_HP, :h, :p, + "Set specific enthalpy `h` [J/kg] and pressure `p` [Pa]."), + (:set_UV!, :thermo_setState_UV, :u, :v, + "Set specific internal energy `u` [J/kg] and specific volume `v` [m^3/kg]."), + (:set_SP!, :thermo_setState_SP, :s, :p, + "Set specific entropy `s` [J/kg/K] and pressure `p` [Pa]."), + (:set_SV!, :thermo_setState_SV, :s, :v, + "Set specific entropy `s` [J/kg/K] and specific volume `v` [m^3/kg]."), + (:set_TD!, :thermo_setState_TD, :T, :rho, + "Set temperature `T` [K] and density `rho` [kg/m^3]."), + (:set_DP!, :thermo_setState_DP, :rho, :p, + "Set density `rho` [kg/m^3] and pressure `p` [Pa]."), + ) + @eval begin + function $jl(g::ThermoLike, $a, $b) + check(LibCantera.$c(_thermo_handle(g), Float64($a), Float64($b))) + return g + end + @doc $doc $jl + end +end + +""" + equilibrate!(gas, XY; solver="auto", rtol=1e-9, max_steps=1000, + max_iter=100, estimate_equil=0) + +Set the phase to a state of chemical equilibrium holding the property pair `XY` +(e.g. `"TP"`, `"HP"`, `"UV"`) fixed. +""" +function equilibrate!(g::ThermoLike, XY::AbstractString; solver="auto", + rtol=1e-9, max_steps=1000, max_iter=100, estimate_equil=0) + check(LibCantera.thermo_equilibrate(_thermo_handle(g), XY, solver, Float64(rtol), + Int32(max_steps), Int32(max_iter), + Int32(estimate_equil))) + return g +end + +""" + report(gas; show_thermo=true, threshold=1e-14) -> String + +Return Cantera's formatted state report for the phase. +""" +function report(g::ThermoLike; show_thermo::Bool=true, threshold=1e-14) + h = _thermo_handle(g) + return get_string((n, b) -> LibCantera.thermo_report(h, Int32(show_thermo), + Float64(threshold), n, b)) +end + +export report + +export n_species, n_elements, species_name, species_names, species_index, + element_name, element_names, + mole_fractions, mole_fractions!, mass_fractions, mass_fractions!, + set_mole_fractions!, set_mass_fractions!, molecular_weights, + molecular_weights!, concentrations, mean_molecular_weight + +export temperature, pressure, density, molar_density, + enthalpy_mass, enthalpy_mole, internal_energy_mass, internal_energy_mole, + entropy_mass, entropy_mole, gibbs_mass, gibbs_mole, + cp_mass, cp_mole, cv_mass, cv_mole, + isothermal_compressibility, thermal_expansion_coeff, + reference_pressure, electric_potential, set_electric_potential!, + standard_enthalpies_RT, standard_entropies_R, standard_gibbs_RT, + standard_int_energies_RT, standard_cp_R, + sound_speed, volume_mass, volume_mole, + critical_temperature, critical_pressure, critical_density, + vapor_fraction, sat_temperature, sat_pressure, + element_index, atomic_weights, atomic_weight, charges, + elemental_mole_fraction, elemental_mass_fraction, equivalence_ratio, + partial_molar_enthalpies, partial_molar_entropies, + partial_molar_int_energies, partial_molar_cp, partial_molar_volumes, + chemical_potentials, electrochemical_potentials, + partial_molar_enthalpies!, partial_molar_entropies!, partial_molar_cp!, + chemical_potentials! + +export set_TP!, set_TPX!, set_TPY!, set_HP!, set_UV!, set_SP!, set_SV!, + set_TD!, set_DP!, equilibrate!, set_equivalence_ratio! diff --git a/interfaces/julia/src/transport.jl b/interfaces/julia/src/transport.jl new file mode 100644 index 0000000000..e92ec39d93 --- /dev/null +++ b/interfaces/julia/src/transport.jl @@ -0,0 +1,61 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Transport-property accessors. Accept a `Solution` or a `Transport`. + +"Dynamic (mixture) viscosity [Pa*s]." +viscosity(g::TransportLike) = checkd(LibCantera.trans_viscosity(_transport_handle(g))) + +"Mixture thermal conductivity [W/m/K]." +thermal_conductivity(g::TransportLike) = + checkd(LibCantera.trans_thermalConductivity(_transport_handle(g))) + +"Electrical conductivity [S/m]." +electrical_conductivity(g::TransportLike) = + checkd(LibCantera.trans_electricalConductivity(_transport_handle(g))) + +"Mixture-averaged diffusion coefficients [m^2/s], one per species." +mix_diff_coeffs(g::Solution) = + get_array(n_species(g), + (n, b) -> LibCantera.trans_getMixDiffCoeffs(_transport_handle(g), n, b)) + +"In-place [`mix_diff_coeffs`](@ref)." +mix_diff_coeffs!(out, g::TransportLike) = + get_array!(out, + (n, b) -> LibCantera.trans_getMixDiffCoeffs(_transport_handle(g), n, b)) + +"Thermal diffusion coefficients [kg/m/s], one per species." +thermal_diff_coeffs(g::Solution) = + get_array(n_species(g), + (n, b) -> LibCantera.trans_getThermalDiffCoeffs( + _transport_handle(g), n, b)) + +""" + binary_diff_coeffs(gas) -> Matrix{Float64} + +`nsp x nsp` matrix of binary diffusion coefficients [m^2/s]. +""" +function binary_diff_coeffs(g::Solution) + nsp = n_species(g) + buf = Vector{Float64}(undef, nsp * nsp) + check(LibCantera.trans_getBinaryDiffCoeffs(_transport_handle(g), Int32(nsp), + Int32(length(buf)), pointer(buf))) + return reshape(buf, nsp, nsp) +end + +""" + multi_diff_coeffs(gas) -> Matrix{Float64} + +`nsp x nsp` matrix of multicomponent diffusion coefficients [m^2/s]. +""" +function multi_diff_coeffs(g::Solution) + nsp = n_species(g) + buf = Vector{Float64}(undef, nsp * nsp) + check(LibCantera.trans_getMultiDiffCoeffs(_transport_handle(g), Int32(nsp), + Int32(length(buf)), pointer(buf))) + return reshape(buf, nsp, nsp) +end + +export viscosity, thermal_conductivity, electrical_conductivity, + mix_diff_coeffs, mix_diff_coeffs!, thermal_diff_coeffs, + binary_diff_coeffs, multi_diff_coeffs diff --git a/interfaces/julia/src/utils.jl b/interfaces/julia/src/utils.jl new file mode 100644 index 0000000000..d65b407422 --- /dev/null +++ b/interfaces/julia/src/utils.jl @@ -0,0 +1,54 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Library-level utility functions and physical constants. + +"One standard atmosphere [Pa]." +const one_atm = 101325.0 + +"Universal gas constant [J/kmol/K] (Cantera's molar convention)." +const gas_constant = 8314.46261815324 + +"Avogadro's number [1/kmol]." +const avogadro = 6.02214076e26 + +""" + cantera_version() -> String + +Version string of the linked Cantera library. +""" +cantera_version() = get_string((n, b) -> LibCantera.ct_version(n, b)) + +""" + git_commit() -> String + +Git commit of the linked Cantera library. +""" +git_commit() = get_string((n, b) -> LibCantera.ct_gitCommit(n, b)) + +""" + add_data_directory(dir) + +Add `dir` to the list of directories searched for Cantera input files. +""" +function add_data_directory(dir::AbstractString) + check(LibCantera.ct_addDataDirectory(dir)) + return nothing +end + +""" + suppress_thermo_warnings(flag=true) + +Suppress (or re-enable) Cantera's thermodynamic warnings. +""" +function suppress_thermo_warnings(flag::Bool=true) + check(LibCantera.ct_suppressThermoWarnings(Int32(flag))) + return nothing +end + +"Reset all Cantera CLib storage, invalidating every existing handle." +reset_storage() = (check(LibCantera.ct_resetStorage()); nothing) + +export one_atm, gas_constant, avogadro +export cantera_version, git_commit, add_data_directory, suppress_thermo_warnings, + reset_storage diff --git a/interfaces/julia/test/reference_values.md b/interfaces/julia/test/reference_values.md new file mode 100644 index 0000000000..601503dda3 --- /dev/null +++ b/interfaces/julia/test/reference_values.md @@ -0,0 +1,41 @@ +# Regenerating test reference values + +The numeric reference values hard-coded in `runtests.jl` are produced with the +**Python** Cantera interface (the reference implementation), so that the Julia +interface is validated against Cantera's own results. + +They were generated with Cantera 3.2.0 and `gri30.yaml`. To regenerate: + +```python +import cantera as ct + +gas = ct.Solution("gri30.yaml") +gas.TPX = 1200.0, ct.one_atm, "CH4:1.0, O2:2.0, N2:7.52" +print("density ", repr(gas.density)) +print("cp_mass ", repr(gas.cp_mass)) +print("cv_mass ", repr(gas.cv_mass)) +print("enthalpy_mass ", repr(gas.enthalpy_mass)) +print("entropy_mass ", repr(gas.entropy_mass)) +print("mean_mw ", repr(gas.mean_molecular_weight)) +print("viscosity ", repr(gas.viscosity)) +print("conductivity ", repr(gas.thermal_conductivity)) +print("wdot(CH4) ", repr(gas.net_production_rates[gas.species_index("CH4")])) +print("kf[0] ", repr(gas.forward_rate_constants[0])) +print("eq[0] ", gas.reaction(0).equation) + +# constant-volume H2/air ignition reference +g2 = ct.Solution("gri30.yaml") +g2.TPX = 1000.0, ct.one_atm, "H2:2, O2:1, N2:4" +r = ct.IdealGasReactor(g2, clone=False) +net = ct.ReactorNet([r]) +net.advance(1e-3) +print("reactor T@1ms ", repr(r.T)) +``` + +Notes: + +* Python species indices are 0-based; the Julia interface is 1-based, so + `gas.species_index("OH") == 4` in Python corresponds to `species_index(gas, + "OH") == 5` in Julia. +* Tolerances in `runtests.jl` are chosen to accommodate benign floating-point + differences between the C++ core called directly (Julia) and through Python. diff --git a/interfaces/julia/test/runtests.jl b/interfaces/julia/test/runtests.jl new file mode 100644 index 0000000000..48e7b1215c --- /dev/null +++ b/interfaces/julia/test/runtests.jl @@ -0,0 +1,204 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + +# Reference values were generated with the Python Cantera interface (v3.2.0) for +# gri30.yaml and are checked in below. See test/reference_values.md for the +# script that regenerates them. + +const MECH = "gri30.yaml" + +@testset "Cantera.jl" begin + + @testset "Library loading & version" begin + v = cantera_version() + @test v isa String + @test !isempty(v) + @test occursin(r"^\d+\.\d+", v) + end + + @testset "Solution construction & finalization" begin + gas = Solution(MECH) + @test n_species(gas) == 53 + @test n_reactions(gas) == 325 + gas2 = Solution(MECH, "gri30") + @test n_species(gas2) == 53 + # close! is idempotent and must not crash. + close!(gas) + close!(gas) + @test occursin("closed", sprint(show, gas)) + close!(gas2) + end + + @testset "Sub-phase views keep their Solution alive" begin + # The views returned by `thermo`/`kinetics`/`transport` borrow handles + # owned by the Solution, so they hold a reference to it: collecting the + # Solution here would run its finalizer and free the handles the views + # are still using. Each view is built from its own Solution and is the + # only thing referencing it, otherwise the views would keep each other's + # owner alive and a regression in any one of them would go unnoticed. + view(f) = let gas = Solution(MECH) + set_TPX!(gas, 1200.0, one_atm, "CH4:1, O2:2, N2:7.52") + f(gas) + end + + tp = view(thermo) + kin = view(kinetics) + tr = view(transport) + GC.gc(); GC.gc() + + for v in (tp, kin, tr) + @test v.parent isa Solution + @test !v.parent.closed + end + @test temperature(tp) ≈ 1200.0 + @test n_reactions(kin) == 325 + @test viscosity(tr) > 0 + end + + @testset "Error handling" begin + @test_throws CanteraError Solution("this_file_does_not_exist.yaml") + gas = Solution(MECH) + # An unknown species name yields index 0 (not found), not an error. + @test species_index(gas, "NOTASPECIES") == 0 + close!(gas) + end + + @testset "Thermo properties vs reference" begin + gas = Solution(MECH) + set_TPX!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") + + @test temperature(gas) ≈ 1200.0 + @test pressure(gas) ≈ one_atm rtol=1e-10 + @test density(gas) ≈ 0.2806317906177915 rtol=1e-8 + @test cp_mass(gas) ≈ 1397.2506879911464 rtol=1e-8 + @test cv_mass(gas) ≈ 1096.3671002332637 rtol=1e-8 + @test enthalpy_mass(gas) ≈ 861934.8781246373 rtol=1e-8 + @test entropy_mass(gas) ≈ 8914.227317028952 rtol=1e-8 + @test mean_molecular_weight(gas) ≈ 27.633486692015207 rtol=1e-8 + close!(gas) + end + + @testset "Composition & species" begin + gas = Solution(MECH) + names = species_names(gas) + @test length(names) == 53 + @test names[1] == "H2" + @test "CH4" in names + @test species_index(gas, "OH") == 5 # 1-based + @test species_name(gas, 5) == "OH" + + mw = molecular_weights(gas) + @test length(mw) == 53 + @test mw[species_index(gas, "CH4")] ≈ 16.043 rtol=1e-4 + + # setters via string and via vector agree + set_TPX!(gas, 300.0, one_atm, "CH4:1, O2:2") + X = mole_fractions(gas) + @test sum(X) ≈ 1.0 rtol=1e-10 + gas2 = Solution(MECH) + set_TP!(gas2, 300.0, one_atm) + set_mole_fractions!(gas2, X) + @test mole_fractions(gas2) ≈ X rtol=1e-10 + + # mass-fraction round trip + Y = mass_fractions(gas) + @test sum(Y) ≈ 1.0 rtol=1e-10 + close!(gas); close!(gas2) + end + + @testset "Kinetics vs reference" begin + gas = Solution(MECH) + set_TPX!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") + + eqs = reaction_equations(gas) + @test length(eqs) == 325 + @test eqs[1] == "2 O + M <=> O2 + M" + + w = net_production_rates(gas) + @test length(w) == 53 + @test w[species_index(gas, "CH4")] ≈ -8.650790742316334e-6 rtol=1e-6 + + kf = forward_rate_constants(gas) + @test length(kf) == 325 + @test kf[1] ≈ 1.0e8 rtol=1e-8 + + # net = fwd - rev rates of progress + rop_net = net_rates_of_progress(gas) + rop_f = forward_rates_of_progress(gas) + rop_r = reverse_rates_of_progress(gas) + @test rop_net ≈ rop_f .- rop_r rtol=1e-10 + + # net production = creation - destruction + @test net_production_rates(gas) ≈ + creation_rates(gas) .- destruction_rates(gas) rtol=1e-8 + close!(gas) + end + + @testset "Transport vs reference" begin + gas = Solution(MECH) + set_TPX!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") + @test transport_model(gas) == "mixture-averaged" + @test viscosity(gas) ≈ 4.6872182846913704e-5 rtol=1e-6 + @test thermal_conductivity(gas) ≈ 0.08965525008638553 rtol=1e-6 + + d = mix_diff_coeffs(gas) + @test length(d) == 53 + @test all(d .> 0) + B = binary_diff_coeffs(gas) + @test size(B) == (53, 53) + close!(gas) + end + + @testset "In-place array getters" begin + gas = Solution(MECH) + set_TPX!(gas, 1500.0, 2 * one_atm, "CH4:1, O2:2, N2:7.52") + out = Vector{Float64}(undef, n_species(gas)) + net_production_rates!(out, gas) + @test out == net_production_rates(gas) + mole_fractions!(out, gas) + @test out ≈ mole_fractions(gas) + close!(gas) + end + + @testset "Reactor time integration vs reference" begin + gas = Solution(MECH) + set_TPX!(gas, 1000.0, one_atm, "H2:2, O2:1, N2:4") + r = IdealGasReactor(gas) + net = ReactorNet(r) + @test time(net) ≈ 0.0 + advance!(net, 1e-3) + @test time(net) ≈ 1e-3 + # ignition to equilibrium-ish temperature + @test temperature(r) ≈ 2869.724616612829 rtol=1e-3 + close!(net); close!(r); close!(gas) + end + + @testset "Analytic kinetic derivatives" begin + gas = Solution(MECH) + set_TPX!(gas, 1400.0, one_atm, "CH4:1, O2:2, N2:7.52") + d = net_production_rates_ddT(gas) + @test length(d) == 53 + @test all(isfinite, d) + # value check against the Python Cantera reference + @test d[species_index(gas, "CH4")] ≈ -3.3322745901816125e-6 rtol=1e-6 + # in-place variant agrees + out = similar(d); net_production_rates_ddT!(out, gas) + @test out == d + close!(gas) + end + + include("test_mechanisms.jl") + include("test_thermo.jl") + include("test_kinetics.jl") + include("test_transport.jl") + include("test_reactor.jl") + include("test_solutionarray.jl") + include("test_func1.jl") + include("test_multiphase.jl") + include("test_connectors.jl") + include("test_rdiag.jl") + include("test_onedim.jl") +end diff --git a/interfaces/julia/test/test_connectors.jl b/interfaces/julia/test/test_connectors.jl new file mode 100644 index 0000000000..5408b41e98 --- /dev/null +++ b/interfaces/julia/test/test_connectors.jl @@ -0,0 +1,105 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Tests for reactor connectors, reservoirs, surfaces and sensitivities. + +using Cantera, Test + +# Names not necessarily exported yet are reached through the module. +import Cantera: Reservoir, Wall, MassFlowController, Valve, PressureController, + ReactorSurface, mass_flow_rate, set_mass_flow_rate!, + device_coefficient, set_device_coefficient!, + heat_rate, expansion_rate, area, set_area!, + set_heat_transfer_coeff!, add_sensitivity_reaction!, + n_sens_params, sensitivity, rtol, atol, + set_sensitivity_tolerances!, connector_type, name + +const CONN_MECH = "h2o2.yaml" + +@testset "MassFlowController" begin + gas = Solution(CONN_MECH) + set_TPX!(gas, 1000.0, one_atm, "H2:2,O2:1") + gas2 = Solution(CONN_MECH) + set_TPX!(gas2, 1000.0, one_atm, "H2:2,O2:1") + up = Reservoir(gas2) + r = IdealGasReactor(gas) + mfc = MassFlowController(up, r; mdot=0.1) + net = ReactorNet([r]) + advance!(net, 0.0) # initialize the network so the device is ready + @test mass_flow_rate(mfc) ≈ 0.1 + @test device_coefficient(mfc) ≈ 0.1 + @test connector_type(mfc) == "MassFlowController" +end + +@testset "Wall heat transfer" begin + g1 = Solution(CONN_MECH); set_TPX!(g1, 1000.0, one_atm, "H2:2,O2:1") + g2 = Solution(CONN_MECH); set_TPX!(g2, 300.0, one_atm, "H2:2,O2:1") + r1 = IdealGasReactor(g1); r2 = IdealGasReactor(g2) + set_initial_volume!(r1, 1.0); set_initial_volume!(r2, 1.0) + set_chemistry_enabled!(r1, false); set_chemistry_enabled!(r2, false) + w = Wall(r1, r2; A=1.0, U=100.0) + @test area(w) ≈ 1.0 + net = ReactorNet([r1, r2]) + T1_0 = temperature(r1); T2_0 = temperature(r2) + advance!(net, 1.0) + T1 = temperature(r1); T2 = temperature(r2) + # temperatures move toward each other + @test T1 < T1_0 + @test T2 > T2_0 + @test T1 > T2 + q = heat_rate(w) + @test isfinite(q) && q > 0 +end + +@testset "Valve" begin + g1 = Solution(CONN_MECH); set_TPX!(g1, 500.0, 3 * one_atm, "H2:2,O2:1") + g2 = Solution(CONN_MECH); set_TPX!(g2, 500.0, 1 * one_atm, "H2:2,O2:1") + r1 = IdealGasReactor(g1); r2 = IdealGasReactor(g2) + set_initial_volume!(r1, 1.0); set_initial_volume!(r2, 1.0) + set_chemistry_enabled!(r1, false); set_chemistry_enabled!(r2, false) + K = 1e-5 + v = Valve(r1, r2; K=K) + @test connector_type(v) == "Valve" + net = ReactorNet([r1, r2]) + advance!(net, 0.0) + p1_0 = pressure(r1); p2_0 = pressure(r2) + @test mass_flow_rate(v) ≈ K * (p1_0 - p2_0) rtol=1e-6 # flow high -> low + @test mass_flow_rate(v) > 0 + advance!(net, 5.0) + @test pressure(r1) < p1_0 # high side drops + @test pressure(r2) > p2_0 # low side rises + @test pressure(r1) ≈ pressure(r2) rtol=1e-2 # pressures equilibrate +end + +@testset "PressureController" begin + up = Solution(CONN_MECH); set_TPX!(up, 800.0, one_atm, "H2:2,O2:1") + rg = Solution(CONN_MECH); set_TPX!(rg, 800.0, one_atm, "H2:2,O2:1") + down = Solution(CONN_MECH); set_TPX!(down, 800.0, one_atm, "H2:2,O2:1") + upstream = Reservoir(up) + r = IdealGasReactor(rg); set_chemistry_enabled!(r, false) + downstream = Reservoir(down) + + mfc = MassFlowController(upstream, r; mdot=0.05) + pc = PressureController(r, downstream; primary=mfc, K=1e-4) + @test connector_type(pc) == "PressureController" + net = ReactorNet([r]) + advance!(net, 0.0) + + @test mass_flow_rate(pc) ≈ mass_flow_rate(mfc) rtol=1e-6 + @test isfinite(mass_flow_rate(pc)) +end + +@testset "Sensitivity" begin + gas = Solution(CONN_MECH) + set_TPX!(gas, 1000.0, one_atm, "H2:2,O2:1") + r = IdealGasReactor(gas) + net = ReactorNet([r]) + add_sensitivity_reaction!(r, 1) # 1-based -> reaction 0 + @test n_sens_params(r) == 1 + set_sensitivity_tolerances!(net; rtol=1e-6, atol=1e-6) + advance!(net, 2e-4) + s = sensitivity(net, "temperature", 1, r) + @test isfinite(s) + @test rtol(net) > 0 + @test atol(net) > 0 +end diff --git a/interfaces/julia/test/test_func1.jl b/interfaces/julia/test/test_func1.jl new file mode 100644 index 0000000000..04f7bced90 --- /dev/null +++ b/interfaces/julia/test/test_func1.jl @@ -0,0 +1,52 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + +# Reference values cross-checked against Python `cantera.Func1` (v3.2.0). + +const CT = Cantera + +@testset "Func1" begin + @testset "basic evaluation" begin + f = CT.Func1("sin", 2.0) + @test f(pi / 4) ≈ sin(2 * pi / 4) + @test CT.evaluate(f, pi / 4) ≈ sin(2 * pi / 4) + @test CT.func_type(f) == "sin" + + c = CT.constant_function(3.5) + @test c(123.0) ≈ 3.5 + end + + @testset "advanced (polynomial)" begin + # polynomial3 with [1,2,3,4] -> t^3 + 2t^2 + 3t + 4 + p = CT.Func1("polynomial3", [1.0, 2.0, 3.0, 4.0]) + @test p(2.0) ≈ 26.0 # Python: 26.0 + @test p(0.0) ≈ 4.0 + end + + @testset "arithmetic" begin + s = CT.Func1("constant", 2.0) + CT.Func1("constant", 3.0) + @test s(1.0) ≈ 5.0 + d = CT.Func1("constant", 5.0) - CT.Func1("constant", 2.0) + @test d(1.0) ≈ 3.0 + pr = CT.Func1("constant", 2.0) * CT.Func1("constant", 4.0) + @test pr(1.0) ≈ 8.0 + r = CT.Func1("constant", 6.0) / CT.Func1("constant", 3.0) + @test r(1.0) ≈ 2.0 + end + + @testset "derivative" begin + # d/dt sin(w t) = w cos(w t); at t=0 -> w + f = CT.Func1("sin", 2.0) + df = CT.derivative(f) + @test df(0.0) ≈ 2.0 + end + + @testset "show / write" begin + f = CT.Func1("sin", 2.0) + @test CT.write(f) isa String + @test occursin("sin", sprint(show, f)) + end +end diff --git a/interfaces/julia/test/test_kinetics.jl b/interfaces/julia/test/test_kinetics.jl new file mode 100644 index 0000000000..faad897240 --- /dev/null +++ b/interfaces/julia/test/test_kinetics.jl @@ -0,0 +1,50 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test +using LinearAlgebra + +# Reference values from Python Cantera 3.2.0, gri30 at (1200 K, 1 atm, +# CH4:1, O2:2, N2:7.52), unless noted otherwise. +@testset "Kinetics (Python parity)" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") + + @testset "stoichiometry and rates" begin + kCH4 = species_index(gas, "CH4") + @test sum(reactant_stoich_coeffs(gas)[kCH4, :]) ≈ 6.0 + @test size(product_stoich_coeffs(gas)) == (n_species(gas), n_reactions(gas)) + @test delta_standard_enthalpy(gas)[1] ≈ -506662565.3788193 rtol=1e-8 + @test forward_rates_of_progress_ddT(gas)[1] ≈ 0.0 atol=1e-30 + @test multiplier(gas, 1) ≈ 1.0 + set_multiplier!(gas, 1, 2.0); @test multiplier(gas, 1) ≈ 2.0 + end + + @testset "heat release / production rates" begin + @test heat_release_rate(gas) ≈ -2068.555280246291 rtol=1e-8 + hpr = heat_production_rates(gas) + @test length(hpr) == n_reactions(gas) + # sum over reactions equals the total heat release rate + @test sum(hpr) ≈ heat_release_rate(gas) rtol=1e-10 + @test sum(hpr) ≈ -2068.5552802462907 rtol=1e-8 + end + + close!(gas) +end + +@testset "Composition Jacobian (_ddX)" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 1400.0, one_atm, "CH4:1, O2:2, N2:7.52") + + J = net_production_rates_ddX(gas) + @test size(J) == (53, 53) + @test J[species_index(gas, "CH4"), species_index(gas, "CH4")] ≈ + -0.0022867583538256886 rtol=1e-8 + @test norm(J) ≈ 7289124.658794 rtol=1e-8 + + R = net_rates_of_progress_ddX(gas) + @test size(R) == (n_reactions(gas), n_species(gas)) + @test all(isfinite, R) + close!(gas) +end diff --git a/interfaces/julia/test/test_mechanisms.jl b/interfaces/julia/test/test_mechanisms.jl new file mode 100644 index 0000000000..f433c651f4 --- /dev/null +++ b/interfaces/julia/test/test_mechanisms.jl @@ -0,0 +1,32 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + +# A second mechanism (h2o2.yaml, 10 species / 29 reactions) to complement the +# gri30-based reference checks and guard against mechanism-specific assumptions. +# Reference values from Python Cantera 3.x at 1000 K, 1 atm, H2:2/O2:1/N2:4. + +@testset "h2o2 mechanism vs reference" begin + gas = Solution("h2o2.yaml") + set_TPX!(gas, 1000.0, one_atm, "H2:2, O2:1, N2:4") + + @test n_species(gas) == 10 + @test n_reactions(gas) == 29 + @test density(gas) ≈ 0.25780918724920693 rtol=1e-8 + @test cp_mass(gas) ≈ 1527.8760405617072 rtol=1e-8 + @test enthalpy_mass(gas) ≈ 1012650.3485874726 rtol=1e-8 + @test viscosity(gas) ≈ 4.199970957871177e-5 rtol=1e-6 + @test thermal_conductivity(gas) ≈ 0.1317893909081434 rtol=1e-6 + + @test species_index(gas, "OH") == 5 # 1-based (Python 0-based 4) + @test reaction_equations(gas)[1] == "2 O + M <=> O2 + M" + + # consistency identities that must hold for any mechanism + @test net_production_rates(gas) ≈ + creation_rates(gas) .- destruction_rates(gas) rtol=1e-8 + @test net_rates_of_progress(gas) ≈ + forward_rates_of_progress(gas) .- reverse_rates_of_progress(gas) rtol=1e-10 + close!(gas) +end diff --git a/interfaces/julia/test/test_multiphase.jl b/interfaces/julia/test/test_multiphase.jl new file mode 100644 index 0000000000..dc60f18cd3 --- /dev/null +++ b/interfaces/julia/test/test_multiphase.jl @@ -0,0 +1,51 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + +# Reference values cross-checked against Python `cantera.Mixture` (v3.2.0) for a +# single gri30 gas phase at 1200 K, 1 atm, CH4:1 O2:2 N2:7.52. + +const CT = Cantera + +@testset "MultiPhase" begin + gas = Solution(MECH) + set_TPX!(gas, 1200.0, one_atm, "CH4:1, O2:2, N2:7.52") + + mix = CT.MultiPhase([gas => 1.0]) + @test CT.n_phases(mix) == 1 + @test CT.n_species(mix) == 53 + @test CT.n_elements(mix) == 5 + + CT.set_temperature!(mix, 1200.0) + CT.set_pressure!(mix, one_atm) + @test CT.temperature(mix) ≈ 1200.0 + @test CT.pressure(mix) ≈ one_atm rtol = 1e-10 + @test CT.phase_moles(mix, 1) ≈ 1.0 rtol = 1e-8 + + ci = CT.element_index(mix, "C") + @test ci > 0 + @test CT.element_moles(mix, ci) ≈ 0.09505703422053233 rtol = 1e-6 + + CT.equilibrate!(mix, "TP") + @test CT.temperature(mix) ≈ 1200.0 rtol = 1e-8 + @test CT.phase_moles(mix, 1) ≈ 1.0000017646000376 rtol = 1e-6 + + # gibbs == sum(mu_k * moles_k) + mu = CT.chemical_potentials(mix) + sm = [CT.species_moles(mix, k) for k in 1:CT.n_species(mix)] + @test CT.gibbs(mix) ≈ sum(mu .* sm) rtol = 1e-8 + @test CT.gibbs(mix) ≈ -347848546.95592207 rtol = 1e-5 + + @test isfinite(CT.enthalpy(mix)) + @test isfinite(CT.entropy(mix)) + @test isfinite(CT.cp(mix)) + @test isfinite(CT.volume(mix)) + @test CT.min_temp(mix) < CT.max_temp(mix) + + @test occursin("phases", sprint(show, mix)) + + CT.close!(mix) + close!(gas) +end diff --git a/interfaces/julia/test/test_onedim.jl b/interfaces/julia/test/test_onedim.jl new file mode 100644 index 0000000000..b603f9e271 --- /dev/null +++ b/interfaces/julia/test/test_onedim.jl @@ -0,0 +1,120 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + + +const CT = Cantera + +@testset "OneDim / FreeFlame" begin + gas = CT.Solution("gri30.yaml") + CT.set_TPX!(gas, 300.0, 101325.0, "CH4:1, O2:2, N2:7.52") + + flame = CT.FreeFlame(gas; width=0.03) + npts0 = CT.n_points(flame) + rho_u = flame.rho_u # unburned density [kg/m^3] + + @testset "construction" begin + @test CT.domain_type(flame.inlet) == "inlet" + @test CT.domain_type(flame.flow) == "free-flow" + @test CT.domain_type(flame.outlet) == "outlet" + comps = CT.component_names(flame.flow) + @test "T" in comps + @test "velocity" in comps + # Python's initial FreeFlame grid is 8 points; setting the fixed + # temperature inserts one anchor node, giving 9. + @test npts0 in (8, 9) + @test occursin("FreeFlame", sprint(show, flame)) + end + + CT.set_refine_criteria!(flame; ratio=3.0, slope=0.06, curve=0.12, prune=0.0) + + t0 = time() + CT.solve!(flame; loglevel=0, auto=true) + dt = time() - t0 + Su = CT.flame_speed(flame) + z = CT.grid(flame) + T = CT.flame_T(flame) + Tmax = maximum(T) + + ad = CT.Solution("gri30.yaml") + CT.set_TPX!(ad, 300.0, 101325.0, "CH4:1, O2:2, N2:7.52") + CT.equilibrate!(ad, "HP") + Tad = CT.temperature(ad) + rho_b = CT.density(ad) # burned density at the adiabatic state + + println("Su=$(round(Su; digits=4)) m/s Tmax=$(round(Tmax; digits=1)) K " * + "Tad=$(round(Tad; digits=1)) K npts=$(length(z)) " * + "z_end=$(round(z[end]; digits=4)) m solve=$(round(dt; digits=1)) s") + + @testset "flame is physical (behavioral)" begin + # Burned gas reaches the adiabatic flame temperature (Python: rel=2e-2). + @test isapprox(Tmax, Tad; rtol=2e-2) + # Laminar flame speed in the physically expected range (Python asserts + # flame speed to rel=1e-1 for its mixtures). + @test isapprox(Su, 0.38; rtol=0.15) + # Mass flux ρu is conserved across the flame: ρ_b·u_b ≈ ρ_u·S_u. + u_b = CT.flame_velocity(flame)[end] + @test isapprox(rho_b * u_b, rho_u * Su; rtol=0.05) + # Grid is a valid, refined, strictly increasing mesh. + @test all(diff(z) .> 0) + @test length(z) > npts0 + # Cold reactants, hot products. + @test T[1] < 400.0 + @test Tmax > 1800.0 + # Fuel consumed, products formed across the front. + Xch4 = CT.flame_X(flame, "CH4") + Xco2 = CT.flame_X(flame, "CO2") + @test Xch4[1] > Xch4[end] + @test Xco2[end] > Xco2[1] + end + + CT.close!(ad) + CT.close!(flame) + @test occursin("closed", sprint(show, flame)) + CT.close!(gas) +end + +@testset "OneDim / BurnerFlame" begin + gas = CT.Solution("gri30.yaml") + CT.set_TPX!(gas, 300.0, 101325.0, "CH4:0.9, O2:2, N2:7.52") # lean CH4/air + + mdot = 0.06 # kg/m^2/s (burner-stabilized flat flame) + bf = CT.BurnerFlame(gas; width=0.03, mdot=mdot) + + @testset "construction" begin + @test CT.domain_type(bf.burner) == "inlet" + @test CT.domain_type(bf.flow) == "unstrained-flow" + @test CT.domain_type(bf.outlet) == "outlet" + comps = CT.component_names(bf.flow) + @test "T" in comps + @test "velocity" in comps + @test CT.n_points(bf) == 7 # [0,0.1,0.2,0.3,0.5,0.7,1.0]*width + @test CT.burner_mdot(bf) ≈ mdot + @test occursin("BurnerFlame", sprint(show, bf)) + end + + CT.solve!(bf; loglevel=0, auto=false, refine_grid=false) + z = CT.grid(bf) + T = CT.flame_T(bf) + + @testset "coarse solve (smoke)" begin + @test all(diff(z) .> 0) # strictly increasing grid + @test z[1] ≈ 0.0 atol = 1e-9 + @test length(T) == length(z) + @test all(isfinite, T) + @test T[1] < 400.0 # cold burner side + @test maximum(T) > 1500.0 # hot downstream + Xch4 = CT.flame_X(bf, "CH4") + Xco2 = CT.flame_X(bf, "CO2") + @test all(isfinite, Xch4) + @test all(isfinite, Xco2) + @test Xch4[1] > Xch4[end] # fuel consumed downstream + @test Xco2[end] > Xco2[1] # products formed downstream + end + + CT.close!(bf) + @test occursin("closed", sprint(show, bf)) + CT.close!(gas) +end diff --git a/interfaces/julia/test/test_rdiag.jl b/interfaces/julia/test/test_rdiag.jl new file mode 100644 index 0000000000..d7d18cbc77 --- /dev/null +++ b/interfaces/julia/test/test_rdiag.jl @@ -0,0 +1,45 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera, Test + +@testset "ReactionPathDiagram" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 2000.0, one_atm, "CH4:1, O2:2, N2:7.52") + + d = ReactionPathDiagram(gas, "C") + set_threshold!(d, 0.01) + set_flow_type!(d, "NetFlow") + build!(d) + + dot = get_dot(d) + @test occursin("digraph", dot) # valid graphviz + @test !isempty(get_data(d)) + @test threshold(d) ≈ 0.01 + @test flow_type(d) == "NetFlow" + + # a few more property round-trips + set_bold_threshold!(d, 0.5) + @test bold_threshold(d) ≈ 0.5 + set_scale!(d, 2.0) + @test scale(d) ≈ 2.0 + set_arrow_width!(d, -1.0) + @test arrow_width(d) ≈ -1.0 + set_title!(d, "carbon flow") + @test title(d) == "carbon flow" + set_font!(d, "Courier") + @test font(d) == "Courier" + set_show_details!(d, true) + @test show_details(d) == true + set_bold_color!(d, "red") + @test bold_color(d) == "red" + + # display_only!/all should not error + display_only!(d, 1) + display_only!(d, :all) + build!(d) + @test occursin("digraph", get_dot(d)) + + close!(d) + @test occursin("closed", sprint(show, d)) +end diff --git a/interfaces/julia/test/test_reactor.jl b/interfaces/julia/test/test_reactor.jl new file mode 100644 index 0000000000..2b47612a41 --- /dev/null +++ b/interfaces/julia/test/test_reactor.jl @@ -0,0 +1,64 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + +@testset "ReactorNet state & components" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 1000.0, one_atm, "H2:2, O2:1, N2:4") + r = IdealGasReactor(gas) + net = ReactorNet(r) + + # mass, volume, temperature + species + @test n_components(net) == 3 + n_species(gas) + y = state(net) + @test length(y) == n_components(net) + @test all(isfinite, y) + + names = component_names(net) + @test length(names) == n_components(net) + @test occursin("mass", names[1]) + @test occursin("temperature", names[3]) + close!(net); close!(r); close!(gas) +end + +@testset "ConstPressureReactor ignition" begin + for ctor in (ConstPressureReactor, IdealGasConstPressureReactor) + gas = Solution("gri30.yaml") + set_TPX!(gas, 1200.0, one_atm, "H2:2, O2:1, N2:4") + p0 = pressure(gas) + + # Independent target: HP equilibrium of the same initial mixture. + ad = Solution("gri30.yaml") + set_TPX!(ad, 1200.0, one_atm, "H2:2, O2:1, N2:4") + equilibrate!(ad, "HP") + Tad = temperature(ad) + + r = ctor(gas) + net = ReactorNet([r]) + advance!(net, 1.0) # well past ignition + + @test pressure(r) ≈ p0 rtol=1e-4 # constant pressure + @test temperature(r) ≈ Tad rtol=2e-2 # adiabatic flame temperature + close!(net); close!(r); close!(ad); close!(gas) + end +end + +@testset "Reservoir holds state" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 800.0, one_atm, "CH4:1, O2:2, N2:7.52") + res = Reservoir(gas) + T0 = temperature(res); p0 = pressure(res) + + driver = Solution("gri30.yaml") + set_TPX!(driver, 1200.0, one_atm, "H2:2, O2:1, N2:4") + r = IdealGasReactor(driver) + MassFlowController(res, r; mdot=0.01) # reservoir feeds the reactor + net = ReactorNet([r]) + advance!(net, 0.5) + + @test temperature(res) ≈ T0 rtol=1e-10 # unchanged + @test pressure(res) ≈ p0 rtol=1e-10 + close!(net); close!(r); close!(res); close!(gas); close!(driver) +end diff --git a/interfaces/julia/test/test_solutionarray.jl b/interfaces/julia/test/test_solutionarray.jl new file mode 100644 index 0000000000..ddac3280fd --- /dev/null +++ b/interfaces/julia/test/test_solutionarray.jl @@ -0,0 +1,37 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + +@testset "SolutionArray" begin + gas = Solution("gri30.yaml") + # snapshot a short reactor history + set_TPX!(gas, 1200.0, one_atm, "H2:2, O2:1, N2:4") + r = IdealGasReactor(gas); net = ReactorNet(r) + states = SolutionArray(gas) + for t in range(0, 1e-3; length=5) + advance!(net, t) + append!(states, reactor_phase(r)) + end + @test length(states) == 5 + T = temperature(states) + @test length(T) == 5 + @test T[end] > T[1] + 500 # constant-volume ignition heats up + @test maximum(T) > 2000 + @test size(mass_fractions(states)) == (n_species(gas), 5) + @test all(density(states) .> 0) + @test extract(states, enthalpy_mass) isa Vector{Float64} + + # fixed-size construction + set_state! + CSV round-trip + sa = SolutionArray(gas, 3) + @test length(sa) == 3 + set_state!(sa, 1; T=300.0, P=one_atm, X="CH4:1, O2:2") + @test temperature(sa)[1] ≈ 300.0 + path = tempname() * ".csv" + write_csv(sa, path) + @test isfile(path) + @test occursin("T,P,", readline(path)) + rm(path; force=true) + close!(net); close!(r); close!(gas) +end diff --git a/interfaces/julia/test/test_thermo.jl b/interfaces/julia/test/test_thermo.jl new file mode 100644 index 0000000000..f802f2279b --- /dev/null +++ b/interfaces/julia/test/test_thermo.jl @@ -0,0 +1,131 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test + +# Reference values from Python Cantera 3.2.0, gri30 at (1200 K, 1 atm, +# CH4:1, O2:2, N2:7.52), unless noted otherwise. +@testset "ThermoPhase (Python parity)" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") + + @testset "scalar properties" begin + @test isothermal_compressibility(gas) ≈ 9.869232667160127e-06 rtol=1e-8 + @test thermal_expansion_coeff(gas) ≈ 0.0008333333333333334 rtol=1e-8 + @test max_temp(gas) ≈ 3000.0 + @test min_temp(gas) ≈ 300.0 + @test sound_speed(gas) ≈ 678.3425211864376 rtol=1e-8 + @test volume_mass(gas) ≈ 3.5633881599749238 rtol=1e-10 + @test volume_mole(gas) ≈ 98.46883929715162 rtol=1e-10 + @test reference_pressure(gas) ≈ 101325.0 + @test electric_potential(gas) ≈ 0.0 + @test equivalence_ratio(gas) ≈ 1.0 rtol=1e-10 + end + + @testset "elements and atoms" begin + iC = element_index(gas, "C") + @test elemental_mass_fraction(gas, iC) ≈ 0.04131690114779184 rtol=1e-10 + @test elemental_mole_fraction(gas, iC) ≈ 0.0415973377703827 rtol=1e-10 + @test atomic_weights(gas) ≈ [15.999, 1.008, 12.011, 14.007, 39.95] rtol=1e-8 + @test all(charges(gas) .== 0.0) + end + + @testset "set_equivalence_ratio! round-trip" begin + set_equivalence_ratio!(gas, 0.5, "CH4:1", "O2:1, N2:3.76") + @test equivalence_ratio(gas) ≈ 0.5 rtol=1e-8 + end + + @testset "standard-state properties" begin + set_TPX!(gas, 1200.0, one_atm, "CH4:1.0, O2:2.0, N2:7.52") + k = species_index(gas, "CH4") + @test standard_cp_R(gas)[k] ≈ 9.790762829272001 rtol=1e-8 + @test standard_enthalpies_RT(gas)[k] ≈ -2.0465605875189334 rtol=1e-8 + @test standard_entropies_R(gas)[k] ≈ 31.561123052390712 rtol=1e-8 + @test standard_gibbs_RT(gas)[k] ≈ -33.607683639909645 rtol=1e-8 + @test standard_int_energies_RT(gas)[k] ≈ -3.0465605875189334 rtol=1e-8 + # G/RT = H/RT - S/R identity + @test standard_gibbs_RT(gas)[k] ≈ + standard_enthalpies_RT(gas)[k] - standard_entropies_R(gas)[k] rtol=1e-10 + end + + close!(gas) +end + + +@testset "two-property state setters (round-trip)" begin + gas = Solution("gri30.yaml") + comp = "CH4:1.0, O2:2.0, N2:7.52" + set_TPX!(gas, 1200.0, one_atm, comp) + T0 = temperature(gas); p0 = pressure(gas); rho0 = density(gas) + h0 = enthalpy_mass(gas); u0 = internal_energy_mass(gas) + s0 = entropy_mass(gas); v0 = volume_mass(gas) + + perturb() = set_TPX!(gas, 600.0, 2 * one_atm, comp) # composition unchanged + + @testset "set_HP!" begin + perturb(); set_HP!(gas, h0, p0) + @test temperature(gas) ≈ T0 rtol=1e-6 + @test pressure(gas) ≈ p0 rtol=1e-10 + end + @testset "set_UV!" begin + perturb(); set_UV!(gas, u0, v0) + @test temperature(gas) ≈ T0 rtol=1e-6 + @test density(gas) ≈ rho0 rtol=1e-8 + end + @testset "set_SP!" begin + perturb(); set_SP!(gas, s0, p0) + @test temperature(gas) ≈ T0 rtol=1e-6 + @test pressure(gas) ≈ p0 rtol=1e-10 + end + @testset "set_SV!" begin + perturb(); set_SV!(gas, s0, v0) + @test temperature(gas) ≈ T0 rtol=1e-6 + @test density(gas) ≈ rho0 rtol=1e-8 + end + @testset "set_DP!" begin + perturb(); set_DP!(gas, rho0, p0) + @test density(gas) ≈ rho0 rtol=1e-10 + @test pressure(gas) ≈ p0 rtol=1e-10 + end + @testset "set_TD!" begin + perturb(); set_TD!(gas, T0, rho0) + @test temperature(gas) ≈ T0 rtol=1e-10 + @test density(gas) ≈ rho0 rtol=1e-10 + end + + close!(gas) +end + +@testset "equilibrate! modes" begin + @testset "TP is a fixed point" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 2000.0, one_atm, "CH4:1, O2:2, N2:7.52") + equilibrate!(gas, "TP") + @test temperature(gas) ≈ 2000.0 rtol=1e-8 + @test pressure(gas) ≈ one_atm rtol=1e-8 + X1 = mole_fractions(gas) + equilibrate!(gas, "TP") # idempotent + @test mole_fractions(gas) ≈ X1 atol=1e-8 + close!(gas) + end + @testset "HP conserves enthalpy & pressure (adiabatic combustion)" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 1500.0, one_atm, "H2:2, O2:1, N2:4") + h0 = enthalpy_mass(gas) + equilibrate!(gas, "HP") + @test enthalpy_mass(gas) ≈ h0 rtol=1e-6 + @test pressure(gas) ≈ one_atm rtol=1e-8 + @test temperature(gas) > 1500.0 # exothermic: burned gas is hotter + close!(gas) + end + @testset "UV conserves internal energy & volume" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 1500.0, one_atm, "H2:2, O2:1, N2:4") + u0 = internal_energy_mass(gas); v0 = volume_mass(gas) + equilibrate!(gas, "UV") + @test internal_energy_mass(gas) ≈ u0 rtol=1e-6 + @test volume_mass(gas) ≈ v0 rtol=1e-8 + close!(gas) + end +end diff --git a/interfaces/julia/test/test_transport.jl b/interfaces/julia/test/test_transport.jl new file mode 100644 index 0000000000..2af4d54a9b --- /dev/null +++ b/interfaces/julia/test/test_transport.jl @@ -0,0 +1,54 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +using Cantera +using Test +using LinearAlgebra: diag + + +@testset "Transport / mixture-averaged" begin + gas = Solution("gri30.yaml") + set_TPX!(gas, 1200.0, one_atm, "CH4:1, O2:2, N2:7.52") + nsp = n_species(gas) + + @test viscosity(gas) > 0 + @test thermal_conductivity(gas) > 0 + + D = mix_diff_coeffs(gas) + @test length(D) == nsp + @test all(D .> 0) + + Dbin = binary_diff_coeffs(gas) + @test size(Dbin) == (nsp, nsp) + @test all(isfinite, Dbin) + @test Dbin ≈ transpose(Dbin) # binary coefficients are symmetric + @test all(diag(Dbin) .> 0) + + # Diffusion coefficients grow with temperature at fixed pressure. + D1200 = mix_diff_coeffs(gas) + set_TPX!(gas, 1800.0, one_atm, "CH4:1, O2:2, N2:7.52") + @test all(mix_diff_coeffs(gas) .> D1200) + + close!(gas) +end + +@testset "Transport / multicomponent" begin + gas = Solution("gri30.yaml"; transport="multicomponent") + set_TPX!(gas, 1200.0, one_atm, "CH4:1, O2:2, N2:7.52") + nsp = n_species(gas) + + @test transport_model(gas) == "multicomponent" + @test viscosity(gas) > 0 + @test thermal_conductivity(gas) > 0 + + Dmulti = multi_diff_coeffs(gas) + @test size(Dmulti) == (nsp, nsp) + @test all(isfinite, Dmulti) # multicomponent coeffs may be negative + + DT = thermal_diff_coeffs(gas) + @test length(DT) == nsp + @test all(isfinite, DT) + @test any(!=(0.0), DT) + + close!(gas) +end diff --git a/interfaces/sourcegen/src/sourcegen/api.py b/interfaces/sourcegen/src/sourcegen/api.py index 00f7016244..cdd0fc3a5c 100644 --- a/interfaces/sourcegen/src/sourcegen/api.py +++ b/interfaces/sourcegen/src/sourcegen/api.py @@ -80,7 +80,7 @@ def create_argparser(): help="show additional logging output", ) parser.add_argument( - "--api", choices=["clib", "csharp", "yaml"], + "--api", choices=["clib", "csharp", "julia", "yaml"], help="language of generated Cantera API code", ) parser.add_argument( diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctkin.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctkin.yaml index 34e0cb1130..ecb99eca23 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctkin.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctkin.yaml @@ -44,6 +44,76 @@ recipes: - name: getCreationRates - name: getDestructionRates - name: getNetProductionRates +- name: getFwdRatesOfProgress_ddT +- name: getFwdRatesOfProgress_ddP +- name: getFwdRatesOfProgress_ddC +- name: getRevRatesOfProgress_ddT +- name: getRevRatesOfProgress_ddP +- name: getRevRatesOfProgress_ddC +- name: getNetRatesOfProgress_ddT +- name: getNetRatesOfProgress_ddP +- name: getNetRatesOfProgress_ddC +- name: getCreationRates_ddT +- name: getCreationRates_ddP +- name: getCreationRates_ddC +- name: getDestructionRates_ddT +- name: getDestructionRates_ddP +- name: getDestructionRates_ddC +- name: getNetProductionRates_ddT +- name: getNetProductionRates_ddP +- name: getNetProductionRates_ddC +- name: getNetProductionRates_ddX + brief: Dense composition Jacobian d(net production rate)/dX, column-major (n x n). + what: method + declaration: int32_t kin_getNetProductionRates_ddX(int32_t handle, int32_t bufLen, double* buf); + parameters: + handle: Handle to queried Kinetics object. + bufLen: Length of the reserved array; must be at least nTotalSpecies squared. + buf: Output buffer, filled column-major with the dense Jacobian. + code: |- + auto& obj = KineticsCabinet::at(handle); + auto mat = obj->netProductionRates_ddX(); + int32_t rows = static_cast(mat.rows()); + int32_t cols = static_cast(mat.cols()); + int32_t need = rows * cols; + if (bufLen < need) { + throw CanteraError("kin_getNetProductionRates_ddX", "bufLen is too small"); + } + for (int32_t idx = 0; idx < need; ++idx) { + buf[idx] = 0.0; + } + for (int k = 0; k < mat.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(mat, k); it; ++it) { + buf[it.col() * rows + it.row()] = it.value(); + } + } + return 0; +- name: getNetRatesOfProgress_ddX + brief: Dense reaction composition Jacobian d(net rate of progress)/dX, column-major. + what: method + declaration: int32_t kin_getNetRatesOfProgress_ddX(int32_t handle, int32_t bufLen, double* buf); + parameters: + handle: Handle to queried Kinetics object. + bufLen: Length of the reserved array; must be at least nReactions times nTotalSpecies. + buf: Output buffer, filled column-major with the dense Jacobian. + code: |- + auto& obj = KineticsCabinet::at(handle); + auto mat = obj->netRatesOfProgress_ddX(); + int32_t rows = static_cast(mat.rows()); + int32_t cols = static_cast(mat.cols()); + int32_t need = rows * cols; + if (bufLen < need) { + throw CanteraError("kin_getNetRatesOfProgress_ddX", "bufLen is too small"); + } + for (int32_t idx = 0; idx < need; ++idx) { + buf[idx] = 0.0; + } + for (int k = 0; k < mat.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(mat, k); it; ++it) { + buf[it.col() * rows + it.row()] = it.value(); + } + } + return 0; - name: multiplier - name: setMultiplier - name: isReversible diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctreactornet.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctreactornet.yaml index 131eff9185..b0b4a94332 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctreactornet.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctreactornet.yaml @@ -23,6 +23,10 @@ recipes: - name: atol - name: sensitivity wraps: sensitivity(const string&, size_t, int) +- name: neq # number of state variables / network components +- name: getState # current state vector, length neq() +- name: componentName # name of state-vector component i + uses: neq # service functions - name: del - name: cabinetSize diff --git a/interfaces/sourcegen/src/sourcegen/headers/ctthermo.yaml b/interfaces/sourcegen/src/sourcegen/headers/ctthermo.yaml index e4d961e646..3e47d64208 100644 --- a/interfaces/sourcegen/src/sourcegen/headers/ctthermo.yaml +++ b/interfaces/sourcegen/src/sourcegen/headers/ctthermo.yaml @@ -102,6 +102,11 @@ recipes: - name: getPartialMolarCp - name: getPartialMolarCv_TV - name: getPartialMolarVolumes +- name: getEnthalpy_RT # standard-state species enthalpies h^0_k/RT +- name: getEntropy_R # standard-state species entropies s^0_k/R +- name: getGibbs_RT # standard-state species Gibbs energies g^0_k/RT +- name: getIntEnergy_RT # standard-state species internal energies u^0_k/RT +- name: getCp_R # standard-state species heat capacities cp^0_k/R - name: setState_TPX # New in Cantera 3.2 wraps: setState_TPX(double, double, span) - name: setState_TPX_byName # New in Cantera 3.2 diff --git a/interfaces/sourcegen/src/sourcegen/julia/__init__.py b/interfaces/sourcegen/src/sourcegen/julia/__init__.py new file mode 100644 index 0000000000..3f07046979 --- /dev/null +++ b/interfaces/sourcegen/src/sourcegen/julia/__init__.py @@ -0,0 +1,4 @@ +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +from .generator import JuliaSourceGenerator diff --git a/interfaces/sourcegen/src/sourcegen/julia/config.yaml b/interfaces/sourcegen/src/sourcegen/julia/config.yaml new file mode 100644 index 0000000000..6b69a94e12 --- /dev/null +++ b/interfaces/sourcegen/src/sourcegen/julia/config.yaml @@ -0,0 +1,28 @@ +# Configuration for Julia code generation. + +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +# Ignore these files entirely: +ignore_files: [] + +# Ignore these specific functions: +ignore_funcs: {} + +# C type crosswalks; +# map C type names to the Julia types used in the `ccall` signature +c_type_crosswalk: + void: "Cvoid" + int32_t: "Int32" + const int32_t: "Int32" + int64_t: "Int64" + const int64_t: "Int64" + double: "Float64" + const double: "Float64" + const char*: "Cstring" + char*: "Ptr{UInt8}" + const double*: "Ptr{Float64}" + double*: "Ptr{Float64}" + const int32_t*: "Ptr{Int32}" + int32_t*: "Ptr{Int32}" + LogCallback: "Ptr{Cvoid}" diff --git a/interfaces/sourcegen/src/sourcegen/julia/generator.py b/interfaces/sourcegen/src/sourcegen/julia/generator.py new file mode 100644 index 0000000000..2162f46d5f --- /dev/null +++ b/interfaces/sourcegen/src/sourcegen/julia/generator.py @@ -0,0 +1,91 @@ +"""Generator for Julia source files.""" + +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +from pathlib import Path +import sys +import logging +from dataclasses import dataclass + +from jinja2 import Environment, BaseLoader + +from ..dataclasses import Func, HeaderFile +from ..generator import SourceGenerator + +from .._helpers import with_unpack_iter + + +_LOGGER = logging.getLogger() +_LOADER = Environment(loader=BaseLoader, trim_blocks=True, lstrip_blocks=True) + + +@dataclass(frozen=True) +@with_unpack_iter +class Config: + """Provides configuration info for the JuliaSourceGenerator class.""" + + c_type_crosswalk: dict[str, str] #: C type crosswalks + + +class JuliaSourceGenerator(SourceGenerator): + """The SourceGenerator for scaffolding Julia files for the Julia interface""" + + def __init__(self, out_dir: str, config: dict, templates: dict) -> None: + if not out_dir: + _LOGGER.critical("Non-empty string identifying output path required.") + sys.exit(1) + self._out_dir = Path(out_dir) / "src" / "generated" + + # use the typed config + self._config = Config(**config) + self._templates = templates + + def _crosswalk(self, c_type: str) -> str: + """Crosswalk of C/Julia types.""" + if c_type not in self._config.c_type_crosswalk: + msg = (f"Unable to crosswalk C type {c_type!r}: add an entry to the " + "'c_type_crosswalk' section of sourcegen/julia/config.yaml.") + _LOGGER.critical(msg) + sys.exit(1) + return self._config.c_type_crosswalk[c_type] + + def _scaffold_func(self, func: Func) -> str: + """Render the `ccall` wrapper for a single CLib function.""" + arg_names = [par.name for par in func.arglist] + arg_types = [self._crosswalk(par.p_type) for par in func.arglist] + + template = _LOADER.from_string(self._templates["julia-ccall-func"]) + return template.render( + name=func.name, + ret_type=self._crosswalk(func.ret_type), + arg_names=", ".join(arg_names), + arg_types="({},)".format(", ".join(arg_types)) if arg_types else "()", + call_args="".join(f", {name}" for name in arg_names)) + + def _write_file(self, file_name: str, template_name: str, **kwargs) -> None: + _LOGGER.info(f" writing {file_name!r}") + t_file = Path(__file__).parent / template_name + template = _LOADER.from_string(t_file.read_text(encoding="utf-8")) + contents = template.render(file_name=file_name, **kwargs) + + self._out_dir.joinpath(file_name).write_text(contents, encoding="utf-8") + + def generate_source(self, headers_files: list[HeaderFile]) -> None: + self._out_dir.mkdir(parents=True, exist_ok=True) + # delete any existing files, so that bindings for a header that is no + # longer generated cannot linger and shadow the current CLib + for f in self._out_dir.iterdir(): + f.unlink() + + generated_files = [] + for header_file in headers_files: + file_name = f"lib{header_file.path.stem}.jl" + self._write_file( + file_name, "template_bindings.jl.j2", + header_file=f"cantera_clib/{header_file.path.stem}.h", + functions=[self._scaffold_func(func) for func in header_file.funcs]) + generated_files.append(file_name) + + self._write_file("_manifest.jl", "template_manifest.jl.j2", + files=generated_files) diff --git a/interfaces/sourcegen/src/sourcegen/julia/template_bindings.jl.j2 b/interfaces/sourcegen/src/sourcegen/julia/template_bindings.jl.j2 new file mode 100644 index 0000000000..60199c7ec0 --- /dev/null +++ b/interfaces/sourcegen/src/sourcegen/julia/template_bindings.jl.j2 @@ -0,0 +1,16 @@ +# file: {{ file_name }} +# +# This file was generated by sourcegen. It will be re-generated by the +# Cantera build process. Do not manually edit. +# +# This file contains `ccall` wrappers related to the {{ header_file }} header +# for calling the native Cantera library. The wrappers are internal: the public +# Julia API lives in the hand-written `src/*.jl` files. +# +# Warning: This module is an experimental part of the Cantera API and +# may be changed or removed without notice. + +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +{{ functions | join("\n\n") }} diff --git a/interfaces/sourcegen/src/sourcegen/julia/template_manifest.jl.j2 b/interfaces/sourcegen/src/sourcegen/julia/template_manifest.jl.j2 new file mode 100644 index 0000000000..901f70009f --- /dev/null +++ b/interfaces/sourcegen/src/sourcegen/julia/template_manifest.jl.j2 @@ -0,0 +1,18 @@ +# file: {{ file_name }} +# +# This file was generated by sourcegen. It will be re-generated by the +# Cantera build process. Do not manually edit. +# +# Manifest of the generated binding files, included by `LibCantera.jl`. +# +# Warning: This module is an experimental part of the Cantera API and +# may be changed or removed without notice. + +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +const GENERATED_FILES = [ +{% for file in files %} + "{{ file }}", +{% endfor %} +] diff --git a/interfaces/sourcegen/src/sourcegen/julia/templates.yaml b/interfaces/sourcegen/src/sourcegen/julia/templates.yaml new file mode 100644 index 0000000000..54bbc8c9a1 --- /dev/null +++ b/interfaces/sourcegen/src/sourcegen/julia/templates.yaml @@ -0,0 +1,9 @@ +# Definitions used for Jinja template replacement. + +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +julia-ccall-func: |- + function {{ name }}({{ arg_names }}) + ccall((:{{ name }}, libcantera[]), {{ ret_type }}, {{ arg_types }}{{ call_args }}) + end