Add the five-equation two-phase model

CMakeLists.txt

@@ -31,6 +31,14 @@ target_link_libraries(euler_2d PRIVATE samurai::samurai common_defs)
 add_executable(euler_3d main_3d.cpp)
 target_link_libraries(euler_3d PRIVATE samurai::samurai common_defs)
 
+# The five-equation two-phase model. One time loop, written once and instantiated
+# per dimension, rather than one main per dimension as above.
+add_executable(two_phase_1d main_two_phase_1d.cpp)
+target_link_libraries(two_phase_1d PRIVATE samurai::samurai common_defs)
+
+add_executable(two_phase_2d main_two_phase_2d.cpp)
+target_link_libraries(two_phase_2d PRIVATE samurai::samurai common_defs)
+
 # ---------------------------------------------------------------------------
 # Tests
 # ---------------------------------------------------------------------------

euler/bc.hpp

@@ -48,6 +48,16 @@ namespace detail
     {
         return stencil_size / 2;
     }
+
+    // Where the momentum sits, in any of the models this repository solves: the
+    // monofluid Euler system and the five-equation two-phase one put it at the
+    // same place on purpose, and each layout holds itself to it with a
+    // static_assert. The solid wall is then one boundary condition rather than
+    // one per model, since mirroring the normal momentum is all it does.
+    inline constexpr std::size_t momentum(std::size_t d)
+    {
+        return 2 + d;
+    }
 }
 
 template <std::size_t StencilSize, class Field>
@@ -98,8 +108,8 @@ struct ReflectiveImpl : public samurai::Bc<Field>
                 const auto& inside = cells[ghost0 - 1 - k];
                 const auto& ghost  = cells[ghost0 + k];
 
-                u[ghost]                                       = u[inside];
-                u[ghost][EulerLayout<Field::dim>::mom(normal)] = -u[inside][EulerLayout<Field::dim>::mom(normal)];
+                u[ghost]                           = u[inside];
+                u[ghost][detail::momentum(normal)] = -u[inside][detail::momentum(normal)];
             }
         };
     }
@@ -209,15 +219,22 @@ namespace bc
         return wide() ? samurai::make_bc<Imposed<4>>(u, f) : samurai::make_bc<Imposed<2>>(u, f);
     }
 
-    // A uniform state imposed on every boundary.
-    template <class Field, class Eos>
-    auto imposed(Field& u, const PrimState<Field::dim>& state, Eos eos)
+    // A uniform state imposed on every boundary, given as the conservative
+    // vector itself. Which components those are is the model's business, so this
+    // is the one both of them go through.
+    template <class Field, class Array>
+    auto imposed_state(Field& u, const Array& cons)
     {
-        const auto cons = prim2cons<Field::dim>(state, eos);
-
         return [&]<std::size_t... I>(std::index_sequence<I...>)
         {
             return wide() ? samurai::make_bc<Imposed<4>>(u, cons[I]...) : samurai::make_bc<Imposed<2>>(u, cons[I]...);
         }(std::make_index_sequence<Field::n_comp>{});
     }
+
+    // The same, from a monofluid primitive state.
+    template <class Field, class Eos>
+    auto imposed(Field& u, const PrimState<Field::dim>& state, Eos eos)
+    {
+        return imposed_state(u, prim2cons<Field::dim>(state, eos));
+    }
 }

euler/reconstruction.hpp

@@ -99,12 +99,15 @@ double limited_slope(double dm, double dp)
     }
 }
 
-// The same, component by component, on a packed primitive state.
-template <SlopeLimiter limiter, std::size_t Dim, class Array>
-ConsArray<Dim> limited_slope(const Array& dm, const Array& dp)
+// The same, component by component, on a packed state of N components. N rather
+// than the dimension because both models go through here: the monofluid state
+// has dim + 2 components and the two-phase one dim + 4, and a limiter has no
+// opinion on which component is which.
+template <SlopeLimiter limiter, std::size_t N, class Array>
+xt::xtensor_fixed<double, xt::xshape<N>> limited_slope(const Array& dm, const Array& dp)
 {
-    ConsArray<Dim> slope;
-    for (std::size_t i = 0; i < EulerLayout<Dim>::size; ++i)
+    xt::xtensor_fixed<double, xt::xshape<N>> slope;
+    for (std::size_t i = 0; i < N; ++i)
     {
         slope[i] = limited_slope<limiter>(dm[i], dp[i]);
     }
@@ -162,18 +165,18 @@ ConsArray<Dim> primitive_jacobian_times(const PrimState<Dim>& prim, const Array&
 // costs one branch per slope, against one Riemann solve per interface: the
 // limiter stays a command line option without being a template parameter of
 // every scheme.
-template <std::size_t Dim, class Array>
-ConsArray<Dim> limited_slope(const Array& dm, const Array& dp, SlopeLimiter limiter)
+template <std::size_t N, class Array>
+xt::xtensor_fixed<double, xt::xshape<N>> limited_slope(const Array& dm, const Array& dp, SlopeLimiter limiter)
 {
     switch (limiter)
     {
         case SlopeLimiter::none:
-            return limited_slope<SlopeLimiter::none, Dim>(dm, dp);
+            return limited_slope<SlopeLimiter::none, N>(dm, dp);
         case SlopeLimiter::minmod:
-            return limited_slope<SlopeLimiter::minmod, Dim>(dm, dp);
+            return limited_slope<SlopeLimiter::minmod, N>(dm, dp);
         case SlopeLimiter::vanleer:
-            return limited_slope<SlopeLimiter::vanleer, Dim>(dm, dp);
+            return limited_slope<SlopeLimiter::vanleer, N>(dm, dp);
         default:
-            return limited_slope<SlopeLimiter::moncen, Dim>(dm, dp);
+            return limited_slope<SlopeLimiter::moncen, N>(dm, dp);
     }
 }

euler/schemes/muscl.hpp

@@ -89,8 +89,10 @@ auto make_muscl_scheme(Eos eos, const MusclOptions& options)
                 const auto w2 = pack<dim>(cons2prim<dim>(field[2], eos));
                 const auto w3 = pack<dim>(cons2prim<dim>(field[3], eos));
 
-                const auto slopeL = limited_slope<dim>(w1 - w0, w2 - w1, options.limiter);
-                const auto slopeR = limited_slope<dim>(w2 - w1, w3 - w2, options.limiter);
+                static constexpr std::size_t n_comp = EulerLayout<dim>::size;
+
+                const auto slopeL = limited_slope<n_comp>(w1 - w0, w2 - w1, options.limiter);
+                const auto slopeR = limited_slope<n_comp>(w2 - w1, w3 - w2, options.limiter);
 
                 ConsArray<dim> wL = w1 + 0.5 * slopeL;
                 ConsArray<dim> wR = w2 - 0.5 * slopeR;

euler/two_phase/eos.hpp

@@ -0,0 +1,99 @@
+// Copyright 2025 the samurai team
+// SPDX-License-Identifier:  BSD-3-Clause
+
+#pragma once
+
+#include <cmath>
+
+#include "../eos.hpp"
+
+// =============================================================================
+//  Mixture equation of state of the five-equation model
+// -----------------------------------------------------------------------------
+//  Two stiffened gases sharing a cell, in pressure and velocity equilibrium.
+//  Their volume fractions alpha_0 and alpha_1 = 1 - alpha_0 say how much of the
+//  cell each one occupies, and the mixture behaves as a single stiffened gas
+//  whose coefficients depend on alpha:
+//
+//      1 / (gamma_m - 1) = sum_i alpha_i / (gamma_i - 1)                    (G)
+//      gamma_m pi_m / (gamma_m - 1) = sum_i alpha_i gamma_i pi_i / (gamma_i - 1)
+//
+//  Those two sums are the whole equation of state. Written as G(alpha) and
+//  P(alpha) they give the internal energy and the pressure directly,
+//
+//      rho e = G p + P          p = (rho e - P) / G
+//
+//  and the sound speed is the stiffened-gas one at the mixture coefficients,
+//
+//      c^2 = gamma_m (p + pi_m) / rho,   gamma_m = 1 + 1/G,  pi_m = P / (G gamma_m)
+//
+//  which is the frozen speed the Allaire model is hyperbolic for. A single phase
+//  (alpha = 0 or 1) gives back exactly its own gamma and pi_inf, which is worth
+//  keeping in mind: a two-phase run on a pure fluid must reproduce the monofluid
+//  solver bit for bit up to the scheme, and that is a test.
+//
+//      G. Allaire, S. Clerc, S. Kokh, "A five-equation model for the simulation
+//      of interfaces between compressible fluids", J. Comput. Phys. 181 (2)
+//      (2002) 577-616, https://doi.org/10.1006/jcph.2002.7143
+// =============================================================================
+
+namespace EOS
+{
+    struct Mixture
+    {
+        StiffenedGas phase[2] = {
+            {1.4, 0., 0.},
+            {1.4, 0., 0.}
+        };
+
+        // sum_i alpha_i / (gamma_i - 1), the inverse of gamma_m - 1.
+        double G(double alpha0) const
+        {
+            const double alpha1 = 1. - alpha0;
+            return alpha0 / (phase[0].gamma - 1.) + alpha1 / (phase[1].gamma - 1.);
+        }
+
+        // sum_i alpha_i gamma_i pi_i / (gamma_i - 1).
+        double P(double alpha0) const
+        {
+            const double alpha1 = 1. - alpha0;
+            return alpha0 * phase[0].gamma * phase[0].pi_inf / (phase[0].gamma - 1.)
+                 + alpha1 * phase[1].gamma * phase[1].pi_inf / (phase[1].gamma - 1.);
+        }
+
+        double gamma(double alpha0) const
+        {
+            return 1. + 1. / G(alpha0);
+        }
+
+        double pi_inf(double alpha0) const
+        {
+            const double g = G(alpha0);
+            return P(alpha0) / (g * (1. + 1. / g));
+        }
+
+        // Pressure from the volumic internal energy rho e.
+        double p(double alpha0, double rho_e) const
+        {
+            return (rho_e - P(alpha0)) / G(alpha0);
+        }
+
+        // The inverse: the volumic internal energy at that pressure.
+        double rho_e(double alpha0, double p) const
+        {
+            return G(alpha0) * p + P(alpha0);
+        }
+
+        double c(double alpha0, double rho, double p) const
+        {
+            return std::sqrt(gamma(alpha0) * (p + pi_inf(alpha0)) / rho);
+        }
+    };
+
+    inline constexpr Mixture two_ideal_gases(double gamma0, double gamma1)
+    {
+        return {
+            {{gamma0, 0., 0.}, {gamma1, 0., 0.}}
+        };
+    }
+}

euler/two_phase/init/advected_interface.hpp

@@ -0,0 +1,93 @@
+// Copyright 2025 the samurai team
+// SPDX-License-Identifier:  BSD-3-Clause
+
+#pragma once
+
+#include <samurai/box.hpp>
+
+#include "../registry.hpp"
+#include "../variables.hpp"
+
+// =============================================================================
+//  An interface carried by a uniform flow
+// -----------------------------------------------------------------------------
+//  A slab of water in air, everything at one atmosphere and moving at the same
+//  speed. The exact solution is the initial state translated: the pressure and
+//  the velocity stay uniform for ever, and only the volume fraction moves.
+//
+//  This is the test the five-equation model is designed to pass and the reason
+//  it exists. A model that advects the volume fraction inconsistently with the
+//  masses and the energy produces a pressure spike at the interface out of
+//  nothing -- the classic failure of a conservative two-material scheme, where
+//  averaging two gases in one cell gives a mixture whose pressure is not the
+//  pressure of either. The interface here is between fluids whose stiffened-gas
+//  coefficients could hardly be further apart, so nothing hides.
+//
+//      P. Kapila et al. / R. Abgrall, "How to prevent pressure oscillations in
+//      multicomponent flow calculations: a quasi conservative approach",
+//      J. Comput. Phys. 125 (1996) 150-160,
+//      https://doi.org/10.1006/jcph.1996.0085
+//
+//  Domain [0,1], periodic, water in [0.3, 0.7], p = 1e5 and u = 100 everywhere.
+// =============================================================================
+
+namespace test_case::advected_interface
+{
+    inline constexpr double left  = 0.3;
+    inline constexpr double right = 0.7;
+
+    inline constexpr double pressure = 1e5;
+    inline constexpr double speed    = 100.;
+
+    inline constexpr double water_rho = 1e3;
+    inline constexpr double air_rho   = 1.;
+
+    inline const EOS::Mixture eos{
+        {{4.4, 6e8, 0.}, {1.4, 0., 0.}}
+    };
+
+    template <std::size_t dim>
+    auto box_fn()
+    {
+        xt::xtensor_fixed<double, xt::xshape<dim>> min_corner;
+        xt::xtensor_fixed<double, xt::xshape<dim>> max_corner;
+        min_corner.fill(0.);
+        max_corner.fill(1.);
+
+        return samurai::Box<double, dim>(min_corner, max_corner);
+    }
+
+    template <class Field>
+    void init_fn(Field& u, const typename Field::cell_t& cell, EOS::Mixture eos_)
+    {
+        static constexpr std::size_t dim = Field::dim;
+
+        const auto x     = cell.center();
+        const bool water = x[0] > left && x[0] < right;
+
+        auto state = water ? two_phase::mixture_state<dim>(1., water_rho, 0., pressure)
+                           : two_phase::mixture_state<dim>(0., 0., air_rho, pressure);
+        state.v[0] = speed;
+
+        u[cell] = two_phase::prim2cons<dim>(state, eos_);
+    }
+
+    template <class Field>
+    void bc_fn(Field& /*u*/, double& /*t*/, EOS::Mixture /*eos*/)
+    {
+    }
+
+    template <class Field>
+    two_phase::test_case_t<Field::dim> definition()
+    {
+        auto definition     = two_phase::test_case_t<Field::dim>{.box  = &box_fn<Field::dim>,
+                                                                 .init = &init_fn<Field>,
+                                                                 .bc   = &bc_fn<Field>,
+                                                                 .eos  = eos};
+        definition.periodic = {};
+        definition.periodic.fill(true);
+        return definition;
+    }
+}
+
+REGISTER_TWO_PHASE_CASE(advected_interface, test_case::advected_interface, 1, 2)

euler/two_phase/init/cases.hpp

@@ -0,0 +1,13 @@
+// Copyright 2025 the samurai team
+// SPDX-License-Identifier:  BSD-3-Clause
+
+#pragma once
+
+// Every two-phase test case, in one include. A case registers itself, so adding
+// one is a file and a line here.
+
+#include "advected_interface.hpp"
+#include "shock_bubble.hpp"
+#include "sod_x_pure.hpp"
+#include "triple_point.hpp"
+#include "water_air_shock_tube.hpp"