samurai gallery

2D Riemann problem (config 3)

Four constant states meet at a point and interact through shocks, contacts and a rolled-up jet.

ρ 0.138 to 2.07
5 mesh levels 29,200 cells 41 frames

The unit square is split into four quadrants, each holding a different constant state of a compressible gas. When released, the four states interact through a rich pattern of shocks, contact discontinuities and a rolled-up central jet - one of the classic two-dimensional Riemann configurations (Schulz-Rinne). It is a favorite stress test for high-resolution and adaptive schemes.

This case is powered by the samurai-euler solver engine; the code shown here is the scenario (the four quadrant states and boundary conditions).

Equations

The compressible Euler equations for an ideal gas (γ=1.4\gamma = 1.4), with the initial state piecewise constant on the four quadrants of [0,1]2[0, 1]^2 meeting at (0.5,0.5)(0.5, 0.5).

Numerical method

  • Flux: HLLC approximate Riemann solver.
  • Time: explicit Euler.
  • Adaptation: multiresolution with a positivity-preserving prediction.

The image shows the density ρ\rho.

What to look for

Curved shocks bound the interaction region, and a thin, unstable jet develops along the diagonal and rolls up. The adaptive mesh follows every discontinuity, which is what makes the fine-scale roll-up affordable.

Source

scenario.hpp · 228 lines

This case runs on the samurai-euler solver. The code below is its scenario: initial and boundary conditions. Fluxes, time integration and mesh adaptation come from the shared solver. 

// Copyright 2025 the samurai team
// SPDX-License-Identifier:  BSD-3-Clause

#pragma once

#include <numbers>

#include <samurai/bc.hpp>

#include "../variables.hpp"
#include "registry.hpp"

namespace test_case::riemann_2d_config_3
{
    double x0 = 0.5;
    double y0 = 0.5;

    PrimState<2> quad1_state{
        1.5,
        1.5,
        xt::xtensor_fixed<double, xt::xshape<2>>{0., 0.}
    };

    PrimState<2> quad2_state{
        0.5323,
        0.3,
        xt::xtensor_fixed<double, xt::xshape<2>>{1.206, 0.}
    };

    PrimState<2> quad3_state{
        0.138,
        0.29,
        xt::xtensor_fixed<double, xt::xshape<2>>{1.206, 1.206}
    };

    PrimState<2> quad4_state{
        0.5323,
        0.3,
        xt::xtensor_fixed<double, xt::xshape<2>>{0, 1.206}
    };

    auto init_fn = [](auto& u, auto& cell)
    {
        auto x = cell.center();

        if (x[0] >= x0 && x[1] >= y0)
        {
            u[cell] = prim2cons<2>(quad1_state);
        }
        else if (x[0] < x0 && x[1] >= y0)
        {
            u[cell] = prim2cons<2>(quad2_state);
        }
        else if (x[0] < x0 && x[1] < y0)
        {
            u[cell] = prim2cons<2>(quad3_state);
        }
        else // (x[0] >= x0 && x[1] < y0)
        {
            u[cell] = prim2cons<2>(quad4_state);
        }
    };

    void bc_fn(auto& u, double /*t*/)
    {
        samurai::make_bc<samurai::Neumann<1>>(u, 0., 0., 0., 0.);
    }

    template <std::size_t dim>
    auto box_fn()
    {
        xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
        xt::xtensor_fixed<double, xt::xshape<dim>> max_corner = {1., 1.};

        return samurai::Box<double, dim>(min_corner, max_corner);
    }

}

REGISTER_TEST_CASE(riemann2d_config3,
                   test_case::riemann_2d_config_3::box_fn,
                   test_case::riemann_2d_config_3::init_fn,
                   test_case::riemann_2d_config_3::bc_fn)

namespace test_case::riemann_2d_config_4
{
    double x0 = 0.5;
    double y0 = 0.5;

    PrimState<2> quad1_state{
        1.1,
        1.1,
        xt::xtensor_fixed<double, xt::xshape<2>>{0., 0.}
    };

    PrimState<2> quad2_state{
        0.5065,
        0.35,
        xt::xtensor_fixed<double, xt::xshape<2>>{0.8939, 0.}
    };

    PrimState<2> quad3_state{
        1.1,
        1.1,
        xt::xtensor_fixed<double, xt::xshape<2>>{0.8939, 0.89396}
    };

    PrimState<2> quad4_state{
        0.5065,
        0.35,
        xt::xtensor_fixed<double, xt::xshape<2>>{0, 0.89396}
    };

    auto init_fn = [](auto& u, auto& cell)
    {
        auto x = cell.center();

        if (x[0] >= x0 && x[1] >= y0)
        {
            u[cell] = prim2cons<2>(quad1_state);
        }
        else if (x[0] < x0 && x[1] >= y0)
        {
            u[cell] = prim2cons<2>(quad2_state);
        }
        else if (x[0] < x0 && x[1] < y0)
        {
            u[cell] = prim2cons<2>(quad3_state);
        }
        else // (x[0] >= x0 && x[1] < y0)
        {
            u[cell] = prim2cons<2>(quad4_state);
        }
    };

    void bc_fn(auto& u, double /*t*/)
    {
        samurai::make_bc<samurai::Neumann<1>>(u, 0., 0., 0., 0.);
    }

    template <std::size_t dim>
    auto box_fn()
    {
        xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
        xt::xtensor_fixed<double, xt::xshape<dim>> max_corner = {1., 1.};

        return samurai::Box<double, dim>(min_corner, max_corner);
    }

}

REGISTER_TEST_CASE(riemann2d_config4,
                   test_case::riemann_2d_config_4::box_fn,
                   test_case::riemann_2d_config_4::init_fn,
                   test_case::riemann_2d_config_4::bc_fn)

namespace test_case::riemann_2d_config_12
{
    double x0 = 0.5;
    double y0 = 0.5;

    PrimState<2> quad1_state{
        0.5197,
        0.4,
        xt::xtensor_fixed<double, xt::xshape<2>>{0., 0.}
    };

    PrimState<2> quad2_state{
        1,
        1,
        xt::xtensor_fixed<double, xt::xshape<2>>{-0.6259, 0.}
    };

    PrimState<2> quad3_state{
        0.8,
        1,
        xt::xtensor_fixed<double, xt::xshape<2>>{-0.6259, -0.6259}
    };

    PrimState<2> quad4_state{
        1,
        1,
        xt::xtensor_fixed<double, xt::xshape<2>>{0, -0.6259}
    };

    auto init_fn = [](auto& u, auto& cell)
    {
        auto x = cell.center();

        if (x[0] >= x0 && x[1] >= y0)
        {
            u[cell] = prim2cons<2>(quad1_state);
        }
        else if (x[0] < x0 && x[1] >= y0)
        {
            u[cell] = prim2cons<2>(quad2_state);
        }
        else if (x[0] < x0 && x[1] < y0)
        {
            u[cell] = prim2cons<2>(quad3_state);
        }
        else // (x[0] >= x0 && x[1] < y0)
        {
            u[cell] = prim2cons<2>(quad4_state);
        }
    };

    void bc_fn(auto& u, double /*t*/)
    {
        samurai::make_bc<samurai::Neumann<1>>(u, 0., 0., 0., 0.);
    }

    template <std::size_t dim>
    auto box_fn()
    {
        xt::xtensor_fixed<double, xt::xshape<dim>> min_corner = {0., 0.};
        xt::xtensor_fixed<double, xt::xshape<dim>> max_corner = {1., 1.};

        return samurai::Box<double, dim>(min_corner, max_corner);
    }

}

REGISTER_TEST_CASE(riemann2d_config12,
                   test_case::riemann_2d_config_12::box_fn,
                   test_case::riemann_2d_config_12::init_fn,
                   test_case::riemann_2d_config_12::bc_fn)