IdealGasMixtureFluidPropertiesTest.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "IdealGasMixtureFluidPropertiesTest.h"
#include "VaporMixtureFluidPropertiesUtils.h"
#include "FluidProperties.h"
 
TEST_F(IdealGasMixtureFluidPropertiesTest, test)
{
  Real T = 400;
  Real p = 100000;
  std::vector<Real> x = {0.7};
 
  const Real v = 1.38490747936373;
  const Real rho = 0.7220698962933115;
  const Real e = 1.995566878921797e+06;
  const Real s = 3362.64537354375;
  const Real c = 442.127817187648;
  const Real cp = 1187.64406714542;
  const Real cv = 841.417197304492;
  const Real mu = 1.863687329029527e-05;
  const Real k = 0.0304142253602899;
 
  REL_TEST(_fp_mix->p_from_v_e(v, e, x), p, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->T_from_v_e(v, e, x), T, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->v_from_p_T(p, T, x), v, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->rho_from_p_T(p, T, x), rho, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->e_from_p_T(p, T, x), e, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->e_from_p_rho(p, rho, x), e, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->s_from_p_T(p, T, x), s, REL_TOL_SAVED_VALUE);
 
  // The mixture entropy exceeds the mass-weighted component entropies at the mixture pressure
  // by the ideal entropy of mixing, -R / M sum_i psi_i ln(psi_i), with psi_i the mole fractions
  {
    const Real M_steam = 0.01801488;
    const Real M_nitrogen = 0.028012734746133888;
    const Real M_mix = 1.0 / ((1.0 - x[0]) / M_steam + x[0] / M_nitrogen);
    const Real psi_steam = (1.0 - x[0]) * M_mix / M_steam;
    const Real psi_nitrogen = x[0] * M_mix / M_nitrogen;
    const Real s_unmixed =
        (1.0 - x[0]) * _fp_steam->s_from_p_T(p, T) + x[0] * _fp_nitrogen->s_from_p_T(p, T);
    const Real s_mixing = -FluidProperties::_R / M_mix *
                          (psi_steam * std::log(psi_steam) + psi_nitrogen * std::log(psi_nitrogen));
    REL_TEST(_fp_mix->s_from_p_T(p, T, x), s_unmixed + s_mixing, REL_TOL_CONSISTENCY);
  }
  REL_TEST(_fp_mix->c_from_p_T(p, T, x), c, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->cp_from_p_T(p, T, x), cp, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->cv_from_p_T(p, T, x), cv, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->mu_from_p_T(p, T, x), mu, REL_TOL_SAVED_VALUE);
  REL_TEST(_fp_mix->k_from_p_T(p, T, x), k, REL_TOL_SAVED_VALUE);
 
  // Test sound speed against thermodynamics definition:
  //   c = sqrt((dp/drho)_s) = v sqrt(-(dp/dv)_s)
  // To compute (dp/dv)_s, we will use the following identity:
  //   (dp/dv)_s = 1 / [(dv/dp)_T - (dv/dT)_p (ds/dp)_T / (ds/dT)_p]
  const Real rel_pert = REL_PERTURBATION;
  const Real fpert1 = 1 + rel_pert;
  const Real fpert2 = 1 - rel_pert;
  const Real dvdp_T =
      (_fp_mix->v_from_p_T(p * fpert1, T, x) - _fp_mix->v_from_p_T(p * fpert2, T, x)) /
      (2 * rel_pert * p);
  const Real dvdT_p =
      (_fp_mix->v_from_p_T(p, T * fpert1, x) - _fp_mix->v_from_p_T(p, T * fpert2, x)) /
      (2 * rel_pert * T);
  const Real dsdp_T =
      (_fp_mix->s_from_p_T(p * fpert1, T, x) - _fp_mix->s_from_p_T(p * fpert2, T, x)) /
      (2 * rel_pert * p);
  const Real dsdT_p =
      (_fp_mix->s_from_p_T(p, T * fpert1, x) - _fp_mix->s_from_p_T(p, T * fpert2, x)) /
      (2 * rel_pert * T);
  const Real dpdv_s = 1.0 / (dvdp_T - dvdT_p * dsdp_T / dsdT_p);
  const Real c_gold = v * std::sqrt(-dpdv_s);
  REL_TEST(_fp_mix->c_from_p_T(p, T, x), c_gold, REL_TOL_DERIVATIVE);
 
  const Real x_steam = 0.3;
  const Real x_nitrogen = 0.7;
 
  const Real M_steam = 0.01801488;
  const Real M_nitrogen = 0.028012734746133888;
  const Real M = 1.0 / (x_steam / M_steam + x_nitrogen / M_nitrogen);
 
  const Real v_gold = FluidProperties::_R * T / (M * p);
  REL_TEST(v, v_gold, REL_TOL_SAVED_VALUE);
}
 
namespace
{
template <typename F>
void
expectADDerivatives(const F & f, const Real a, const Real b, const Real x)
{
  DNDerivativeType da, db, dx;
  Moose::derivInsert(da, 0, 1);
  Moose::derivInsert(db, 1, 1);
  Moose::derivInsert(dx, 2, 1);
  const ADReal f_ad = f(ADReal(a, da), ADReal(b, db), std::vector<ADReal>{ADReal(x, dx)});
  REL_TEST(f_ad, f(a, b, std::vector<Real>{x}), REL_TOL_CONSISTENCY);
 
  const Real ha = REL_PERTURBATION * a;
  const Real hb = REL_PERTURBATION * b;
  const Real hx = REL_PERTURBATION * x;
  REL_TEST(f_ad.derivatives()[0],
           (f(a + ha, b, std::vector<Real>{x}) - f(a - ha, b, std::vector<Real>{x})) / (2 * ha),
           REL_TOL_DERIVATIVE);
  REL_TEST(f_ad.derivatives()[1],
           (f(a, b + hb, std::vector<Real>{x}) - f(a, b - hb, std::vector<Real>{x})) / (2 * hb),
           REL_TOL_DERIVATIVE);
  REL_TEST(f_ad.derivatives()[2],
           (f(a, b, std::vector<Real>{x + hx}) - f(a, b, std::vector<Real>{x - hx})) / (2 * hx),
           REL_TOL_DERIVATIVE);
}
 
void
expectMixtureADDerivatives(const IdealGasMixtureFluidProperties & fp,
                           const Real p,
                           const Real T,
                           const Real x)
{
  const Real v = fp.v_from_p_T(p, T, {x});
  const Real e = fp.e_from_p_T(p, T, {x});
 
  expectADDerivatives([&fp](const auto & v, const auto & e, const auto & x)
                      { return fp.T_from_v_e(v, e, x); },
                      v,
                      e,
                      x);
  expectADDerivatives([&fp](const auto & v, const auto & e, const auto & x)
                      { return fp.p_from_v_e(v, e, x); },
                      v,
                      e,
                      x);
  expectADDerivatives([&fp](const auto & p, const auto & rho, const auto & x)
                      { return fp.e_from_p_rho(p, rho, x); },
                      p,
                      1 / v,
                      x);
  expectADDerivatives([&fp](const auto & p, const auto & T, const auto & x)
                      { return fp.c_from_p_T(p, T, x); },
                      p,
                      T,
                      x);
}
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, adDerivatives)
{
  expectMixtureADDerivatives(*_fp_mix, 100000, 400, 0.7);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, deprecatedSpecificHeatRatio)
{
  // The deprecated specific heat ratio without pressure and temperature arguments still works for
  // components with constant specific heats
  const std::vector<Real> x = {0.3, 0.7};
  REL_TEST(_fp_mix->mixtureSpecificHeatRatio(x),
           _fp_mix->mixtureSpecificHeatRatio(100000.0, 400.0, x),
           REL_TOL_CONSISTENCY);
 
  // It reports an error for components with temperature-dependent specific heats
  const std::vector<std::string> species = {"N2", "O2"};
  const std::string name = "air";
  const auto & air = buildNASA7Mixture(species, name);
  EXPECT_THROW(air.mixtureSpecificHeatRatio(std::vector<Real>{0.767, 0.233}), std::exception);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, nasa7Cantera)
{
  // Reformed fuel of a solid oxide cell, with mole fractions H2O 0.35, H2 0.30, CO 0.10, CO2 0.15,
  // CH4 0.05, and N2 0.05, at 1073.15 K (800 C). The secondary mass fractions and the reference
  // values were computed with Cantera 3.2.0 from gri30.yaml. Cantera uses the gas constant
  // 8.31446261815324 J/(mol K), whereas FluidProperties::_R is 8.3144598 J/(mol K). Properties
  // proportional to the gas constant are rescaled by the ratio of the two, density by its inverse,
  // and sound speed by its square root, so that only roundoff separates the results.
  const std::vector<std::string> species = {"H2O", "H2", "CO", "CO2", "CH4", "N2"};
  const std::string name = "reformate";
  const auto & fp = buildNASA7Mixture(species, name);
  const std::vector<Real> x = {0.03266496536638716,
                               0.1512807010437342,
                               0.3565358285737432,
                               0.04332374664128218,
                               0.0756511524284036};
  const Real p = 1.5e5;
  const Real T = 1073.15;
  const Real R_ratio = FluidProperties::_R / 8.31446261815324;
  const Real tol = 1e-12;
 
  const Real rho = 0.311262086619508 / R_ratio;
  const Real e = -7502142.694955811 * R_ratio;
  REL_TEST(fp.rho_from_p_T(p, T, x), rho, tol);
  REL_TEST(fp.e_from_p_T(p, T, x), e, tol);
  REL_TEST(fp.cp_from_p_T(p, T, x), 2216.7883592581443 * R_ratio, tol);
  REL_TEST(fp.cv_from_p_T(p, T, x), 1767.7281187454187 * R_ratio, tol);
  REL_TEST(fp.c_from_p_T(p, T, x), 777.386327867082 * std::sqrt(R_ratio), tol);
 
  // The (v, e) and (p, rho) variable sets describe the same state
  REL_TEST(fp.T_from_v_e(1 / rho, e, x), T, REL_TOL_CONSISTENCY);
  REL_TEST(fp.p_from_v_e(1 / rho, e, x), p, REL_TOL_CONSISTENCY);
  REL_TEST(fp.e_from_p_rho(p, rho, x), e, REL_TOL_CONSISTENCY);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, nasa7Inversion)
{
  // The internal energy and enthalpy of ideal gases do not depend on pressure or specific volume
  const Real p = 1e5;
  const Real v = 1;
 
  // Computes the temperature from the internal energy at temperatures spread geometrically between
  // T_low and T_high, 60 intervals placing consecutive temperatures about 7% apart
  const auto sweep = [p, v](const IdealGasMixtureFluidProperties & fp,
                            const std::vector<Real> & x,
                            const Real T_low,
                            const Real T_high)
  {
    const unsigned int n = 60;
    for (const auto i : make_range(n + 1))
    {
      const Real T = T_low * std::pow(T_high / T_low, Real(i) / n);
      REL_TEST(fp.T_from_v_e(v, fp.e_from_p_T(p, T, x), x), T, REL_TOL_CONSISTENCY);
    }
  };
 
  // From half the highest lower end to 1.5 times the lowest upper end of the temperature ranges of
  // the data of the components: 300 K for N2 and 3500 K for the others
  const std::vector<std::string> reformate_species = {"H2O", "H2", "CO", "CO2", "CH4", "N2"};
  const std::string reformate_name = "reformate";
  const auto & reformate = buildNASA7Mixture(reformate_species, reformate_name);
  sweep(reformate, {0.033, 0.151, 0.357, 0.043, 0.076}, 150, 5250);
 
  // Air, with the O2 mass fraction 0.233
  const std::vector<std::string> air_species = {"N2", "O2"};
  const std::string air_name = "air";
  const auto & air = buildNASA7Mixture(air_species, air_name);
  const std::vector<Real> x_air = {0.233};
  sweep(air, x_air, 150, 5250);
 
  // CH4 and N2 in equal mass fractions. Started below the solution, Newton's method can overshoot
  // to a spurious solution of the extrapolated CH4 polynomials.
  const std::vector<std::string> methane_nitrogen_species = {"N2", "CH4"};
  const std::string methane_nitrogen_name = "methane_nitrogen";
  const auto & methane_nitrogen =
      buildNASA7Mixture(methane_nitrogen_species, methane_nitrogen_name);
  sweep(methane_nitrogen, {0.5}, 150, 5250);
 
  // The internal energy of air, which includes the zero enthalpies of formation of N2 and O2,
  // passes through zero near 417 K, where a convergence test relative to the target energy would
  // demand accuracy below roundoff
  for (const Real e : {-1e-3, -1e-6, 0.0, 1e-6, 1e-3})
  {
    const Real T = air.T_from_v_e(v, e, x_air);
    ABS_TEST(
        air.e_from_p_T(p, T, x_air), e, REL_TOL_CONSISTENCY * air.cv_from_p_T(p, T, x_air) * T);
  }
 
  // A NASA 7-coefficient gas, CO2, as a secondary component of nitrogen with constant specific
  // heats; the latter has no temperature range
  {
    InputParameters params = _factory.getValidParams("NASA7IdealGasFluidProperties");
    params.set<MooseEnum>("species") = "CO2";
    _fe_problem->addUserObject("NASA7IdealGasFluidProperties", "co2", params);
  }
  {
    InputParameters params = _factory.getValidParams("IdealGasMixtureFluidProperties");
    params.set<std::vector<UserObjectName>>("component_fluid_properties") = {"fp_nitrogen", "co2"};
    _fe_problem->addUserObject("IdealGasMixtureFluidProperties", "nitrogen_co2", params);
  }
  sweep(
      _fe_problem->getUserObject<IdealGasMixtureFluidProperties>("nitrogen_co2"), {0.3}, 100, 5250);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, nasa7AdDerivatives)
{
  // Air, with the O2 mass fraction 0.233, on each side of the midpoint temperature of the data. The
  // temperature from (v, e) comes from a nonlinear Newton solve here, whose derivatives with
  // respect to the internal energy and the mass fraction the automatic differentiation must carry.
  const std::vector<std::string> species = {"N2", "O2"};
  const std::string name = "air";
  const auto & air = buildNASA7Mixture(species, name);
  for (const Real T : {700.0, 1500.0})
    expectMixtureADDerivatives(air, 100000, T, 0.233);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, wilkeCantera)
{
  // Mixture-averaged viscosity of N2, O2, H2, and CO2 with mole fractions 0.4, 0.1, 0.3, and 0.2
  // at 1000 K, computed with Cantera 3.2.0 from gri30.yaml, which applies the rule of Wilke to the
  // species viscosities below. Components with these viscosities and the molar masses of Cantera
  // reproduce it.
  const std::vector<std::string> species = {"N2", "O2", "H2", "CO2"};
  const std::vector<Real> molar_masses = {0.028014, 0.031998, 0.002016, 0.044009};
  const std::vector<Real> viscosities = {
      4.149814387750263e-05, 4.791169430728531e-05, 1.967636836557897e-05, 4.099132511194919e-05};
  const std::vector<Real> x = {0.1343776247270284, 0.025398958508315132, 0.3696371577355955};
  const Real mu_mix = 4.1726130812941113e-05;
 
  std::vector<UserObjectName> component_names;
  for (const auto i : index_range(species))
  {
    const std::string name = "wilke_" + species[i];
    InputParameters params = _factory.getValidParams("IdealGasFluidProperties");
    params.set<Real>("molar_mass") = molar_masses[i];
    params.set<Real>("mu") = viscosities[i];
    _fe_problem->addUserObject("IdealGasFluidProperties", name, params);
    component_names.push_back(name);
  }
  InputParameters params = _factory.getValidParams("IdealGasMixtureFluidProperties");
  params.set<std::vector<UserObjectName>>("component_fluid_properties") = component_names;
  params.set<MooseEnum>("transport_mixing_rule") = "wilke";
  _fe_problem->addUserObject("IdealGasMixtureFluidProperties", "wilke_mixture", params);
  const auto & fp = _fe_problem->getUserObject<IdealGasMixtureFluidProperties>("wilke_mixture");
 
  const Real p = 1e5;
  const Real T = 1000;
  REL_TEST(fp.mu_from_p_T(p, T, x), mu_mix, REL_TOL_SAVED_VALUE);
 
  // A single component keeps its own viscosity
  REL_TEST(fp.mu_from_p_T(p, T, {0, 0, 0}), viscosities[0], REL_TOL_SAVED_VALUE);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, masonSaxena)
{
  // The rule of Mason and Saxena for the thermal conductivity uses the interaction factors of the
  // rule of Wilke for the viscosity, which depend on the component viscosities and molar masses
  // only. Two consequences are checked: components whose thermal conductivities are proportional
  // to their viscosities give a mixture thermal conductivity proportional to the mixture viscosity
  // with the same factor, and components with equal viscosities and molar masses give the average
  // of their thermal conductivities weighted by molar fraction.
  const auto build_mixture = [this](const std::string & name,
                                    const std::vector<Real> & molar_masses,
                                    const std::vector<Real> & viscosities,
                                    const std::vector<Real> & conductivities)
  {
    std::vector<UserObjectName> component_names;
    for (const auto i : index_range(molar_masses))
    {
      const std::string component_name = name + "_" + std::to_string(i);
      InputParameters params = _factory.getValidParams("IdealGasFluidProperties");
      params.set<Real>("molar_mass") = molar_masses[i];
      params.set<Real>("mu") = viscosities[i];
      params.set<Real>("k") = conductivities[i];
      _fe_problem->addUserObject("IdealGasFluidProperties", component_name, params);
      component_names.push_back(component_name);
    }
    InputParameters params = _factory.getValidParams("IdealGasMixtureFluidProperties");
    params.set<std::vector<UserObjectName>>("component_fluid_properties") = component_names;
    params.set<MooseEnum>("transport_mixing_rule") = "wilke";
    _fe_problem->addUserObject("IdealGasMixtureFluidProperties", name, params);
    return &_fe_problem->getUserObject<IdealGasMixtureFluidProperties>(name);
  };
 
  const Real p = 1e5;
  const Real T = 1000;
 
  // Molar masses and viscosities of N2, O2, H2, and CO2 near 1000 K, and an arbitrary factor
  const std::vector<Real> molar_masses = {0.028014, 0.031998, 0.002016, 0.044009};
  const std::vector<Real> viscosities = {4.15e-05, 4.79e-05, 1.97e-05, 4.10e-05};
  const Real factor = 1500;
  std::vector<Real> conductivities;
  for (const auto mu : viscosities)
    conductivities.push_back(factor * mu);
  const auto * proportional =
      build_mixture("proportional", molar_masses, viscosities, conductivities);
  const std::vector<Real> x = {0.13, 0.03, 0.37};
  REL_TEST(proportional->k_from_p_T(p, T, x),
           factor * proportional->mu_from_p_T(p, T, x),
           REL_TOL_CONSISTENCY);
 
  const auto * equal = build_mixture("equal", {0.028, 0.028}, {4e-5, 4e-5}, {0.06, 0.08});
  REL_TEST(equal->k_from_p_T(p, T, {0.25}), 0.75 * 0.06 + 0.25 * 0.08, REL_TOL_CONSISTENCY);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, wilkeAdDerivatives)
{
  // Air, with the O2 mass fraction 0.233, and a fuel of H2O and H2, with the H2 mass fraction
  // 0.05, near the temperature of a solid oxide cell. The viscosity of H2O comes from the IAPWS
  // formulation and the others from kinetic theory.
  const std::vector<std::pair<std::vector<std::string>, Real>> mixtures = {{{"N2", "O2"}, 0.233},
                                                                           {{"H2O", "H2"}, 0.05}};
  for (const auto & [species, x] : mixtures)
  {
    const std::string name = "wilke_" + species[0] + "_" + species[1];
    const auto & fp = buildNASA7Mixture(species, name, "wilke");
    expectADDerivatives([&fp](const auto & p, const auto & T, const auto & x)
                        { return fp.mu_from_p_T(p, T, x); },
                        1e5,
                        1073.15,
                        x);
    expectADDerivatives([&fp](const auto & p, const auto & T, const auto & x)
                        { return fp.k_from_p_T(p, T, x); },
                        1e5,
                        1073.15,
                        x);
  }
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, binaryDiffusion)
{
  // Binary diffusion coefficients at 1e5 Pa and 1000 K computed with Cantera 3.2.0 from
  // gri30.yaml, which the Chapman-Enskog coefficients reproduce within 0.1% for pairs without H2O
  const auto & fp = buildNASA7Mixture({"N2", "O2", "H2"}, "diffusion");
  const Real p = 1e5;
  const Real T = 1000;
  const Real tol = 2e-3;
  REL_TEST(fp.binaryDiffusionCoefficient(p, T, 0, 1), 0.00016508453019480936, tol);
  REL_TEST(fp.binaryDiffusionCoefficient(p, T, 2, 0), 0.0005928256767875003, tol);
 
  // The coefficient is symmetric in the pair
  REL_TEST(fp.binaryDiffusionCoefficient(p, T, 0, 2),
           fp.binaryDiffusionCoefficient(p, T, 2, 0),
           REL_TOL_CONSISTENCY);
 
  // The derivatives with respect to pressure (index 0) and temperature (index 1): the coefficient
  // is inversely proportional to pressure
  DNDerivativeType dp_dp, dT_dT;
  Moose::derivInsert(dp_dp, 0, 1);
  Moose::derivInsert(dT_dT, 1, 1);
  const ADReal D_ad = fp.binaryDiffusionCoefficient(ADReal(p, dp_dp), ADReal(T, dT_dT), 0, 2);
  const Real D = fp.binaryDiffusionCoefficient(p, T, 0, 2);
  const Real h = REL_PERTURBATION * T;
  REL_TEST(D_ad, D, REL_TOL_CONSISTENCY);
  REL_TEST(D_ad.derivatives()[0], -D / p, REL_TOL_CONSISTENCY);
  REL_TEST(D_ad.derivatives()[1],
           (fp.binaryDiffusionCoefficient(p, T + h, 0, 2) -
            fp.binaryDiffusionCoefficient(p, T - h, 0, 2)) /
               (2 * h),
           REL_TOL_DERIVATIVE);
 
  // IdealGasFluidProperties components have no Lennard-Jones parameters
  EXPECT_THROW(_fp_mix->binaryDiffusionCoefficient(p, T, 0, 1), std::exception);
}
 
TEST_F(IdealGasMixtureFluidPropertiesTest, entropySingleComponent)
{
  // A mixture with only one component present has the entropy of that component
  // Derivative index 0 is pressure and index 1 is temperature
  const Real p = 1e5;
  const Real T = 400;
  DNDerivativeType dp_dp, dT_dT;
  Moose::derivInsert(dp_dp, 0, 1);
  Moose::derivInsert(dT_dT, 1, 1);
  const ADReal p_ad(p, dp_dp);
  const ADReal T_ad(T, dT_dT);
  for (const auto & [x, fp] : {std::pair<Real, const IdealGasFluidProperties *>(0, _fp_steam),
                               std::pair<Real, const IdealGasFluidProperties *>(1, _fp_nitrogen)})
  {
    Real s, ds_dp, ds_dT;
    fp->s_from_p_T(p, T, s, ds_dp, ds_dT);
    REL_TEST(_fp_mix->s_from_p_T(p, T, {x}), s, REL_TOL_CONSISTENCY);
    const ADReal s_ad = _fp_mix->s_from_p_T(p_ad, T_ad, std::vector<ADReal>{ADReal(x)});
    REL_TEST(s_ad.derivatives()[0], ds_dp, REL_TOL_CONSISTENCY);
    REL_TEST(s_ad.derivatives()[1], ds_dT, REL_TOL_CONSISTENCY);
  }
}