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.transformation.to_amrex as trafo

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[ \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} & - 1.0 n_{0} \left(n_{0} q_{1} + n_{1} q_{2}\right) - n_{0} \left(1.0 n_{0} q_{1} + 1.0 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} & - 1.0 n_{0} \left(n_{0} q_{1} + n_{1} q_{2}\right) - n_{0} \left(1.0 n_{0} q_{1} + 1.0 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} & - 1.0 n_{0} \left(n_{0} q_{1} + n_{1} q_{2}\right) - n_{0} \left(1.0 n_{0} q_{1} + 1.0 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} & - 1.0 n_{0} \left(n_{0} q_{1} + n_{1} q_{2}\right) - n_{0} \left(1.0 n_{0} q_{1} + 1.0 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-05 09:56:08.648 | WARNING  | library.core.misc.misc:__init__:146 - No 'clean_directory' attribute found in output Zstruct. Default: False

Check the output

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

class Model {
public:
    static constexpr int n_dof_q    = 3;
    static constexpr int n_dof_qaux = 2;
    static constexpr int dimension  = 2;



    AMREX_GPU_HOST_DEVICE
    AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,1>
    flux_x ( 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>{};
        amrex::Real t0 = (1.0 / amrex::Math::powi<1>(Q(0)));
        res(0,0) = Q(1);
        res(1,0) = (1.0/2.0)*9.81*amrex::Math::powi<2>(Q(0)) + amrex::Math::powi<2>(Q(1))*t0;
        res(2,0) = Q(1)*Q(2)*t0;
        return res;
    }
        


    AMREX_GPU_HOST_DEVICE
    AMREX_FORCE_INLINE
    static amrex::SmallMatrix<amrex::Real,3,1>
    flux_y ( 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>{};
        amrex::Real t0 = (1.0 / amrex::Math::powi<1>(Q(0)));
        res(0,0) = Q(2);
        res(1,0) = Q(1)*Q(2)*t0;
        res(2,0) = (1.0/2.0)*9.81*amrex::Math::powi<2>(Q(0)) + amrex::Math::powi<2>(Q(2))*t0;
        return res;
    }
        


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


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


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


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

    }
        


    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) 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) 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>
    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,3,1>
    residual ( 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,3,1>
    source_implicit ( 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,5,1>
    project_2d_to_3d ( 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,5,1>{};
        amrex::Real t0 = (1.0 / amrex::Math::powi<1>(Q(0)));
        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(2));
        return res;
    }

        


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

    }
        
};