Triangular Mesh Guide

This section provides a comprehensive guide to using unstructured triangular meshes in FLUXOS.

Overview

FLUXOS supports unstructured triangular meshes as an alternative to the default regular Cartesian grid. The triangular mesh solver uses an edge-based finite volume formulation with a rotated Roe/HLL flux solver, providing second-order accuracy on arbitrary mesh geometries.

When to use triangular meshes:

  • Irregular domain boundaries (coastlines, river banks, infrastructure)

  • Local mesh refinement needed (e.g., finer resolution near channels or structures)

  • Complex channel geometries that are poorly represented by Cartesian grids

  • Domains where uniform resolution is wasteful (large flat areas with small features)

When to use regular meshes:

  • Rectangular domains with uniform features

  • Maximum computational efficiency on simple geometries

  • Compatibility with existing Cartesian-based workflows

Build Requirements

Triangular mesh support requires the USE_TRIMESH CMake option:

cmake -DMODE_release=ON ..
make -j$(nproc)

For GPU acceleration of the triangular mesh solver:

cmake -DMODE_release=ON -DUSE_CUDA=ON ..
make -j$(nproc)

Mesh Formats

FLUXOS supports two standard triangular mesh formats:

Gmsh .msh v2.2

Gmsh is a popular open-source mesh generator. FLUXOS reads the .msh v2.2 ASCII format.

Key sections parsed:

  • $Nodes: Vertex coordinates (x, y, z)

  • $Elements: Triangles (type 2) and boundary lines (type 1)

  • Physical group tags map to boundary condition types

Creating a mesh in Gmsh:

  1. Define the domain geometry in Gmsh (points, lines, surfaces)

  2. Assign physical groups to boundary lines (e.g., group 1 = wall, group 2 = outflow)

  3. Generate the mesh: Mesh > 2D

  4. Save as .msh format version 2.2

// Example Gmsh .geo file
Point(1) = {0, 0, 0, 1.0};
Point(2) = {100, 0, 0, 1.0};
Point(3) = {100, 50, 0, 1.0};
Point(4) = {0, 50, 0, 1.0};

Line(1) = {1, 2};    // bottom wall
Line(2) = {2, 3};    // outflow
Line(3) = {3, 4};    // top wall
Line(4) = {4, 1};    // inflow wall

Line Loop(1) = {1, 2, 3, 4};
Plane Surface(1) = {1};

Physical Line("wall") = {1, 3, 4};   // tag 1
Physical Line("outflow") = {2};       // tag 2
Physical Surface("domain") = {1};

Triangle .node/.ele

The Triangle mesh generator produces three files:

  • .node: Vertex coordinates

  • .ele: Triangle connectivity

  • .edge: Edge information (optional)

Specify the base filename (without extension) in the configuration.

JSON Configuration

To use a triangular mesh, add these settings to your master JSON file:

{
    "MESH_TYPE": "triangular",
    "MESH_FILE": "path/to/domain.msh",
    "MESH_FORMAT": "gmsh",
    "BOUNDARY_CONDITIONS": {
        "1": {"type": "wall"},
        "2": {"type": "outflow"}
    }
}

Configuration options:

Key

Values

Description

MESH_TYPE

"regular" or "triangular"

Selects the mesh type (default: "regular")

MESH_FILE

file path

Path to the mesh file (relative to the JSON config directory)

MESH_FORMAT

"gmsh" or "triangle"

Mesh file format

BOUNDARY_CONDITIONS

JSON object

Maps physical group tags to boundary types

Boundary condition types:

  • "wall" (default): Reflective boundary - normal velocity is reflected, tangential velocity preserved

  • "outflow": Transmissive boundary - zero-gradient free outflow

Note

The DEM_FILE entry is still required in the JSON config even with triangular meshes. However, if the .msh file contains non-zero vertex z-coordinates (as generated by the Python configuration template), the solver uses vertex elevations directly – computing cell bed elevation as the mean of its three vertex z values. If vertex z values are all zero (e.g., manually created mesh), the solver falls back to nearest-neighbor interpolation from the DEM grid.

Adaptive Mesh Generation from GeoTIFF

The Python configuration template (supporting_scripts/1_Model_Config/model_config_template.py) can automatically generate adaptive triangular meshes from GeoTIFF DEM data. The mesh generator creates finer elements in steep terrain and coarser elements in flat areas, optimizing computational resources.

Edit the mesh block of the template’s _config dict:

mesh_type            = "triangular",
trimesh_min_size     = 5.0,
trimesh_max_size     = 50.0,
trimesh_slope_factor = 1.0,
mesh_output_msh      = "domain.msh",

then run python model_config_template.py to produce bin/domain.msh (plus the matching .asc and modset.json).

Parameters:

  • trimesh_min_size: Minimum triangle edge length in meters (used at steepest slopes)

  • trimesh_max_size: Maximum triangle edge length in meters (used in flat areas)

  • trimesh_slope_factor: Aggressiveness of slope-based refinement (0 = uniform, higher values = more refinement in steep terrain)

How it works:

  1. The tool reads the GeoTIFF and computes slope: \(|\nabla z| = \sqrt{(\partial z/\partial x)^2 + (\partial z/\partial y)^2}\)

  2. Slope is mapped to element size: size = max_size - slope_norm * (max_size - min_size) * slope_factor

  3. The domain boundary is extracted from valid (non-NODATA) DEM cells

  4. Gmsh generates the mesh using the slope-based size field as a background mesh

  5. DEM elevations are interpolated (bilinear) to each mesh vertex and stored as the z-coordinate

  6. The output .msh v2.2 file contains vertex (x, y, z) where z = DEM elevation

The C++ solver automatically detects non-zero vertex z values and uses them for cell bed elevation (mean of 3 vertex z values per cell), bypassing the DEM grid interpolation. A companion .asc file is also generated for the DEM_FILE config entry.

See Digital Elevation Model (DEM) for full documentation of the preprocessing tool.

Numerical Methods

The triangular mesh solver employs the following algorithms:

6-Phase Hydrodynamics Driver:

Each timestep performs:

  1. Wet/dry classification: Identify dry cells and update status

  2. Gradient computation: Least-squares gradient using neighbor cell values with precomputed 2x2 inverse matrix

  3. Slope limiting: Barth-Jespersen limiter applied per cell to ensure monotonicity

  4. Edge flux computation: MUSCL-reconstructed left/right states at each edge, rotated Roe/HLL solver in edge-normal frame

  5. Flux accumulation: Sum edge fluxes into cell residuals

  6. State update: Advance state variables, apply friction, enforce wet/dry constraints

CFL Condition:

The time step uses the cell inradius (inscribed circle radius) for stability:

\[\Delta t = CFL \cdot \min_i \frac{r_i}{|\mathbf{u}_i| + \sqrt{g h_i}}, \quad r_i = \frac{2 A_i}{P_i}\]

This naturally accounts for varying cell sizes across the mesh.

Edge-Based ADE Solver:

Solute transport uses:

  • Upwind advection based on hydrodynamic mass flux direction at each edge

  • Fick’s law dispersion with effective coefficient

  • Mass conservation limiting (outgoing flux bounded by available mass)

  • Concentration bounding by neighbor min/max

Output Format

Triangular mesh results are output as VTK XML Unstructured Grid files (.vtu):

  • One .vtu file per output timestep

  • A .pvd time series collection file for ParaView animation

  • Cell-centered data arrays: h, z, zb, ux, uy, velocity_magnitude, qx, qy, us, conc_SW, twetimetracer, cell_area

See the Output files page for detailed output format documentation and the Supporting Scripts page for visualization tools, including the KML exporter (fluxos_viewer.py) that can export triangular mesh results as animated KMZ files for Google Earth.

GPU Acceleration

When built with both USE_TRIMESH and USE_CUDA, the triangular mesh solver uses 7 specialized CUDA kernels:

Kernel

Parallelism

Description

kernel_tri_wetdry

1 thread/cell

Wet/dry classification and status update

kernel_tri_gradient

1 thread/cell

Least-squares gradient computation

kernel_tri_limiter

1 thread/cell

Barth-Jespersen slope limiter

kernel_tri_edge_flux

1 thread/edge

Rotated Roe/HLL flux (main compute kernel)

kernel_tri_accumulate

1 thread/edge

Flux accumulation with atomicAdd

kernel_tri_update

1 thread/cell

State update + wet/dry clamping

kernel_tri_courant

1 thread/cell

CFL time step reduction

MPI Parallelization

For distributed computing with triangular meshes:

  • Partitioning: METIS graph partitioning (if available) or naive block partitioning

  • Halo exchange: Cells adjacent to partition boundaries are exchanged between MPI ranks

  • Communication: MPI_Isend/MPI_Irecv per neighbor rank

Build with MPI:

cmake -DMODE_release=ON -DUSE_MPI=ON ..
make -j$(nproc)
mpirun -np 4 ./bin/fluxos_mpi ./input/modset_triangular.json

Source Code Organization

All triangular mesh code is contained in src/fluxos/trimesh/:

File

Purpose

TriMesh.h/.cpp

Mesh topology, geometry, edge building, LSQ matrices

TriSolution.h/.cpp

Flat-array solution storage and allocation

tri_mesh_io.h/.cpp

Gmsh and Triangle format readers

tri_initiate.h/.cpp

Initialization, DEM interpolation, boundary conditions

tri_gradient.h/.cpp

LSQ gradient + Barth-Jespersen limiter

tri_solver_wet.h/.cpp

Rotated Roe/HLL flux with MUSCL reconstruction

tri_solver_dry.h/.cpp

Ritter dry-front solution

tri_hydrodynamics.h/.cpp

6-phase driver + CFL condition

tri_ade_solver.h/.cpp

Edge-based advection-dispersion

tri_forcing.h/.cpp

Meteo and inflow forcing

tri_write_results.h/.cpp

VTK .vtu output + .pvd time series

tri_mpi_domain.h/.cpp

MPI decomposition + halo exchange

tri_cuda_memory.h/.cu

GPU memory management

tri_cuda_kernels.h/.cu

CUDA kernels (7 kernels)