Transformation to C Code

Author

Ingo Steldermann

Published

August 18, 2025

Transformation to C 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_amrex import AmrexModel

Model definition

bcs = BC.BoundaryConditions(
    [
        BC.Wall(tag="top"),
        BC.Wall(tag="bottom"),
        BC.Wall(tag="left"),
        BC.Wall(tag="right"),
    ]
)

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,
)

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 \vee bc_{idx} = 3 \end{cases}\)

Code transformation

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

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

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

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

Check the output

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

#pragma once
#include <AMReX_Array4.H>
#include <AMReX_Vector.H>
#include <AMReX_SmallMatrix.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 = 4;

    // --- Helpers ---
    static const std::vector<std::string> get_boundary_tags() { return { "bottom", "left", "right", "top" }; }

    // --- Kernels ---

    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,2> flux(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,2>{};
    amrex::Real t0 = 0.5*1.0*9.81*amrex::Math::powi<2>(Q(0));
    amrex::Real t1 = (1.0 / amrex::Math::powi<1>(Q(0)));
    amrex::Real t2 = Q(1)*Q(2)*t1;
    res(0,0) = Q(1);
    res(0,1) = Q(2);
    res(1,0) = amrex::Math::powi<2>(Q(1))*t1 + t0;
    res(1,1) = t2;
    res(2,0) = t2;
    res(2,1) = amrex::Math::powi<2>(Q(2))*t1 + t0;
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,2> dflux(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,2>{};
    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;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,18,1> nonconservative_matrix(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,18,1>{};
    res(0) = 0;
    res(1) = 0;
    res(2) = 0;
    res(3) = 0;
    res(4) = 0;
    res(5) = 0;
    res(6) = 0;
    res(7) = 0;
    res(8) = 0;
    res(9) = 0;
    res(10) = 0;
    res(11) = 0;
    res(12) = 0;
    res(13) = 0;
    res(14) = 0;
    res(15) = 0;
    res(16) = 0;
    res(17) = 0;
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,18,1> quasilinear_matrix(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,18,1>{};
    amrex::Real t0 = 1.0*1.0*9.81*Q(0);
    amrex::Real t1 = (1.0 / amrex::Math::powi<2>(Q(0)));
    amrex::Real t2 = -Q(1)*Q(2)*t1;
    amrex::Real t3 = (1.0 / amrex::Math::powi<1>(Q(0)));
    amrex::Real t4 = Q(1)*t3;
    amrex::Real t5 = Q(2)*t3;
    res(0) = 0;
    res(1) = 0;
    res(2) = 1;
    res(3) = 0;
    res(4) = 0;
    res(5) = 1;
    res(6) = -amrex::Math::powi<2>(Q(1))*t1 + t0;
    res(7) = t2;
    res(8) = 2*t4;
    res(9) = t5;
    res(10) = 0;
    res(11) = t4;
    res(12) = t2;
    res(13) = -amrex::Math::powi<2>(Q(2))*t1 + t0;
    res(14) = t5;
    res(15) = 0;
    res(16) = t4;
    res(17) = 2*t5;
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,1> source(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,1>{};
    res(0,0) = 0;
    res(1,0) = 0;
    res(2,0) = 0;
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,9,1> source_jacobian_wrt_variables(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,9,1>{};
    res(0) = 0;
    res(1) = 0;
    res(2) = 0;
    res(3) = 0;
    res(4) = 0;
    res(5) = 0;
    res(6) = 0;
    res(7) = 0;
    res(8) = 0;
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,6,1> source_jacobian_wrt_aux_variables(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,6,1>{};
    res(0) = 0;
    res(1) = 0;
    res(2) = 0;
    res(3) = 0;
    res(4) = 0;
    res(5) = 0;
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,1> eigenvalues(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux,
        amrex::SmallMatrix<amrex::Real,2,1> const& n) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,1>{};
    amrex::Real t0 = 1.0*n(0)*Q(1);
    amrex::Real t1 = 1.0*n(1)*Q(2);
    amrex::Real t2 = (1.0 / amrex::Math::powi<2>(Q(0)));
    amrex::Real t3 = 1.0*amrex::Math::pow(1.0*9.81*amrex::Math::powi<5>(Q(0)), 1.0/2.0)*amrex::Math::pow(amrex::Math::powi<2>(n(0)) + amrex::Math::powi<2>(n(1)), 1.0/2.0);
    amrex::Real 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;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,3> left_eigenvectors(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux,
        amrex::SmallMatrix<amrex::Real,2,1> const& n) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,3>{};
    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;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,3> right_eigenvectors(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux,
        amrex::SmallMatrix<amrex::Real,2,1> const& n) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,3>{};
    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;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,1> residual(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux,
        amrex::SmallMatrix<amrex::Real,3,1> const& X,
        amrex::Real const& time,
        amrex::Real const& dX) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,1>{};
    res(0) = 0;
    res(1) = 0;
    res(2) = 0;
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,6,1> interpolate(
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux,
        amrex::SmallMatrix<amrex::Real,3,1> const& X) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,6,1>{};
    amrex::Real t0 = (1.0 / amrex::Math::powi<1>(Q(0)));
    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(2));
        return res;
    }


    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,1> boundary_conditions(
        const int bc_idx,
        amrex::SmallMatrix<amrex::Real,3,1> const& Q,
        amrex::SmallMatrix<amrex::Real,2,1> const& Qaux,
        amrex::SmallMatrix<amrex::Real,2,1> const& n,
        amrex::SmallMatrix<amrex::Real,3,1> const& X,
        amrex::Real const& time,
        amrex::Real const& dX) noexcept
    {
        auto res = amrex::SmallMatrix<amrex::Real,3,1>{};
    if (bc_idx == 0 || bc_idx == 1 || bc_idx == 2 || bc_idx == 3) {
        amrex::Real t0 = n(1)*Q(2);
        amrex::Real t1 = n(0)*Q(1) + t0;
        amrex::Real t2 = 1.0*n(0);
        amrex::Real t3 = 1.0*n(0)*Q(1) + 1.0*t0;
        amrex::Real t4 = 1.0*n(1);
        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;
    }

};