GammaLawThermo_test.cpp

Go to the documentation of this file.

 
#include "thermo/GammaLawThermo.h"
 
#include <gtest/gtest.h>
 
#include <cmath>
#include <stdexcept>
#include <vector>
 
namespace
{
 
    using rizer::ElectronGamma;
    using rizer::GammaLawThermo;
 
    constexpr double kGasConstant = 8314.462618153242; // [J/kmol/K]
    constexpr double kGamma = 1.4;
    constexpr double kHeavyWeight = 28.96;        // [kg/kmol]
    constexpr double kElectronWeight = 5.485799090e-4; // [kg/kmol]

    GammaLawThermo twoPool()
    {
        return GammaLawThermo(kGamma, {kHeavyWeight, kElectronWeight},
                              /*electron_index=*/1);
    }

    GammaLawThermo neutral()
    {
        return GammaLawThermo(kGamma, {kHeavyWeight}, /*electron_index=*/-1);
    }
 
    // ------------------------------------------------------------------------------ //
    //   pressures                                                                    //
    // ------------------------------------------------------------------------------ //

    TEST(GammaLawThermoPressures, EachPoolUsesItsOwnTemperature)
    {
        const std::vector<double> Y = {1.0 - 1.0e-6, 1.0e-6};
        const double rho = 1.0;
        const double Tg = 3000.0;
        const double Te = 12000.0;
        double p_heavy = 0.0;
        double p_electron = 0.0;
        twoPool().pressures(0.0, 0.0, rho, Tg, Te, Y.data(), p_heavy, p_electron);
 
        EXPECT_DOUBLE_EQ(p_heavy, rho * (Y[0] / kHeavyWeight) * kGasConstant * Tg);
        EXPECT_DOUBLE_EQ(p_electron,
                         rho * (Y[1] / kElectronWeight) * kGasConstant * Te);
    }

    TEST(GammaLawThermoPressures, TheElectronPartialPressureIsNotNegligible)
    {
        const std::vector<double> Y = {1.0 - 1.0e-6, 1.0e-6};
        double p_heavy = 0.0;
        double p_electron = 0.0;
        twoPool().pressures(0.0, 0.0, 1.0, 3000.0, 12000.0, Y.data(), p_heavy,
                            p_electron);
        EXPECT_GT(p_electron / p_heavy, 0.1);
    }

    TEST(GammaLawThermoPressures, ANeutralMixtureHasNoElectronPressure)
    {
        const std::vector<double> Y = {1.0};
        double p_heavy = 0.0;
        double p_electron = 0.0;
        neutral().pressures(0.0, 0.0, 1.2, 300.0, 300.0, Y.data(), p_heavy,
                            p_electron);
        EXPECT_DOUBLE_EQ(p_electron, 0.0);
        EXPECT_GT(p_heavy, 0.0);
    }
 
    // ------------------------------------------------------------------------------ //
    //   soundSpeed                                                                   //
    // ------------------------------------------------------------------------------ //

    TEST(GammaLawThermoSoundSpeed, ReducesToTheTextbookFormWithNoElectron)
    {
        const std::vector<double> Y = {1.0};
        const double expected =
            std::sqrt(kGamma * kGasConstant * 300.0 / kHeavyWeight);
        EXPECT_DOUBLE_EQ(neutral().soundSpeed(0.0, 0.0, 1.2, 300.0, 300.0, Y.data()),
                         expected);
    }

    TEST(GammaLawThermoSoundSpeed, AddsTheTwoPoolsWithTheirOwnGammas)
    {
        const std::vector<double> Y = {1.0 - 1.0e-6, 1.0e-6};
        const double rho = 1.0;
        const double Tg = 3000.0;
        const double Te = 12000.0;
        const GammaLawThermo thermo = twoPool();
        double p_heavy = 0.0;
        double p_electron = 0.0;
        thermo.pressures(0.0, 0.0, rho, Tg, Te, Y.data(), p_heavy, p_electron);
 
        const double expected =
            std::sqrt((kGamma * p_heavy + ElectronGamma * p_electron) / rho);
        EXPECT_NEAR(thermo.soundSpeed(0.0, 0.0, rho, Tg, Te, Y.data()), expected,
                    1.0e-9 * expected);
    }

    TEST(GammaLawThermoSoundSpeed, RefusesAMasslessComposition)
    {
        const std::vector<double> Y = {0.0};
        EXPECT_THROW(neutral().soundSpeed(0.0, 0.0, 1.2, 300.0, 300.0, Y.data()),
                     std::invalid_argument);
    }
 
    // ------------------------------------------------------------------------------ //
    //   Inversions                                                                   //
    // ------------------------------------------------------------------------------ //

    TEST(GammaLawThermoInversion, DensityFromPressureInvertsPressuresExactly)
    {
        const std::vector<double> Y = {1.0 - 1.0e-6, 1.0e-6};
        const GammaLawThermo thermo = twoPool();
        const double Tg = 3000.0;
        const double Te = 12000.0;
        const double rho = 0.37;
        double p_heavy = 0.0;
        double p_electron = 0.0;
        thermo.pressures(0.0, 0.0, rho, Tg, Te, Y.data(), p_heavy, p_electron);
 
        const double recovered =
            thermo.densityFromPressure(p_heavy + p_electron, Tg, Te, Y.data());
        EXPECT_NEAR(recovered, rho, 1.0e-12 * rho);
    }

    TEST(GammaLawThermoInversion, HeavyTemperatureInvertsHeavyEnergyIgnoringTheSeed)
    {
        const std::vector<double> Y = {1.0 - 1.0e-6, 1.0e-6};
        const GammaLawThermo thermo = twoPool();
        const double Tg = 4250.0;
        const double energy = thermo.heavyEnergy(Tg, Y.data());
 
        for (const double seed : {1.0, 300.0, 50000.0})
        {
            EXPECT_NEAR(thermo.heavyTemperature(1.0, Y.data(), energy, seed), Tg,
                        1.0e-9)
                << "seed = " << seed;
        }
    }

    TEST(GammaLawThermoInversion, HeavyEnergyIsPerKilogramOfMixture)
    {
        const GammaLawThermo thermo = twoPool();
        const std::vector<double> full = {1.0, 0.0};
        const std::vector<double> half = {0.5, 0.5};
        EXPECT_DOUBLE_EQ(thermo.heavyEnergy(3000.0, half.data()),
                         0.5 * thermo.heavyEnergy(3000.0, full.data()));
    }
 
    // ------------------------------------------------------------------------------ //
    //   electronCv                                                                   //
    // ------------------------------------------------------------------------------ //

    TEST(GammaLawThermoElectronCv, IsTheMonatomicValue)
    {
        EXPECT_NEAR(twoPool().electronCv(),
                    1.5 * kGasConstant / kElectronWeight,
                    1.0e-6 * (1.5 * kGasConstant / kElectronWeight));
    }

    TEST(GammaLawThermoElectronCv, RefusesWhenThereIsNoElectronSpecies)
    {
        EXPECT_THROW(neutral().electronCv(), std::invalid_argument);
    }
 
} // namespace