HydroStep_test.cpp

Go to the documentation of this file.

 
#include "models/nrp/HydroStep.h"
#include "thermo/GammaLawThermo.h"
 
#include <gtest/gtest.h>
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <vector>
 
namespace
{
 
    using rizer::GammaLawThermo;
    using rizer::HydroStepConfig;
    using rizer::HydroStepState;
    using rizer::NodeState;
 
    constexpr double kGamma = 1.4;
    constexpr double kMolecularWeight = 28.96; // [kg/kmol]
    constexpr double kPressure = 101325.0;     // [Pa]
    constexpr double kTemperature = 300.0;     // [K]
    constexpr double kGasConstant = 8314.462618153242; // [J/kmol/K]

    class RecordingThermo : public rizer::ThermoAccessor
    {
    public:
        RecordingThermo()
            : m_backing(kGamma, {kMolecularWeight}, /*electron_index=*/-1)
        {
        }
 
        void pressures(double r, double t, double rho, double Tg, double Te,
                       const double *Y, double &p_heavy,
                       double &p_electron) const override
        {
            radii.push_back(r);
            m_backing.pressures(r, t, rho, Tg, Te, Y, p_heavy, p_electron);
        }
 
        double soundSpeed(double r, double t, double rho, double Tg, double Te,
                          const double *Y) const override
        {
            return m_backing.soundSpeed(r, t, rho, Tg, Te, Y);
        }
 
        double heavyEnergy(double Tg, const double *Y) const override
        {
            return m_backing.heavyEnergy(Tg, Y);
        }
 
        double electronCv() const override { return m_backing.electronCv(); }
 
        double densityFromPressure(double pressure, double Tg, double Te,
                                   const double *Y) const override
        {
            return m_backing.densityFromPressure(pressure, Tg, Te, Y);
        }
 
        double heavyTemperature(double rho, const double *Y,
                                double heavy_internal_energy_target,
                                double Tg_seed) const override
        {
            return m_backing.heavyTemperature(rho, Y, heavy_internal_energy_target,
                                              Tg_seed);
        }
 
        mutable std::vector<double> radii;
 
    private:
        GammaLawThermo m_backing;
    };

    struct Grid
    {
        std::vector<double> r_node;
        std::vector<double> r_face;
        std::vector<double> V_cell;
    };

    Grid makeGrid(std::size_t n_points, double r_max, double r_min = 0.0,
                  double stretch = 1.0)
    {
        Grid grid;
        grid.r_node.resize(n_points);
        for (std::size_t i = 0; i < n_points; i++)
        {
            const double x =
                static_cast<double>(i) / static_cast<double>(n_points - 1);
            grid.r_node[i] = r_min + (r_max - r_min) * std::pow(x, stretch);
        }
        grid.r_face.resize(n_points + 1);
        grid.r_face[0] = (r_min == 0.0)
                             ? 0.0
                             : grid.r_node[0] -
                                   0.5 * (grid.r_node[1] - grid.r_node[0]);
        for (std::size_t i = 1; i < n_points; i++)
        {
            grid.r_face[i] = 0.5 * (grid.r_node[i - 1] + grid.r_node[i]);
        }
        grid.r_face[n_points] =
            grid.r_node[n_points - 1] +
            0.5 * (grid.r_node[n_points - 1] - grid.r_node[n_points - 2]);
        grid.V_cell.resize(n_points);
        for (std::size_t i = 0; i < n_points; i++)
        {
            grid.V_cell[i] = 0.5 * (grid.r_face[i + 1] * grid.r_face[i + 1] -
                                    grid.r_face[i] * grid.r_face[i]);
        }
        return grid;
    }

    HydroStepState makeState(std::size_t n_points)
    {
        HydroStepState state;
        state.n_species = 1;
        const double rho = kPressure * kMolecularWeight / (kGasConstant * kTemperature);
        state.rho.assign(n_points, rho);
        state.u_r.assign(n_points, 0.0);
        state.Tg.assign(n_points, kTemperature);
        state.Te.assign(n_points, kTemperature);
        state.Y.assign(n_points, 1.0);
        return state;
    }

    HydroStepConfig makeConfig(const Grid &grid, const rizer::ThermoAccessor &thermo)
    {
        HydroStepConfig cfg;
        cfg.r_node = grid.r_node;
        cfg.r_face = grid.r_face;
        cfg.V_cell = grid.V_cell;
        cfg.cylindrical_source = true;
        cfg.p_amb = kPressure;
        cfg.T_amb = kTemperature;
        cfg.Y_amb = {1.0};
        cfg.dt_total = 1.0e-12;
        cfg.n_subcycles = 1;
        cfg.electron_index = -1;
        cfg.thermo = &thermo;
        return cfg;
    }
 
    // ------------------------------------------------------------------------------ //
    //   Ghost placement                                                              //
    // ------------------------------------------------------------------------------ //

    TEST(HydroStepGhostRadius, OuterGhostSitsAtItsOwnCellNode)
    {
        const Grid grid = makeGrid(21, 4.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
 
        std::vector<double> offered;
        cfg.ghost_state = [&](double r_ghost, double, bool is_inner,
                              const NodeState &interior) {
            EXPECT_FALSE(is_inner) << "a full disk must not consult the axis face";
            offered.push_back(r_ghost);
            return interior;
        };
 
        HydroStepState state = makeState(21);
        rizer::hydroStep(cfg, state);
 
        const double expected = 2.0 * grid.r_node[20] - grid.r_node[19];
        ASSERT_FALSE(offered.empty());
        EXPECT_DOUBLE_EQ(offered.front(), expected);
        // Half a spacing away from the face, which is what it must NOT be.
        EXPECT_GT(std::abs(offered.front() - grid.r_face[21]), 1.0e-5);
    }

    TEST(HydroStepGhostRadius, BothGhostsAreOfferedOnAnAnnulus)
    {
        const Grid grid = makeGrid(21, 4.0e-3, /*r_min=*/1.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
 
        double inner = 0.0;
        double outer = 0.0;
        cfg.ghost_state = [&](double r_ghost, double, bool is_inner,
                              const NodeState &interior) {
            (is_inner ? inner : outer) = r_ghost;
            return interior;
        };
 
        HydroStepState state = makeState(21);
        rizer::hydroStep(cfg, state);
 
        EXPECT_DOUBLE_EQ(inner, 2.0 * grid.r_node[0] - grid.r_node[1]);
        EXPECT_DOUBLE_EQ(outer, 2.0 * grid.r_node[20] - grid.r_node[19]);
    }

    TEST(HydroStepGhostRadius, EachEndUsesItsOwnSpacing)
    {
        const Grid grid = makeGrid(21, 4.0e-3, /*r_min=*/1.0e-3, /*stretch=*/2.5);
        const double inner_spacing = grid.r_node[1] - grid.r_node[0];
        const double outer_spacing = grid.r_node[20] - grid.r_node[19];
        // The premise of the test: the mesh really is lopsided.
        ASSERT_GT(outer_spacing / inner_spacing, 5.0);
 
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        double inner = 0.0;
        double outer = 0.0;
        cfg.ghost_state = [&](double r_ghost, double, bool is_inner,
                              const NodeState &interior) {
            (is_inner ? inner : outer) = r_ghost;
            return interior;
        };
 
        HydroStepState state = makeState(21);
        rizer::hydroStep(cfg, state);
 
        EXPECT_DOUBLE_EQ(inner, grid.r_node[0] - inner_spacing);
        EXPECT_DOUBLE_EQ(outer, grid.r_node[20] + outer_spacing);
    }

    TEST(HydroStepGhostRadius, ThermoSeesTheSameGhostRadiiAsTheCallback)
    {
        const Grid grid = makeGrid(21, 4.0e-3, /*r_min=*/1.0e-3, /*stretch=*/2.5);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        double inner = 0.0;
        double outer = 0.0;
        cfg.ghost_state = [&](double r_ghost, double, bool is_inner,
                              const NodeState &interior) {
            (is_inner ? inner : outer) = r_ghost;
            return interior;
        };
 
        HydroStepState state = makeState(21);
        rizer::hydroStep(cfg, state);
 
        ASSERT_FALSE(thermo.radii.empty());
        const auto bounds = std::minmax_element(thermo.radii.begin(), thermo.radii.end());
        EXPECT_DOUBLE_EQ(*bounds.first, inner);
        EXPECT_DOUBLE_EQ(*bounds.second, outer);
    }

    TEST(HydroStepGhostRadius, TheAxisGhostStaysOnTheAxis)
    {
        const Grid grid = makeGrid(21, 4.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        cfg.outer_bc_mode = rizer::OuterBc::AmbientGhost;
 
        HydroStepState state = makeState(21);
        rizer::hydroStep(cfg, state);
 
        ASSERT_FALSE(thermo.radii.empty());
        const double smallest =
            *std::min_element(thermo.radii.begin(), thermo.radii.end());
        // Exactly the axis, not merely non-negative: a ghost placed one spacing
        // below zero would be a negative radius, but one placed at the first
        // face would be positive and still wrong.
        EXPECT_DOUBLE_EQ(smallest, 0.0);
    }
 
    // ------------------------------------------------------------------------------ //
    //   Input guards                                                                 //
    // ------------------------------------------------------------------------------ //

    TEST(HydroStepGuards, RefusesASingleNode)
    {
        Grid grid;
        grid.r_node = {1.0e-3};
        grid.r_face = {0.5e-3, 1.5e-3};
        grid.V_cell = {1.0e-6};
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        HydroStepState state = makeState(1);
        EXPECT_THROW(rizer::hydroStep(cfg, state), std::invalid_argument);
    }

    TEST(HydroStepGuards, RefusesANullThermo)
    {
        const Grid grid = makeGrid(21, 4.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        cfg.thermo = nullptr;
        HydroStepState state = makeState(21);
        EXPECT_THROW(rizer::hydroStep(cfg, state), std::invalid_argument);
    }

    TEST(HydroStepGuards, RefusesACflOutsideTheUnitInterval)
    {
        const Grid grid = makeGrid(21, 4.0e-3);
        RecordingThermo thermo;
        HydroStepState state = makeState(21);
        for (const double cfl : {0.0, -0.5, 1.5})
        {
            HydroStepConfig cfg = makeConfig(grid, thermo);
            cfg.n_subcycles = 0; // so the CFL bound is actually derived
            cfg.cfl = cfl;
            EXPECT_THROW(rizer::hydroStep(cfg, state), std::invalid_argument)
                << "cfl = " << cfl;
        }
    }

    TEST(HydroStepGuards, RefusesAnAnnulusWithNoGhostCallback)
    {
        const Grid grid = makeGrid(21, 4.0e-3, /*r_min=*/1.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        HydroStepState state = makeState(21);
        EXPECT_THROW(rizer::hydroStep(cfg, state), std::invalid_argument);
    }

    TEST(HydroStepGuards, RefusesAGhostWhoseCompositionDoesNotSumToOne)
    {
        const Grid grid = makeGrid(21, 4.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        cfg.ghost_state = [](double, double, bool, const NodeState &interior) {
            NodeState ghost = interior;
            ghost.Y[0] = 0.5;
            return ghost;
        };
        HydroStepState state = makeState(21);
        EXPECT_THROW(rizer::hydroStep(cfg, state), std::invalid_argument);
    }

    TEST(HydroStepGuards, RefusesANonFiniteGhost)
    {
        const Grid grid = makeGrid(21, 4.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        cfg.ghost_state = [](double, double, bool, const NodeState &interior) {
            NodeState ghost = interior;
            ghost.u_r = std::nan("");
            return ghost;
        };
        HydroStepState state = makeState(21);
        EXPECT_THROW(rizer::hydroStep(cfg, state), std::invalid_argument);
    }

    TEST(HydroStepGuards, AZeroLengthCallAdvancesNothing)
    {
        const Grid grid = makeGrid(21, 4.0e-3);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        cfg.dt_total = 0.0;
        cfg.n_subcycles = 0;
        HydroStepState state = makeState(21);
        const HydroStepState before = state;
        const rizer::HydroStepResult result = rizer::hydroStep(cfg, state);
        EXPECT_EQ(result.n_subcycles_used, 0u);
        EXPECT_EQ(state.rho, before.rho);
        EXPECT_EQ(state.u_r, before.u_r);
        EXPECT_EQ(state.Tg, before.Tg);
    }
 
    // ------------------------------------------------------------------------------ //
    //   Well-balancedness                                                            //
    // ------------------------------------------------------------------------------ //

    TEST(HydroStepWellBalanced, AStateAtRestStaysAtRestOnAStretchedMesh)
    {
        const Grid grid = makeGrid(41, 4.0e-3, /*r_min=*/0.0, /*stretch=*/2.0);
        RecordingThermo thermo;
        HydroStepConfig cfg = makeConfig(grid, thermo);
        cfg.outer_bc_mode = rizer::OuterBc::AmbientGhost;
        cfg.dt_total = 1.0e-9;
        cfg.n_subcycles = 0;
 
        HydroStepState state = makeState(41);
        rizer::hydroStep(cfg, state);
 
        for (std::size_t i = 0; i < state.u_r.size(); i++)
        {
            EXPECT_NEAR(state.u_r[i], 0.0, 1.0e-14) << "node " << i;
        }
    }
 
} // namespace