Shallow Water Tutorials 1D with Numpy

Author

Ingo Steldermann

Published

July 10, 2025

Shallow Water Tutorial 1D with Numpy (Simple)

Reference

The following the model is described in the paper:

 @article{Delestre_2013, 
 title={SWASHES: a compilation of Shallow Water Analytic Solutions for Hydraulic and Environmental Studies}, 
 volume={72}, 
 ISSN={0271-2091, 1097-0363}, DOI={10.1002/fld.3741}, 
 note={arXiv:1110.0288 [physics]}, 
 number={3}, 
 journal={International Journal for Numerical Methods in Fluids}, 
 author={Delestre, Olivier and Lucas, Carine and Ksinant, Pierre-Antoine and Darboux, Frédéric and Laguerre, Christian and Vo, Thi Ngoc Tuoi and James, Francois and Cordier, Stephane}, 
 year={2013}, 
 month=may, 
 pages={269–300} 
}

Imports

Load packages

import os
import numpy as np
import pytest
from sympy import Matrix
import pytest

from zoomy_core.fvm.solver_numpy import HyperbolicSolver, Settings
import zoomy_core.fvm.timestepping as timestepping
import zoomy_core.fvm.nonconservative_flux as nc_flux
from zoomy_core.model.basemodel import Model
from attr import field, define
import zoomy_core.model.initial_conditions as IC
import zoomy_core.model.boundary_conditions as BC
import zoomy_core.misc.io as io
from zoomy_core.mesh.mesh import compute_derivatives
from zoomy_core.misc.misc import Zstruct, ZArray
from zoomy_tests.swashes import plots_paper


import zoomy_core.mesh.mesh as petscMesh
import zoomy_core.postprocessing.postprocessing as postprocessing

@define(frozen=True, slots=True, kw_only=True)
class SWE(Model):
    dimension: int = 1
    variables: Zstruct = field(init=False, default=dimension + 1)
    aux_variables: Zstruct = field(default=1)
    _default_parameters: dict = field(
        init=False,
        factory=lambda: {"g": 9.81, "ex": 0.0, "ey": 0.0, "ez": 1.0}
        )


    def project_2d_to_3d(self):
        out = ZArray.zeros(6)
        dim = self.dimension
        z = self.position[2]
        b = 0
        h = self.variables[0]
        U = [hu / h for hu in self.variables[1:1+dim]]
        rho_w = 1000.
        g = 9.81
        out[0] = b
        out[1] = h
        out[2] = U[0]
        out[3] = 0 if dim == 1 else U[1]
        out[4] = 0
        out[5] = rho_w * g * h * (1-z)
        return out

    def flux(self):
        dim = self.dimension
        h = self.variables[0]
        U = Matrix([hu / h for hu in self.variables[1:1+dim]])
        g = self.parameters.g
        I = Matrix.eye(dim)
        F = Matrix.zeros(self.variables.length(), dim)
        F[0, :] = (h * U)
        F[1:, :] = h * U * U.T + g/2 * h**2 * I
        
        return F


bcs = BC.BoundaryConditions(
    [
        BC.Extrapolation(tag="left"),
        BC.Extrapolation(tag="right"),
    ]
)

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

ic = IC.UserFunction(custom_ic)

model = SWE(
    boundary_conditions=bcs,
    initial_conditions=ic,
)

main_dir = os.getenv("ZOOMY_DIR")
# mesh = petscMesh.Mesh.from_gmsh(
#     os.path.join(main_dir, "meshes/channel_quad_2d/mesh.msh")
# )

mesh = petscMesh.Mesh.create_1d(domain=(0.0, 10.0), n_inner_cells=500)

class SWESolver(HyperbolicSolver):
    def update_qaux(self, Q, Qaux, Qold, Qauxold, mesh, model, parameters, time, dt):
        dudx = compute_derivatives(Q[1]/Q[0], mesh, derivatives_multi_index=[[0, 0]])[:,0]
        Qaux[0]=dudx
        return Qaux
    
settings = Settings(name="ShallowWater", output=Zstruct(directory="outputs/shallow_water_1d", filename="swe", clean_directory=True))

solver = SWESolver(time_end=6.0, settings=settings, compute_dt=timestepping.adaptive(CFL=0.95))
Qnew, Qaux = solver.solve(mesh, model)

2025-11-06 12:12:41.400 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 10, time: 0.718260, dt: 0.067958, next write at time: 0.666667

2025-11-06 12:12:41.440 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 20, time: 1.392979, dt: 0.067263, next write at time: 1.333333

2025-11-06 12:12:41.470 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 30, time: 2.064282, dt: 0.067054, next write at time: 2.000000

2025-11-06 12:12:41.500 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 40, time: 2.734237, dt: 0.066956, next write at time: 3.333333

2025-11-06 12:12:41.529 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 50, time: 3.403455, dt: 0.066898, next write at time: 4.000000

2025-11-06 12:12:41.558 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 60, time: 4.072203, dt: 0.066858, next write at time: 4.666667

2025-11-06 12:12:41.588 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 70, time: 4.740621, dt: 0.066830, next write at time: 5.333333

2025-11-06 12:12:41.619 | INFO     | zoomy_core.fvm.solver_numpy:run:285 - iteration: 80, time: 5.408794, dt: 0.066808, next write at time: 6.000000

2025-11-06 12:12:41.645 | INFO     | zoomy_core.fvm.solver_numpy:solve:295 - Finished simulation with in 0.291 seconds

io.generate_vtk(os.path.join(settings.output.directory, f"{settings.output.filename}.h5"))
postprocessing.vtk_project_2d_to_3d(model, settings, filename='out_3d')

2025-11-06 12:12:43.452 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 0/10

2025-11-06 12:12:43.455 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 1/10

2025-11-06 12:12:43.456 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 2/10

2025-11-06 12:12:43.458 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 3/10

2025-11-06 12:12:43.460 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 4/10

2025-11-06 12:12:43.462 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 5/10

2025-11-06 12:12:43.465 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 6/10

2025-11-06 12:12:43.468 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 7/10

2025-11-06 12:12:43.470 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 8/10

2025-11-06 12:12:43.472 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:68 - Converted snapshot 9/10

2025-11-06 12:12:43.498 | INFO     | zoomy_core.postprocessing.postprocessing:vtk_project_2d_to_3d:71 - Output is written to: /home/ingo/Git/Zoomy/outputs/shallow_water_1d/out_3d.h5/out_3d.*.vtk

fig = plots_paper.plot_swe(os.path.join(settings.output.directory, settings.output.filename + ".h5"))

{'length': 10.0, 'width': None, 'dx': 0.2, 'dy': None, 'ncellx': 50.0, 'ncelly': None}

@pytest.mark.nbworking
def test_working():
    assert True