Transformation to OpenFOAM Code

Author

Ingo Steldermann

Published

October 5, 2025

Transformation to OpenFOAM Code (Simple)

Imports

Load packages

import os
import numpy as np

from zoomy_jax.fvm.solver_jax import Settings
from zoomy_core.model.models.shallow_water import ShallowWaterEquations

import zoomy_core.model.initial_conditions as IC
import zoomy_core.model.boundary_conditions as BC
from zoomy_core.misc.misc import Zstruct
import zoomy_core.misc.misc as misc
from zoomy_core.transformation.to_openfoam import FoamModel
import zoomy_core.mesh.mesh as petscMesh

Model definition

bcs = BC.BoundaryConditions(
    [
        BC.Wall(tag="wall"),
        BC.Wall(tag="inflow"),
        BC.Wall(tag="outflow"),
    ]
)

def custom_ic(x):
    Q = np.zeros(3, dtype=float)
    Q[0] = np.where(x[0] < 5., 0.005, 0.001)
    return Q

ic = IC.UserFunction(custom_ic)

model = ShallowWaterEquations(
    dimension=2,
    boundary_conditions=bcs,
    initial_conditions=ic,
)

# the mesh is necessary to connect the boundary conditions to the mesh tags
main_dir = misc.get_main_directory()
mesh = petscMesh.Mesh.from_gmsh(
    os.path.join(main_dir, "meshes/channel_quad_2d/mesh.msh")
)

Sympy Model

model.flux()

\(\displaystyle \left[\begin{matrix}q_{1} & q_{2}\\0.5 ez g q_{0}^{2} + \frac{q_{1}^{2}}{q_{0}} & \frac{q_{1} q_{2}}{q_{0}}\\\frac{q_{1} q_{2}}{q_{0}} & 0.5 ez g q_{0}^{2} + \frac{q_{2}^{2}}{q_{0}}\end{matrix}\right]\)

model.print_boundary_conditions()

\(\displaystyle \begin{cases} \left[\begin{matrix}q_{0} & - 1.0 n_{0} \left(1.0 n_{0} q_{1} + 1.0 n_{1} q_{2}\right) - 1.0 n_{0} \left(n_{0} q_{1} + n_{1} q_{2}\right) + 1.0 q_{1} & - 1.0 n_{1} \left(1.0 n_{0} q_{1} + 1.0 n_{1} q_{2}\right) - 1.0 n_{1} \left(n_{0} q_{1} + n_{1} q_{2}\right) + 1.0 q_{2}\end{matrix}\right] & \text{for}\: bc_{idx} = 0 \vee bc_{idx} = 1 \vee bc_{idx} = 2 \end{cases}\)

Code transformation

settings = Settings(name="ShallowWater", output=Zstruct(directory="outputs/trafo", filename="swe.h5"))
FoamModel.write_code(model, settings)

2026-01-10 11:32:37.509 | WARNING  | zoomy_core.misc.misc:__init__:194 - No 'clean_directory' attribute found in output Zstruct. Default: False

2026-01-10 11:32:37.510 | WARNING  | zoomy_core.misc.misc:__init__:197 - No 'snapshots' attribute found in output Zstruct. Default: 2

'/home/ingo/Git/Zoomy/outputs/trafo/.foam_interface/Model.H'

Check the output

main_dir = misc.get_main_directory()
path = os.path.join(main_dir, os.path.join(settings.output.directory, '.foam_interface/Model.H'))
with open(path, "r") as f:
    print(f.read())

#pragma once
#include "List.H"
#include "vector.H"
#include "scalar.H"
#include <vector>
#include <string>


struct Model {

    // --- Constants ---
    static constexpr int n_dof_q    = 3;
    static constexpr int n_dof_qaux = 2;
    static constexpr int dimension  = 2;
    static constexpr int n_boundary_tags = 3;

    // --- Helpers ---
    static const std::vector<std::string> get_boundary_tags() { return { "inflow", "outflow", "wall" }; }

    // --- Kernels ---

    static inline Foam::List<Foam::List<Foam::scalar>> flux(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux)
    {
        auto res = Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(2, 0.0));
    Foam::scalar t0 = 0.5*1.0*9.81*Foam::pow(Q[0], 2);
    Foam::scalar t1 = (1.0 / Foam::pow(Q[0], 1));
    Foam::scalar t2 = Q[1]*Q[2]*t1;
    res[0][0] = Q[1];
    res[0][1] = Q[2];
    res[1][0] = Foam::pow(Q[1], 2)*t1 + t0;
    res[1][1] = t2;
    res[2][0] = t2;
    res[2][1] = Foam::pow(Q[2], 2)*t1 + t0;
        return res;
    }


    static inline Foam::List<Foam::List<Foam::scalar>> dflux(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux)
    {
        auto res = Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(2, 0.0));
    res[0][0] = 0;
    res[0][1] = 0;
    res[1][0] = 0;
    res[1][1] = 0;
    res[2][0] = 0;
    res[2][1] = 0;
        return res;
    }


    static inline Foam::List<Foam::List<Foam::List<Foam::scalar>>> nonconservative_matrix(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux)
    {
        auto res = Foam::List<Foam::List<Foam::List<Foam::scalar>>>(3, Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(2, 0.0)));
    res[0][0][0] = 0;
    res[0][0][1] = 0;
    res[0][1][0] = 0;
    res[0][1][1] = 0;
    res[0][2][0] = 0;
    res[0][2][1] = 0;
    res[1][0][0] = 0;
    res[1][0][1] = 0;
    res[1][1][0] = 0;
    res[1][1][1] = 0;
    res[1][2][0] = 0;
    res[1][2][1] = 0;
    res[2][0][0] = 0;
    res[2][0][1] = 0;
    res[2][1][0] = 0;
    res[2][1][1] = 0;
    res[2][2][0] = 0;
    res[2][2][1] = 0;
        return res;
    }


    static inline Foam::List<Foam::List<Foam::List<Foam::scalar>>> quasilinear_matrix(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux)
    {
        auto res = Foam::List<Foam::List<Foam::List<Foam::scalar>>>(3, Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(2, 0.0)));
    Foam::scalar t0 = 1.0*1.0*9.81*Q[0];
    Foam::scalar t1 = (1.0 / Foam::pow(Q[0], 2));
    Foam::scalar t2 = -Q[1]*Q[2]*t1;
    Foam::scalar t3 = (1.0 / Foam::pow(Q[0], 1));
    Foam::scalar t4 = Q[1]*t3;
    Foam::scalar t5 = Q[2]*t3;
    res[0][0][0] = 0;
    res[0][0][1] = 0;
    res[0][1][0] = 1;
    res[0][1][1] = 0;
    res[0][2][0] = 0;
    res[0][2][1] = 1;
    res[1][0][0] = -Foam::pow(Q[1], 2)*t1 + t0;
    res[1][0][1] = t2;
    res[1][1][0] = 2*t4;
    res[1][1][1] = t5;
    res[1][2][0] = 0;
    res[1][2][1] = t4;
    res[2][0][0] = t2;
    res[2][0][1] = -Foam::pow(Q[2], 2)*t1 + t0;
    res[2][1][0] = t5;
    res[2][1][1] = 0;
    res[2][2][0] = t4;
    res[2][2][1] = 2*t5;
        return res;
    }


    static inline Foam::List<Foam::List<Foam::scalar>> source(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux)
    {
        auto res = Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(1, 0.0));
    res[0][0] = 0;
    res[1][0] = 0;
    res[2][0] = 0;
        return res;
    }


    static inline Foam::List<Foam::List<Foam::List<Foam::scalar>>> source_jacobian_wrt_variables(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux)
    {
        auto res = Foam::List<Foam::List<Foam::List<Foam::scalar>>>(3, Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(1, 0.0)));
    res[0][0][0] = 0;
    res[0][1][0] = 0;
    res[0][2][0] = 0;
    res[1][0][0] = 0;
    res[1][1][0] = 0;
    res[1][2][0] = 0;
    res[2][0][0] = 0;
    res[2][1][0] = 0;
    res[2][2][0] = 0;
        return res;
    }


    static inline Foam::List<Foam::List<Foam::List<Foam::scalar>>> source_jacobian_wrt_aux_variables(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux)
    {
        auto res = Foam::List<Foam::List<Foam::List<Foam::scalar>>>(2, Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(1, 0.0)));
    res[0][0][0] = 0;
    res[0][1][0] = 0;
    res[0][2][0] = 0;
    res[1][0][0] = 0;
    res[1][1][0] = 0;
    res[1][2][0] = 0;
        return res;
    }


    static inline Foam::List<Foam::scalar> eigenvalues(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux,
        const Foam::vector& n)
    {
        auto res = Foam::List<Foam::scalar>(3, 0.0);
    Foam::scalar t0 = 1.0*n.x()*Q[1];
    Foam::scalar t1 = 1.0*n.y()*Q[2];
    Foam::scalar t2 = (1.0 / Foam::pow(Q[0], 2));
    Foam::scalar t3 = 1.0*Foam::pow(1.0*9.81*Foam::pow(Q[0], 5), 1.0/2.0)*Foam::pow(Foam::pow(n.x(), 2) + Foam::pow(n.y(), 2), 1.0/2.0);
    Foam::scalar t4 = Q[0]*t0 + Q[0]*t1;
    res[0] = (t0 + t1)/Q[0];
    res[1] = t2*(t3 + t4);
    res[2] = t2*(-t3 + t4);
        return res;
    }


    static inline Foam::List<Foam::List<Foam::scalar>> left_eigenvectors(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux,
        const Foam::vector& n)
    {
        auto res = Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(3, 0.0));
    res[0][0] = 0;
    res[0][1] = 0;
    res[0][2] = 0;
    res[1][0] = 0;
    res[1][1] = 0;
    res[1][2] = 0;
    res[2][0] = 0;
    res[2][1] = 0;
    res[2][2] = 0;
        return res;
    }


    static inline Foam::List<Foam::List<Foam::scalar>> right_eigenvectors(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux,
        const Foam::vector& n)
    {
        auto res = Foam::List<Foam::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(3, 0.0));
    res[0][0] = 0;
    res[0][1] = 0;
    res[0][2] = 0;
    res[1][0] = 0;
    res[1][1] = 0;
    res[1][2] = 0;
    res[2][0] = 0;
    res[2][1] = 0;
    res[2][2] = 0;
        return res;
    }


    static inline Foam::List<Foam::scalar> residual(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux,
        const Foam::vector& X,
        const Foam::scalar& time,
        const Foam::scalar& dX)
    {
        auto res = Foam::List<Foam::scalar>(3, 0.0);
    res[0] = 0;
    res[1] = 0;
    res[2] = 0;
        return res;
    }


    static inline Foam::List<Foam::scalar> interpolate(
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux,
        const Foam::vector& X)
    {
        auto res = Foam::List<Foam::scalar>(6, 0.0);
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 1));
    res[0] = 0.0;
    res[1] = Q[0];
    res[2] = Q[1]*t0;
    res[3] = Q[2]*t0;
    res[4] = 0.0;
    res[5] = 9810.0*Q[0]*(1 - X.z());
        return res;
    }


    static inline Foam::List<Foam::scalar> boundary_conditions(
        const int bc_idx,
        const Foam::List<Foam::scalar>& Q,
        const Foam::List<Foam::scalar>& Qaux,
        const Foam::vector& n,
        const Foam::vector& X,
        const Foam::scalar& time,
        const Foam::scalar& dX)
    {
        auto res = Foam::List<Foam::scalar>(3, 0.0);
    if (bc_idx == 0 || bc_idx == 1 || bc_idx == 2) {
        Foam::scalar t0 = n.y()*Q[2];
        Foam::scalar t1 = n.x()*Q[1] + t0;
        Foam::scalar t2 = 1.0*n.x();
        Foam::scalar t3 = 1.0*n.x()*Q[1] + 1.0*t0;
        Foam::scalar t4 = 1.0*n.y();
        res[0] = Q[0];
        res[1] = 1.0*Q[1] - t1*t2 - t2*t3;
        res[2] = 1.0*Q[2] - t1*t4 - t3*t4;
    }
        return res;
    }

};