2D Riemann problem (config 3)
Four constant states meet at a point and interact through shocks, contacts and a rolled-up jet.
2D Riemann problem (configuration 3)
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 (
Numerical method
- Flux: HLLC approximate Riemann solver.
- Time: explicit Euler.
- Adaptation: multiresolution with a positivity-preserving prediction.
The image shows the density
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 code scenario.hpp
Powered by the samurai-euler solver engine - the code below is this case's scenario (initial & boundary conditions); the flux, time integration and mesh adaptation come from the shared engine.
// 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)