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.transformation.to_openfoam as trafo
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 = os.getenv("ZOOMY_DIR")
mesh = petscMesh.Mesh.from_gmsh(
    os.path.join(main_dir, "meshes/channel_quad_2d/mesh.msh")
)
model.initialize_boundary_conditions(mesh)

Sympy Model

model.flux()

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

model.print_boundary_conditions()

\(\displaystyle \left[\begin{matrix}q_{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} & - 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}\\q_{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} & - 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}\\q_{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} & - 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]\)

Code transformation

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

2025-10-19 14:10:56.211 | WARNING  | library.core.misc.misc:__init__:146 - No 'clean_directory' attribute found in output Zstruct. Default: False

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[5], line 2
      1 settings = Settings(name="ShallowWater", output=Zstruct(directory="outputs/trafo", filename="swe.h5"))
----> 2 trafo.write_code(model, settings)

File ~/Git/Zoomy/library/python/transformation/to_openfoam.py:226, in write_code(model, settings, additional_writes)
    224 def write_code(model, settings, additional_writes=None):
    225     printer = FoamPrinter(model)
--> 226     code = printer.create_model(model, additional_writes=additional_writes)
    227     main_dir = os.getenv("ZOOMY_DIR")
    228     path = os.path.join(main_dir, settings.output.directory, ".foam_interface")

File ~/Git/Zoomy/library/python/transformation/to_openfoam.py:199, in FoamPrinter.create_model(self, model, additional_writes)
    197 funcs += self.create_function('flux', model.flux(), n_dof, n_dof_qaux)
    198 funcs += self.create_function('flux_jacobian', model.flux_jacobian(), n_dof, n_dof_qaux)
--> 199 funcs += self.create_function('dflux', model.dflux(), n_dof, n_dof_qaux)
    200 funcs += self.create_function('dflux_jacobian', model.dflux_jacobian(), n_dof, n_dof_qaux)
    201 funcs += self.create_function('nonconservative_matrix', model.nonconservative_matrix(), n_dof, n_dof_qaux)

File ~/Git/Zoomy/library/python/transformation/to_openfoam.py:123, in FoamPrinter.create_function(self, name, expr, n_dof_q, n_dof_qaux, target)
    121 if isinstance(expr, list):
    122     dim = len(expr)
--> 123     return [self.create_function(f"{name}_{d}", expr[i], n_dof_q, n_dof_qaux)
    124             for i, d in enumerate(['x', 'y', 'z'][:dim])]
    126 rows, cols = expr.shape
    127 body = self.convert_expression_body(expr, target)

File ~/Git/Zoomy/library/python/transformation/to_openfoam.py:126, in FoamPrinter.create_function(self, name, expr, n_dof_q, n_dof_qaux, target)
    122             dim = len(expr)
    123             return [self.create_function(f"{name}_{d}", expr[i], n_dof_q, n_dof_qaux)
    124                     for i, d in enumerate(['x', 'y', 'z'][:dim])]
--> 126         rows, cols = expr.shape
    127         body = self.convert_expression_body(expr, target)
    128         return f"""
    129 inline Foam::List<Foam::List<Foam::scalar>> {name}(
    130     const Foam::List<Foam::scalar>& Q,
   (...)    136 }}
    137         """

AttributeError: 'int' object has no attribute 'shape'

Check the output

main_dir = os.getenv("ZOOMY_DIR")
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"

namespace Model
{
constexpr int n_dof_q    = 3;
constexpr int n_dof_qaux = 2;
constexpr int dimension  = 2;
const Foam::List<Foam::word> map_boundary_tag_to_function_index{ "wall", "inflow", "outflow" };

inline Foam::List<Foam::List<Foam::scalar>> flux_x(
    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));
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 1));
        res[0][0] = Q[1];
        res[1][0] = (1.0/2.0)*9.81*Foam::pow(Q[0], 2) + Foam::pow(Q[1], 2)*t0;
        res[2][0] = Q[1]*Q[2]*t0;
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> flux_y(
    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));
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 1));
        res[0][0] = Q[2];
        res[1][0] = Q[1]*Q[2]*t0;
        res[2][0] = (1.0/2.0)*9.81*Foam::pow(Q[0], 2) + Foam::pow(Q[2], 2)*t0;
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> flux_jacobian_x(
    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>(3, 0.0));
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 2));
        Foam::scalar t1 = (1.0 / Foam::pow(Q[0], 1));
        Foam::scalar t2 = Q[1]*t1;
        res[0][0] = 0;
        res[0][1] = 1;
        res[0][2] = 0;
        res[1][0] = 9.81*Q[0] - Foam::pow(Q[1], 2)*t0;
        res[1][1] = 2*t2;
        res[1][2] = 0;
        res[2][0] = -Q[1]*Q[2]*t0;
        res[2][1] = Q[2]*t1;
        res[2][2] = t2;
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> flux_jacobian_y(
    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>(3, 0.0));
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 2));
        Foam::scalar t1 = (1.0 / Foam::pow(Q[0], 1));
        Foam::scalar t2 = Q[2]*t1;
        res[0][0] = 0;
        res[0][1] = 0;
        res[0][2] = 1;
        res[1][0] = -Q[1]*Q[2]*t0;
        res[1][1] = t2;
        res[1][2] = Q[1]*t1;
        res[2][0] = 9.81*Q[0] - Foam::pow(Q[2], 2)*t0;
        res[2][1] = 0;
        res[2][2] = 2*t2;
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> dflux_x(
    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] = Q[0];
        res[1][0] = Q[1];
        res[2][0] = Q[2];
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> dflux_y(
    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] = Q[0];
        res[1][0] = Q[1];
        res[2][0] = Q[2];
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> dflux_jacobian_x(
    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>(3, 0.0));
    res[0][0] = 1;
        res[0][1] = 0;
        res[0][2] = 0;
        res[1][0] = 0;
        res[1][1] = 1;
        res[1][2] = 0;
        res[2][0] = 0;
        res[2][1] = 0;
        res[2][2] = 1;
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> dflux_jacobian_y(
    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>(3, 0.0));
    res[0][0] = 1;
        res[0][1] = 0;
        res[0][2] = 0;
        res[1][0] = 0;
        res[1][1] = 1;
        res[1][2] = 0;
        res[2][0] = 0;
        res[2][1] = 0;
        res[2][2] = 1;
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> nonconservative_matrix_x(
    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>(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;
}
        

inline Foam::List<Foam::List<Foam::scalar>> nonconservative_matrix_y(
    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>(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;
}
        

inline Foam::List<Foam::List<Foam::scalar>> quasilinear_matrix_x(
    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>(3, 0.0));
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 2));
        Foam::scalar t1 = (1.0 / Foam::pow(Q[0], 1));
        Foam::scalar t2 = Q[1]*t1;
        res[0][0] = 0;
        res[0][1] = 1;
        res[0][2] = 0;
        res[1][0] = 9.81*Q[0] - Foam::pow(Q[1], 2)*t0;
        res[1][1] = 2*t2;
        res[1][2] = 0;
        res[2][0] = -Q[1]*Q[2]*t0;
        res[2][1] = Q[2]*t1;
        res[2][2] = t2;
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> quasilinear_matrix_y(
    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>(3, 0.0));
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 2));
        Foam::scalar t1 = (1.0 / Foam::pow(Q[0], 1));
        Foam::scalar t2 = Q[2]*t1;
        res[0][0] = 0;
        res[0][1] = 0;
        res[0][2] = 1;
        res[1][0] = -Q[1]*Q[2]*t0;
        res[1][1] = t2;
        res[1][2] = Q[1]*t1;
        res[2][0] = 9.81*Q[0] - Foam::pow(Q[2], 2)*t0;
        res[2][1] = 0;
        res[2][2] = 2*t2;
    return res;
}
        

inline Foam::List<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::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(1, 0.0));
    Foam::scalar t0 = n.x()*Q[1];
        Foam::scalar t1 = n.y()*Q[2];
        Foam::scalar t2 = (1.0 / Foam::pow(Q[0], 2));
        Foam::scalar t3 = Foam::pow(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][0] = (t0 + t1)/Q[0];
        res[1][0] = t2*(t3 + t4);
        res[2][0] = t2*(-t3 + t4);
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> left_eigenvectors(
    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>(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;
}
        

inline Foam::List<Foam::List<Foam::scalar>> right_eigenvectors(
    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>(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;
}
        

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;
}
        

inline Foam::List<Foam::List<Foam::scalar>> residual(
    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;
}
        

inline Foam::List<Foam::List<Foam::scalar>> source_implicit(
    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;
}
        

inline Foam::List<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::List<Foam::scalar>>(5, Foam::List<Foam::scalar>(1, 0.0));
    Foam::scalar t0 = (1.0 / Foam::pow(Q[0], 1));
        res[0][0] = Q[0];
        res[1][0] = Q[1]*t0;
        res[2][0] = Q[2]*t0;
        res[3][0] = 0;
        res[4][0] = 9810.0*Q[0]*(1 - X.z());
    return res;
}
        

inline Foam::List<Foam::List<Foam::scalar>> boundary_conditions(
    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::List<Foam::scalar>>(3, Foam::List<Foam::scalar>(3, 0.0));
    Foam::scalar t0 = n.y()*Q[2];
        Foam::scalar t1 = 1.0*n.x()*Q[1] + 1.0*t0;
        Foam::scalar t2 = 1.0*n.x()*Q[1] + 1.0*t0;
        Foam::scalar t3 = -n.x()*t1 - n.x()*t2 + 1.0*Q[1];
        Foam::scalar t4 = -n.y()*t1 - n.y()*t2 + 1.0*Q[2];
        res[0][0] = Q[0];
        res[0][1] = t3;
        res[0][2] = t4;
        res[1][0] = Q[0];
        res[1][1] = t3;
        res[1][2] = t4;
        res[2][0] = Q[0];
        res[2][1] = t3;
        res[2][2] = t4;
    return res;
}
        
} // namespace Model