diff --git a/Documentation/ProjectFile/properties/property/t_IdealGasLawBinaryMixture.md b/Documentation/ProjectFile/properties/property/t_IdealGasLawBinaryMixture.md
new file mode 100644
index 0000000000000000000000000000000000000000..a66be060fca297dbcfc79b3840760efb14d8669f
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/t_IdealGasLawBinaryMixture.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::IdealGasLawBinaryMixture
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 7ca90c0f3600c0bfe4a64eba69fc3306a3081821..df6d8372b08e948f970557541f52234d57a9ede5 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -93,6 +93,11 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
return createIdealGasLaw(config);
}
+ if (boost::iequals(property_type, "IdealGasLawBinaryMixture"))
+ {
+ return createIdealGasLawBinaryMixture(config);
+ }
+
if (boost::iequals(property_type, "StrainDependentPermeability"))
{
return createStrainDependentPermeability(
diff --git a/MaterialLib/MPL/Properties/CreateIdealGasLawBinaryMixture.cpp b/MaterialLib/MPL/Properties/CreateIdealGasLawBinaryMixture.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..14349de233bca5c0dce3b915ca975241258f12cc
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateIdealGasLawBinaryMixture.cpp
@@ -0,0 +1,34 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date Jan, 2021
+ *
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "IdealGasLawBinaryMixture.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<IdealGasLawBinaryMixture> createIdealGasLawBinaryMixture(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "IdealGasLawBinaryMixture");
+
+ // Second access for storage.
+ //! \ogs_file_param{properties__property__name}
+ auto property_name = config.peekConfigParameter<std::string>("name");
+
+ DBUG("Create IdealGasLawBinaryMixture medium property {:s}.",
+ property_name);
+
+ //! \ogs_file_param_special{properties__property__IdealGasLawBinaryMixture}
+ return std::make_unique<IdealGasLawBinaryMixture>(std::move(property_name));
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateIdealGasLawBinaryMixture.h b/MaterialLib/MPL/Properties/CreateIdealGasLawBinaryMixture.h
new file mode 100644
index 0000000000000000000000000000000000000000..2faea27d056c6d99e0e3860431c36d92f7f637a1
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateIdealGasLawBinaryMixture.h
@@ -0,0 +1,31 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date Jan, 2021
+ *
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class IdealGasLawBinaryMixture;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<IdealGasLawBinaryMixture> createIdealGasLawBinaryMixture(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 71d1974f9cf125b138cbef86fab3e354f8f2826a..fbf83202de93fd6a05292f607af916bc505ff687 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -30,6 +30,7 @@
#include "CreateFunction.h"
#include "CreateGasPressureDependentPermeability.h"
#include "CreateIdealGasLaw.h"
+#include "CreateIdealGasLawBinaryMixture.h"
#include "CreateKozenyCarmanModel.h"
#include "CreateLinear.h"
#include "CreateOrthotropicEmbeddedFracturePermeability.h"
diff --git a/MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.cpp b/MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..00983323508a88d9812f3d70a6af64734bfc1daf
--- /dev/null
+++ b/MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.cpp
@@ -0,0 +1,146 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date Jan, 2021
+ * \brief
+ *
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.h"
+
+#include "MaterialLib/PhysicalConstant.h"
+
+namespace MaterialPropertyLib
+{
+IdealGasLawBinaryMixture::IdealGasLawBinaryMixture(std::string name)
+{
+ name_ = std::move(name);
+}
+
+// Enums with indices for the sorting of the partial density and its derivatives
+// in the result vector or the result tensor of the functions value and dValue.
+enum class VariableIndex : int
+{
+ gas_pressure,
+ capillary_pressure,
+ temperature
+};
+enum class DensityIndex : int
+{
+ gas_phase_density,
+ vapour_density,
+ air_density
+};
+
+PropertyDataType IdealGasLawBinaryMixture::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ const double R = MaterialLib::PhysicalConstant::IdealGasConstant;
+ const double pGR = std::get<double>(
+ variable_array[static_cast<int>(Variable::phase_pressure)]);
+ const double T = std::get<double>(
+ variable_array[static_cast<int>(Variable::temperature)]);
+ const double MW = std::get<double>(
+ variable_array[static_cast<int>(Variable::molar_mass_vapour)]);
+ const double MC = std::get<double>(
+ variable_array[static_cast<int>(Variable::molar_mass)]);
+ // mole fraction of the vapour component
+ const double xnWG = std::get<double>(
+ variable_array[static_cast<int>(Variable::molar_fraction)]);
+ const double xnCG = 1. - xnWG;
+
+ const double rhoWGR = xnWG * pGR * MW / R / T;
+ const double rhoCGR = xnCG * pGR * MC / R / T;
+ Eigen::Matrix<double, 3, 1> result;
+
+ result[static_cast<int>(DensityIndex::vapour_density)] = rhoWGR;
+ result[static_cast<int>(DensityIndex::air_density)] = rhoCGR;
+ result[static_cast<int>(DensityIndex::gas_phase_density)] = rhoWGR + rhoCGR;
+
+ return result;
+}
+
+PropertyDataType IdealGasLawBinaryMixture::dValue(
+ VariableArray const& variable_array, Variable const /*primary_variable*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ const double R = MaterialLib::PhysicalConstant::IdealGasConstant;
+ const double pGR = std::get<double>(
+ variable_array[static_cast<int>(Variable::phase_pressure)]);
+ const double T = std::get<double>(
+ variable_array[static_cast<int>(Variable::temperature)]);
+ const double MW = std::get<double>(
+ variable_array[static_cast<int>(Variable::molar_mass_vapour)]);
+ const double MC = std::get<double>(
+ variable_array[static_cast<int>(Variable::molar_mass)]);
+ // mole fraction of the vapour component
+ const double xnWG = std::get<double>(
+ variable_array[static_cast<int>(Variable::molar_fraction)]);
+ const double xnCG = 1. - xnWG;
+ const double dxnWG_dpGR =
+ std::get<double>(variable_array[static_cast<int>(Variable::dxn_dpGR)]);
+ const double dxnWG_dpCap =
+ std::get<double>(variable_array[static_cast<int>(Variable::dxn_dpCap)]);
+ const double dxnWG_dT =
+ std::get<double>(variable_array[static_cast<int>(Variable::dxn_dT)]);
+
+ Eigen::Matrix<double, 3, 3> result;
+
+ const double drhoCGR_dpGR = -MC / R / T * (pGR * dxnWG_dpGR - xnCG);
+ const double drhoCGR_dpCap = -pGR * MC / R / T * dxnWG_dpCap;
+ const double drhoCGR_dT = -pGR * MC / R / T * (dxnWG_dT + xnCG / T);
+
+ const double drhoWGR_dpGR = MW / R / T * (pGR * dxnWG_dpGR + xnWG);
+ const double drhoWGR_dpCap = pGR * MW / R / T * dxnWG_dpCap;
+ const double drhoWGR_dT = pGR * MW / R / T / T * (T * dxnWG_dT - xnWG);
+
+ const double drhoGR_dpGR = drhoCGR_dpGR + drhoWGR_dpGR;
+ const double drhoGR_dpCap = drhoCGR_dpCap + drhoWGR_dpCap;
+ const double drhoGR_dT = drhoCGR_dT + drhoWGR_dT;
+
+ result(static_cast<int>(DensityIndex::air_density),
+ static_cast<int>(VariableIndex::gas_pressure)) = drhoCGR_dpGR;
+ result(static_cast<int>(DensityIndex::air_density),
+ static_cast<int>(VariableIndex::capillary_pressure)) = drhoCGR_dpCap;
+ result(static_cast<int>(DensityIndex::air_density),
+ static_cast<int>(VariableIndex::temperature)) = drhoCGR_dT;
+
+ result(static_cast<int>(DensityIndex::vapour_density),
+ static_cast<int>(VariableIndex::gas_pressure)) = drhoWGR_dpGR;
+ result(static_cast<int>(DensityIndex::vapour_density),
+ static_cast<int>(VariableIndex::capillary_pressure)) = drhoWGR_dpCap;
+ result(static_cast<int>(DensityIndex::vapour_density),
+ static_cast<int>(VariableIndex::temperature)) = drhoWGR_dT;
+
+ result(static_cast<int>(DensityIndex::gas_phase_density),
+ static_cast<int>(VariableIndex::gas_pressure)) = drhoGR_dpGR;
+ result(static_cast<int>(DensityIndex::gas_phase_density),
+ static_cast<int>(VariableIndex::capillary_pressure)) = drhoGR_dpCap;
+ result(static_cast<int>(DensityIndex::gas_phase_density),
+ static_cast<int>(VariableIndex::temperature)) = drhoGR_dT;
+
+ return result;
+}
+
+PropertyDataType IdealGasLawBinaryMixture::d2Value(
+ VariableArray const& /*variable_array*/,
+ Variable const /*primary_variable1*/, Variable const /*primary_variable2*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ OGS_FATAL("IdealGasLawBinaryMixture::d2Value is not implemented.");
+
+ return 0.;
+}
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.h b/MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.h
new file mode 100644
index 0000000000000000000000000000000000000000..220729724e7ee214df4cc4c0a49b7e936058c8bd
--- /dev/null
+++ b/MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.h
@@ -0,0 +1,62 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date Jan, 2021
+ * \brief
+ *
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+#pragma once
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Phase;
+/**
+ * \brief Density function for binary ideal gases.
+ * \details This property must be a phase property, it
+ * computes the density of a mixture of ideal gases as function of phase
+ * pressure, capillary pressure, and temperature.
+ */
+class IdealGasLawBinaryMixture final : public Property
+{
+public:
+ explicit IdealGasLawBinaryMixture(std::string name);
+
+ void checkScale() const override
+ {
+ if (!std::holds_alternative<Phase*>(scale_))
+ {
+ OGS_FATAL(
+ "The property 'IdealGasLawBinaryMixture' is implemented on the "
+ "'phase' scale only.");
+ }
+ }
+
+ /// Those methods override the base class implementations and
+ /// actually compute and set the property values_ and dValues_.
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+ PropertyDataType d2Value(VariableArray const& variable_array,
+ Variable const variable1, Variable const variable2,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+};
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index 52e41f7bf25b91020964e084a6987fcaa719bf6f..da771ea6d3116e3facdc81dfcc67b897e5913495 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -31,6 +31,7 @@
#include "Function.h"
#include "GasPressureDependentPermeability.h"
#include "IdealGasLaw.h"
+#include "IdealGasLawBinaryMixture.h"
#include "Linear.h"
#include "OrthotropicEmbeddedFracturePermeability.h"
#include "Parameter.h"
diff --git a/MaterialLib/MPL/VariableType.h b/MaterialLib/MPL/VariableType.h
index c1288ab742adb6db5105ee95d45511487741d1e7..316196d2251f7dfd42bf2c90f347a808a8ff266b 100644
--- a/MaterialLib/MPL/VariableType.h
+++ b/MaterialLib/MPL/VariableType.h
@@ -47,6 +47,9 @@ enum class Variable : int
concentration,
density,
displacement,
+ dxn_dT, // xn - molar fraction
+ dxn_dpGR, // derivatives w.r.t. temperature, gas pressure, and
+ dxn_dpCap, // capillary pressure
effective_pore_pressure,
enthalpy_of_evaporation,
equivalent_plastic_strain,
@@ -55,6 +58,7 @@ enum class Variable : int
liquid_saturation,
mechanical_strain,
molar_mass,
+ molar_mass_vapour,
molar_fraction,
phase_pressure,
porosity,
@@ -74,6 +78,9 @@ static const std::array<std::string,
variable_enum_to_string{{"capillary_pressure",
"concentration",
"density",
+ "dxn_dT",
+ "dxn_dpGR",
+ "dxn_dpCap",
"displacement",
"effective_pore_pressure",
"enthalpy_of_evaporation",
@@ -83,6 +90,7 @@ static const std::array<std::string,
"liquid_saturation",
"mechanical_strain",
"molar_mass",
+ "molar_mass_vapour",
"molar_fraction",
"phase_pressure",
"porosity",
diff --git a/Tests/MaterialLib/TestMPLIdealGasLawBinaryMixture.cpp b/Tests/MaterialLib/TestMPLIdealGasLawBinaryMixture.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b8c1640baef4d08410371d4fab8784f55d9ac7ab
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLIdealGasLawBinaryMixture.cpp
@@ -0,0 +1,268 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+#include <gtest/gtest.h>
+
+#include "MaterialLib/MPL/Properties/IdealGasLawBinaryMixture.h"
+#include "MaterialLib/PhysicalConstant.h"
+#include "TestMPL.h"
+#include "Tests/TestTools.h"
+
+double molarFractionFantasy(const double pGR, const double pGR_min,
+ const double pGR_max, const double pCap,
+ const double pCap_min, const double pCap_max,
+ const double T, const double T_min,
+ const double T_max)
+{
+ auto const dpGR = pGR_max - pGR_min;
+ auto const dpCap = pCap_max - pCap_min;
+ auto const dT = T_max - T_min;
+
+ auto const a = 1. / (6. * std::pow(dpGR, 2.));
+ auto const b = 1. / (6. * dpGR);
+ auto const c = 1. / (6. * std::pow(dpCap, 2.));
+ auto const d = 1. / (6. * (dpCap));
+ auto const e = 1. / (6. * std::pow(dT, 2.));
+ auto const f = 1. / (6. * (dT));
+
+ return a * (pGR - pGR_min) * (pGR - pGR_min) + b * (pGR - pGR_min) +
+ c * (pCap - pCap_min) * (pCap - pCap_min) + d * (pCap - pCap_min) +
+ e * (T - T_min) * (T - T_min) + f * (T - T_min);
+}
+
+std::array<double, 3> molarMassFantasyDerivatives(
+ const double pGR, const double pGR_min, const double pGR_max,
+ const double pCap, const double pCap_min, const double pCap_max,
+ const double T, const double T_min, const double T_max)
+{
+ auto const dpGR = pGR_max - pGR_min;
+ auto const dpCap = pCap_max - pCap_min;
+ auto const dT = T_max - T_min;
+
+ auto const a = 1. / (6. * std::pow(dpGR, 2.));
+ auto const b = 1. / (6. * dpGR);
+ auto const c = 1. / (6. * std::pow(dpCap, 2.));
+ auto const d = 1. / (6. * (dpCap));
+ auto const e = 1. / (6. * std::pow(dT, 2.));
+ auto const f = 1. / (6. * (dT));
+
+ const double d_dpGR = 2 * a * (pGR - pGR_min) + b;
+ const double d_dpCap = 2 * c * (pCap - pCap_min) + d;
+ const double d_dT = 2 * e * (T - T_min) + f;
+
+ return {d_dpGR, d_dpCap, d_dT};
+}
+
+TEST(MaterialPropertyLib, IdealGasLawBinaryMixture)
+{
+ const double MC = 0.028949;
+ const double MW = 0.018052;
+
+ auto const density_model = MPL::IdealGasLawBinaryMixture("density");
+
+ enum
+ {
+ gas_pressure,
+ capillary_pressure,
+ temperature
+ };
+ enum
+ {
+ gas_phase_density,
+ vapour_density,
+ air_density
+ };
+
+ ParameterLib::SpatialPosition const pos;
+ double const t = std::numeric_limits<double>::quiet_NaN();
+ double const dt = std::numeric_limits<double>::quiet_NaN();
+ MPL::VariableArray vars;
+
+ vars[static_cast<int>(MPL::Variable::molar_mass_vapour)] = MW;
+ vars[static_cast<int>(MPL::Variable::molar_mass)] = MC;
+
+ // Density derivatives of gases are dependent on pressure and temperature.
+ // In the case of binary gases, there is also the dependence of the
+ // composition (e.g. molar fractions).
+
+ // The molar fraction of the vapour component in the gas phase is
+ // represented here by a fantasy function. In Reality, this is actually a
+ // physical quantity that depends on the temperature-dependent vapour
+ // pressure, the capillary pressure-dependent correction factor for curved
+ // menisci and the pressure of the gas phase.
+
+ // In this test, however, the function is merely a polynomial,
+ // which over the entire parameter space pGR-pCap-T assumes
+ // values in the range [0,1].
+
+ // Boundaries of the parameter space
+ auto const pGR_min = 80000.;
+ auto const pGR_max = 320000.;
+ auto const pCap_min = 100.;
+ auto const pCap_max = 1.e8;
+ auto const T_min = 270.;
+ auto const T_max = 520.;
+
+ // travel through the entire parameter space
+ for (double pGR = pGR_min; pGR <= pGR_max; pGR += 200.)
+ {
+ for (double pCap = pCap_min; pCap <= pCap_max; pCap *= 2.)
+ {
+ for (double T = T_min; T <= T_max; T += 0.2)
+ {
+ vars[static_cast<int>(MPL::Variable::phase_pressure)] = pGR;
+ vars[static_cast<int>(MPL::Variable::temperature)] = T;
+
+ // Molar fraction fantasy function
+ auto const xnWG =
+ molarFractionFantasy(pGR, pGR_min, pGR_max, pCap, pCap_min,
+ pCap_max, T, T_min, T_max);
+
+ // ... and its derivatives
+ auto const dxnWG = molarMassFantasyDerivatives(
+ pGR, pGR_min, pGR_max, pCap, pCap_min, pCap_max, T, T_min,
+ T_max);
+
+ vars[static_cast<int>(MPL::Variable::molar_fraction)] = xnWG;
+ vars[static_cast<int>(MPL::Variable::dxn_dpGR)] =
+ dxnWG[gas_pressure];
+ vars[static_cast<int>(MPL::Variable::dxn_dpCap)] =
+ dxnWG[capillary_pressure];
+ vars[static_cast<int>(MPL::Variable::dxn_dT)] =
+ dxnWG[temperature];
+
+ auto const density = std::get<Eigen::Matrix<double, 3, 1>>(
+ density_model.value(vars, pos, t, dt));
+
+ auto const rhoGR = density[gas_phase_density];
+ auto const rhoCGR = density[air_density];
+ auto const rhoWGR = density[vapour_density];
+
+ const double R =
+ MaterialLib::PhysicalConstant::IdealGasConstant;
+
+ auto const rhoCGR_ = (1. - xnWG) * pGR * MC / R / T;
+ auto const rhoWGR_ = xnWG * pGR * MW / R / T;
+ auto const rhoGR_ = rhoCGR_ + rhoWGR_;
+
+ ASSERT_NEAR(rhoGR, rhoGR_, 1.e-18);
+ ASSERT_NEAR(rhoCGR, rhoCGR_, 1.e-18);
+ ASSERT_NEAR(rhoWGR, rhoWGR_, 1.e-18);
+
+ auto const density_derivative =
+ std::get<Eigen::Matrix<double, 3, 3>>(density_model.dValue(
+ vars, MPL::Variable::phase_pressure, pos, t, dt));
+
+ // Check derivatives via central differences
+ // Pertubation of gas pressure
+ auto const eps_pGR = 1.;
+
+ auto const pGR_plus = pGR + eps_pGR;
+ auto const xnWG_plus =
+ molarFractionFantasy(pGR_plus, pGR_min, pGR_max, pCap,
+ pCap_min, pCap_max, T, T_min, T_max);
+
+ vars[static_cast<int>(MPL::Variable::phase_pressure)] =
+ pGR_plus;
+ vars[static_cast<int>(MPL::Variable::molar_fraction)] =
+ xnWG_plus;
+
+ auto density_plus = std::get<Eigen::Matrix<double, 3, 1>>(
+ density_model.value(vars, pos, t, dt));
+
+ auto const pGR_minus = pGR - eps_pGR;
+ vars[static_cast<int>(MPL::Variable::phase_pressure)] =
+ pGR_minus;
+ vars[static_cast<int>(MPL::Variable::molar_fraction)] =
+ molarFractionFantasy(pGR_minus, pGR_min, pGR_max, pCap,
+ pCap_min, pCap_max, T, T_min, T_max);
+
+ auto density_minus = std::get<Eigen::Matrix<double, 3, 1>>(
+ density_model.value(vars, pos, t, dt));
+
+ for (std::size_t d = gas_phase_density; d <= air_density; d++)
+ {
+ auto const drho_dpGR =
+ (density_plus[d] - density_minus[d]) / (2 * eps_pGR);
+
+ ASSERT_NEAR(drho_dpGR, density_derivative(d, gas_pressure),
+ 1.e-10);
+ }
+
+ // Pertubation of capillary pressure
+ auto const eps_pCap = 1.;
+ auto const pCap_plus = pCap + eps_pCap;
+
+ vars[static_cast<int>(MPL::Variable::phase_pressure)] = pGR;
+ vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ pCap_plus;
+ vars[static_cast<int>(MPL::Variable::molar_fraction)] =
+ molarFractionFantasy(pGR, pGR_min, pGR_max, pCap_plus,
+ pCap_min, pCap_max, T, T_min, T_max);
+
+ density_plus = std::get<Eigen::Matrix<double, 3, 1>>(
+ density_model.value(vars, pos, t, dt));
+
+ auto const pCap_minus = pCap - eps_pCap;
+ vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ pCap_minus;
+ vars[static_cast<int>(MPL::Variable::molar_fraction)] =
+ molarFractionFantasy(pGR, pGR_min, pGR_max, pCap_minus,
+ pCap_min, pCap_max, T, T_min, T_max);
+
+ density_minus = std::get<Eigen::Matrix<double, 3, 1>>(
+ density_model.value(vars, pos, t, dt));
+
+ for (std::size_t d = gas_phase_density; d <= air_density; d++)
+ {
+ auto const drho_dpCap =
+ (density_plus[d] - density_minus[d]) / (2 * eps_pCap);
+
+ ASSERT_NEAR(drho_dpCap,
+ density_derivative(d, capillary_pressure),
+ 1.e-10);
+ }
+
+ // Pertubation of temperature
+ auto const eps_T = 0.01;
+ auto const T_plus = T + eps_T;
+
+ vars[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ pCap;
+ vars[static_cast<int>(MPL::Variable::temperature)] = T_plus;
+ vars[static_cast<int>(MPL::Variable::molar_fraction)] =
+ molarFractionFantasy(pGR, pGR_min, pGR_max, pCap, pCap_min,
+ pCap_max, T_plus, T_min, T_max);
+
+ density_plus = std::get<Eigen::Matrix<double, 3, 1>>(
+ density_model.value(vars, pos, t, dt));
+
+ auto const T_minus = T - eps_T;
+ vars[static_cast<int>(MPL::Variable::temperature)] = T_minus;
+ vars[static_cast<int>(MPL::Variable::molar_fraction)] =
+ molarFractionFantasy(pGR, pGR_min, pGR_max, pCap, pCap_min,
+ pCap_max, T_minus, T_min, T_max);
+
+ density_minus = std::get<Eigen::Matrix<double, 3, 1>>(
+ density_model.value(vars, pos, t, dt));
+
+ for (std::size_t d = gas_phase_density; d <= air_density; d++)
+ {
+ auto const drho_dT =
+ (density_plus[d] - density_minus[d]) / (2 * eps_T);
+
+ ASSERT_NEAR(drho_dT, density_derivative(d, temperature),
+ 1.e-10);
+ }
+
+ } // end of T-loop
+ } // end of pCap-loop
+ } // end of pGR-loop
+}