[Mat] Added classes for fluid density model

MaterialLib/CMakeLists.txt

-GET_SOURCE_FILES(SOURCES)
-APPEND_SOURCE_FILES(SOURCES Adsorption)
-APPEND_SOURCE_FILES(SOURCES SolidModels)
+# Source files
+get_source_files(SOURCES)
+append_source_files(SOURCES Adsorption)
+append_source_files(SOURCES SolidModels)
 
-add_library(MaterialLib ${SOURCES})
+append_source_files(SOURCES Fluid)
+append_source_files(SOURCES Fluid/Density)
+
+add_library(MaterialLib ${SOURCES} )
 target_link_libraries(MaterialLib
     BaseLib
 )

MaterialLib/Fluid/ConstantFluidProperty.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   ConstantFluidProperty.h
+ *
+ * Created on August 15, 2016, 12:11 PM
+ */
+
+#ifndef CONSTANTFLUIDPROPERTY_H
+#define CONSTANTFLUIDPROPERTY_H
+
+#include "FluidProperty.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/// Constant fluid density properties
+class ConstantFluidProperty : public FluidProperty
+{
+public:
+    ConstantFluidProperty(const double value) : FluidProperty(), _value(value)
+    {
+    }
+
+    /// Get model name.
+    virtual std::string getName() const final { return "constant"; }
+    /// Get property value.
+    /// \param var_vals Variable values in an array. The order of its elements
+    ///                 is given in enum class PropertyVariable.
+    virtual double getValue(const double /* var_vals*/[]) const final
+    {
+        return _value;
+    }
+
+    /// Get the partial differential of the property value
+    /// \param var_vals  Variable values  in an array. The order of its elements
+    ///                   is given in enum class PropertyVariable.
+    virtual double getdValue(const double /* var_vals*/[],
+                             const PropertyVariable /* var */) const final
+    {
+        return 0.;
+    }
+
+private:
+    double _value;
+};
+
+}  // end namespace
+}  // end namespace
+#endif /* CONSTANTFLUIDPROPERTY_H */

MaterialLib/Fluid/Density/IdealGasLaw.h

+/*!
+   \file  IdealGasLow.h
+   \brief Declaration of class LIdealGasLow for fluid density by the ideal gas
+   law
+          depending on one variable linearly.
+
+   \author Wenqing Wang
+   \date Jan 2015
+
+   \copyright
+    Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+               Distributed under a Modified BSD License.
+               See accompanying file LICENSE.txt or
+               http://www.opengeosys.org/project/license
+ */
+#ifndef IDEAL_GAS_LAW_H_
+#define IDEAL_GAS_LAW_H_
+
+#include <string>
+#include <cassert>
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PhysicalConstant.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/// Fluid density by ideal gas low
+
+class IdealGasLaw : public FluidProperty
+{
+public:
+    ///   \param molar_mass Molar mass of the gas phase.
+    IdealGasLaw(const double molar_mass)
+        : FluidProperty(),
+          _molar_mass(molar_mass),
+          _derivative_functions{&IdealGasLaw::dIdealGasLaw_dT,
+                                &IdealGasLaw::dIdealGasLaw_dp}
+    {
+    }
+
+    /// Get density model name.
+    virtual std::string getName() const final { return "Ideal gas law"; }
+    /// Get density value.
+    /// \param var_vals Variable values in an array. The order of its elements
+    ///                 is given in enum class PropertyVariable.
+    virtual double getValue(const double var_vals[]) const final
+    {
+        return _molar_mass * var_vals[static_cast<int>(PropertyVariable::pg)] /
+               (PhysicalConstant::IdealGasConstant *
+                var_vals[static_cast<int>(PropertyVariable::T)]);
+    }
+
+    /// Get the partial differential of the density with respect to temperature
+    /// or gas pressure.
+    /// \param var_vals  Variable values  in an array. The order of its elements
+    ///                   is given in enum class PropertyVariable.
+    virtual double getdValue(const double var_vals[],
+                             const PropertyVariable var) const final
+    {
+        assert(var == PropertyVariable::T || var == PropertyVariable::pg);
+
+        const int func_id = (var == PropertyVariable::T) ? 0 : 1;
+        return (this->*_derivative_functions[func_id])(
+            var_vals[static_cast<int>(PropertyVariable::T)],
+            var_vals[static_cast<int>(PropertyVariable::pg)]);
+    }
+
+private:
+    /// Molar mass of gas phase.
+    double _molar_mass;
+
+    /// Get the partial differential of density with the respect to temperature
+    /// \param T  Temperature in K.
+    /// \param pg Gas phase pressure in Pa.
+
+    double dIdealGasLaw_dT(const double T, const double pg) const
+    {
+        return -_molar_mass * pg / (PhysicalConstant::IdealGasConstant * T * T);
+    }
+
+    /// Get the partial differential of density with the respect to pressure
+    /// \param T  Temperature in K.
+    /// \param pg Gas phase pressure in Pa.
+    double dIdealGasLaw_dp(const double T, const double /* pg */) const
+    {
+        return _molar_mass / (PhysicalConstant::IdealGasConstant * T);
+    }
+
+    typedef double (IdealGasLaw::*DerivativeFunctionPointer)(const double,
+                                                     const double) const;
+
+    /// An array of pointer to derivative functions.
+    DerivativeFunctionPointer _derivative_functions[2];
+};
+}  // end namespace
+}  // end namespace
+#endif

MaterialLib/Fluid/Density/LinearTemperatureDependentDensity.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   LinearTemperatureDependentDensity.h
+ * \author: wenqing
+ *
+ * Created on August 10, 2016, 11:34 AM
+ */
+
+#ifndef LINEARTEMPERATUREDEPENDENTDENSITY_H
+#define LINEARTEMPERATUREDEPENDENTDENSITY_H
+
+#include <string>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/// Linear temperature dependent density model.
+class LinearTemperatureDependentDensity : public FluidProperty
+{
+public:
+    /// \param config  ConfigTree object which contains the input data
+    ///                including <type>temperature_dependent</type> and it has
+    ///                a tag of <density>
+    LinearTemperatureDependentDensity(BaseLib::ConfigTree const& config)
+        :  //! \ogs_file_param{material__fluid__density__linear_temperature__rho0}
+    _rho0(config.getConfigParameter<double>("rho0")),
+          //! \ogs_file_param{material__fluid__density__linear_temperature__temperature0}
+    _temperature0(config.getConfigParameter<double>("temperature0")),
+          //! \ogs_file_param{material__fluid__density__linear_temperature__beta}
+    _beta(config.getConfigParameter<double>("beta"))
+    {
+    }
+
+    /// Get model name.
+    virtual std::string getName() const final
+    {
+        return "Linear temperature dependent density";
+    }
+
+    /// Get density value.
+    /// \param var_vals Variable values in an array. The order of its elements
+    ///                 is given in enum class PropertyVariable.
+    virtual double getValue(const double var_vals[]) const final
+    {
+        return _rho0 *
+               (1 +
+                _beta * (var_vals[static_cast<int>(PropertyVariable::T)] -
+                         _temperature0));
+    }
+
+    /// Get the partial differential of the density with respect to temperature.
+    /// \param var_vals  Variable values  in an array. The order of its elements
+    ///                   is given in enum class PropertyVariable.
+    virtual double getdValue(const double /* var_vals */[],
+                             const PropertyVariable /* var  */) const final
+    {
+        return _rho0 * _beta;
+        ;
+    }
+
+private:
+    double _rho0;          ///<  Reference density.
+    double _temperature0;  ///<  Reference temperature.
+    double _beta;          ///<  Parameter.
+};
+
+}  // end namespace
+}  // end namespace
+
+#endif /* LINEARTEMPERATUREDEPENDENTDENSITY_H */

MaterialLib/Fluid/Density/LiquidDensity.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \author wenqing
+ * \file   LiquidDensity.h
+ *
+ * Created on August 4, 2016, 10:15 AM
+ */
+
+#ifndef LIQUIDDENSITY_H
+#define LIQUIDDENSITY_H
+
+#include <vector>
+
+#include "BaseLib/Error.h"
+
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/**
+ *  Class of a generic density model for liquid fluids varying with temperature
+ *  or pressure.
+ *
+ *  The formula is given on
+ *  <a href="engineeringtoolbox">
+ * http://www.engineeringtoolbox.com/fluid-density-temperature-pressure-d_309.html</a>
+ *   which reads
+ *   \f[
+ *        \rho_l = \rho_0/(1+\beta(T-T_0))/(1-(p-p_0)/E)
+ *   \f]
+ *   where
+ *    \f{eqnarray*}{
+ *       &\rho_l:& \mbox{liquid density,}\\
+ *       &\rho_0:& \mbox{initial liquid density,}\\
+ *       &\beta: & \mbox{volumetric temperature expansion coefficient,}\\
+ *       &T:&      \mbox{temperature,}\\
+ *       &T_0:&    \mbox{initial temperature,}\\
+ *       &p:&      \mbox{pressure,}\\
+ *       &p_0:&    \mbox{initial pressure,}\\
+ *       &E:&      \mbox{bulk modulus fluid elasticity.}\\
+ *    \f}
+ */
+class LiquidDensity : public FluidProperty
+{
+public:
+    /// \param config  ConfigTree object which contains the input data
+    ///                 including  <type>fluid</type> and it has
+    ///                 a tag of <density>
+    LiquidDensity(BaseLib::ConfigTree const& config)
+        :  //! \ogs_file_param{material__fluid__density__liquid_density__beta}
+    _beta(config.getConfigParameter<double>("beta")),
+          //! \ogs_file_param{material__fluid__density__liquid_density__rho0}
+    _rho0(config.getConfigParameter<double>("rho0")),
+          //! \ogs_file_param{material__fluid__density__liquid_density__temperature0}
+    _temperature0(config.getConfigParameter<double>("temperature0")),
+          //! \ogs_file_param{material__fluid__density__liquid_density__p0}
+    _p0(config.getConfigParameter<double>("p0")),
+          //! \ogs_file_param{material__fluid__density__liquid_density__bulk_modulus}
+    _bulk_moudlus(config.getConfigParameter<double>("bulk_modulus")),
+          _derivative_functions{&LiquidDensity::dLiquidDensity_dT,
+                                &LiquidDensity::dLiquidDensity_dp}
+    {
+    }
+
+    /// Get density model name.
+    virtual std::string getName() const final
+    {
+        return "Temperature or pressure dependent liquid density";
+    }
+
+    /// Get density value.
+    /// \param var_vals Variable values in an array. The order of its elements
+    ///                 is given in enum class PropertyVariable.
+    virtual double getValue(const double var_vals[]) const final
+    {
+        const double T = var_vals[static_cast<int>(PropertyVariable::T)];
+        const double p = var_vals[static_cast<int>(PropertyVariable::pl)];
+        return _rho0 / (1. + _beta * (T - _temperature0)) /
+               (1. - (p - _p0) / _bulk_moudlus);
+    }
+
+    /// Get the partial differential of the density with respect to temperature
+    /// or liquid pressure.
+    /// \param var_vals  Variable values  in an array. The order of its elements
+    ///                   is given in enum class PropertyVariable.
+    virtual double getdValue(const double var_vals[],
+                             const PropertyVariable var) const final
+    {
+        assert(var == PropertyVariable::T || var == PropertyVariable::pl);
+
+        const int func_id = (var == PropertyVariable::T) ? 0 : 1;
+        const double T = var_vals[static_cast<int>(PropertyVariable::T)];
+        const double p = var_vals[static_cast<int>(PropertyVariable::pl)];
+        return (this->*_derivative_functions[func_id])(T, p);
+    }
+
+private:
+    /// Volumetric temperature expansion coefficient.
+    double _beta;
+
+    /// Initial density.
+    double _rho0;
+
+    /// Initial temperature.
+    double _temperature0;
+
+    /// Initial pressure.
+    double _p0;
+
+    /// Bulk modulus.
+    double _bulk_moudlus;
+
+    /**
+     *  Calculate the derivative of density of fluids with respect to
+     * temperature.
+     *   \param T    Temperature (K).
+     *   \param p    Pressure (Pa).
+     */
+    double dLiquidDensity_dT(const double T, const double p) const
+    {
+        const double fac_T = 1. + _beta * (T - _temperature0);
+        return -_beta * _rho0 / (fac_T * fac_T) /
+               (1. - (p - _p0) / _bulk_moudlus);
+    }
+
+    /**
+     *  Calculate the derivative of density of fluids with respect to pressure.
+     *   \param T    Temperature (K).
+     *   \param p    Pressure (Pa).
+     */
+    double dLiquidDensity_dp(const double T, const double p) const
+    {
+        const double fac_p = 1. - (p - _p0) / _bulk_moudlus;
+        return _rho0 / (1. + _beta * (T - _temperature0)) /
+               (fac_p * fac_p * _bulk_moudlus);
+    }
+
+    typedef double (LiquidDensity::*ptr2_derivative_f)(const double,
+                                                       const double) const;
+
+    /// An array of pointer to derivative functions.
+    ptr2_derivative_f _derivative_functions[2];
+};
+
+}  // end of namespace
+}  // end of namespace
+
+#endif /* LIQUIDDENSITY_H */

MaterialLib/Fluid/Density/createFluidDensityModel.cpp

+/*!
+   \file  createFluidDensityModel.cpp
+   \brief create an instance of a fluid density class.
+
+   \copyright
+    Copyright (c) 2012-2016, 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 "createFluidDensityModel.h"
+
+#include "BaseLib/Error.h"
+
+#include "MaterialLib/Fluid/ConstantFluidProperty.h"
+#include "IdealGasLaw.h"
+#include "LinearTemperatureDependentDensity.h"
+#include "LiquidDensity.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+std::unique_ptr<FluidProperty> createFluidDensityModel(BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{material__fluid__density__type}
+    auto const type = config.getConfigParameter<std::string>("type");
+
+    if (type == "constant")
+    {
+        return std::unique_ptr<FluidProperty>(new ConstantFluidProperty(
+            //! \ogs_file_param{material__fluid__density__constant_value}
+            config.getConfigParameter<double>("value")) );
+    }
+    //! \ogs_file_param{material__fluid__density__LiquidDensity}
+    else if (type == "LiquidDensity")
+        return std::unique_ptr<FluidProperty>(new LiquidDensity(config));
+    //! \ogs_file_param{material__fluid__density__TemperatureDependent}
+    else if (type == "TemperatureDependent")
+        return std::unique_ptr<FluidProperty>(
+                new LinearTemperatureDependentDensity(config));
+    //! \ogs_file_param{material__fluid__density__IdealGasLaw}
+    else if (type == "IdealGasLaw")
+    {
+        //! \ogs_file_param{material__fluid__density__IdealGasLaw__molar_mass}
+        return std::unique_ptr<FluidProperty>(new IdealGasLaw(
+            config.getConfigParameter<double>("molar_mass")) );
+    }
+    else
+    {
+        OGS_FATAL(
+            "The density type %s is unavailable.\n", type.data(),
+            "The available types are liquid_density "
+            "temperature_dependent and ideal_gas_law");
+    }
+}
+
+}  // end namespace
+}  // end namespace

MaterialLib/Fluid/Density/createFluidDensityModel.h

+/*!
+   \file  createFluidDensityModel.h
+   \brief create an instance of a fluid density class.
+
+   \copyright
+    Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+               Distributed under a Modified BSD License.
+               See accompanying file LICENSE.txt or
+               http://www.opengeosys.org/project/license
+*/
+#ifndef CREATE_FLUID_DENSITY_MODEL_H_
+#define CREATE_FLUID_DENSITY_MODEL_H_
+
+#include <memory>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/// Create a density model
+/// \param config  ConfigTree object has a tag of <density>
+std::unique_ptr<FluidProperty> createFluidDensityModel(BaseLib::ConfigTree const& config);
+}
+}  // end namespace
+#endif

MaterialLib/Fluid/FluidProperty.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   FluidProperty.h
+ *
+ * Created on August 12, 2016, 3:34 PM
+ */
+
+#ifndef FLUIDPROPERTY_H
+#define FLUIDPROPERTY_H
+
+#include <string>
+
+namespace MaterialLib
+{
+namespace Fluid
+{
+/// Variable that determine the property.
+enum class PropertyVariable
+{
+    T = 0,   ///< temperature
+    pl = 1,  ///< pressure of the liquid phase (1st phase for some cases))
+    pg = 2,  ///< pressure of the gas phase (2nd phase for some cases))
+};
+
+/// Base class of fluid density properties
+class FluidProperty
+{
+public:
+
+    /// Get model name.
+    virtual std::string getName() const = 0;
+
+    /// Get property value.
+    /// \param var_vals Variable values in an array. The order of its elements
+    ///                 is given in enum class PropertyVariable.
+    virtual double getValue(const double /* var_vals*/[]) const = 0;
+
+    /// Get the partial differential of the property value
+    /// \param var_vals  Variable values  in an array. The order of its elements
+    ///                   is given in enum class PropertyVariable.
+    virtual double getdValue(const double /* var_vals*/[],
+                             const PropertyVariable /* var */) const = 0;
+};
+
+}  // end namespace
+}  // end namespace
+
+#endif /* FLUIDPROPERTY_H */

MaterialLib/Fluid/FluidPropertyHeaders.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * File:   FluidPropertyHeaders.h
+ * Author: wenqing
+ *
+ * Created on August 18, 2016, 11:10 AM
+ */
+
+#ifndef FLUIDPROPERTYHEADERS_H
+#define FLUIDPROPERTYHEADERS_H
+
+#include "Density/createFluidDensityModel.h"
+#include "Viscosity/createViscosityModel.h"
+
+#endif /* FLUIDPROPERTYHEADERS_H */
+

Tests/MaterialLib/TestFluidDensity.cpp

+/*!
+   \file  TestFluidDensity.cpp
+   \brief Test classes for fluid density models.
+
+   \author Wenqing Wang
+   \date Jan 2015
+
+   \copyright
+    Copyright (c) 2012-2016, 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 <memory>
+
+#include "TestTools.h"
+
+#include "MaterialLib/Fluid/Density/createFluidDensityModel.h"
+#include "MaterialLib/PhysicalConstant.h"
+
+namespace
+{
+using namespace MaterialLib;
+using namespace MaterialLib::Fluid;
+
+//----------------------------------------------------------------------------
+// Test density models.
+std::unique_ptr<FluidProperty> createTestFluidDensityModel(const char xml[])
+{
+    auto const ptree = readXml(xml);
+    BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
+                             BaseLib::ConfigTree::onwarning);
+    auto const& sub_config = conf.getConfigSubtree("density");
+    return MaterialLib::Fluid::createFluidDensityModel(sub_config);
+}
+
+TEST(Material, checkConstantDensity)
+{
+    const char xml[] =
+        "<density>"
+        "   <type>constant</type>"
+        "   <value> 998.0 </value> "
+        "</density>";
+    const auto rho = createTestFluidDensityModel(xml);
+
+    ASSERT_EQ(998.0, rho->getValue(nullptr));
+    ASSERT_EQ(0.0,
+              rho->getdValue(nullptr, MaterialLib::Fluid::PropertyVariable::T));
+}
+
+TEST(Material, checkIdealGasLaw)
+{
+    const char xml[] =
+        "<density>"
+        "   <type>IdealGasLaw</type>"
+        "   <molar_mass> 28.96 </molar_mass> "
+        "</density>";
+    const auto rho = createTestFluidDensityModel(xml);
+
+    const double molar_air = 28.96;
+    const double T = 290.;
+    const double p = 1.e+5;
+    const double R = PhysicalConstant::IdealGasConstant;
+    const double expected_air_dens = molar_air * p / (R * T);
+    double vars[] = {290, 0, 1.e+5};
+    ASSERT_NEAR(expected_air_dens, rho->getValue(vars), 1.e-10);
+
+    const double expected_d_air_dens_dT = -molar_air * p / (R * T * T);
+    ASSERT_NEAR(expected_d_air_dens_dT,
+                rho->getdValue(vars, Fluid::PropertyVariable::T), 1.e-10);
+
+    const double expected_d_air_dens_dp = molar_air / (R * T);
+    ASSERT_NEAR(expected_d_air_dens_dp,
+                rho->getdValue(vars, Fluid::PropertyVariable::pg), 1.e-10);
+}
+
+TEST(Material, checkLinearTemperatureDependentDensity)
+{
+    const char xml[] =
+        "<density>"
+        "   <type>TemperatureDependent</type>"
+        "   <temperature0> 293.0 </temperature0> "
+        "   <beta> 4.3e-4 </beta> "
+        "   <rho0>1000.</rho0>"
+        "</density>";
+
+    const auto rho = createTestFluidDensityModel(xml);
+
+    double vars[] = {273.1 + 60.0};
+    ASSERT_NEAR(1000.0 * (1 + 4.3e-4 * (vars[0] - 293.0)), rho->getValue(vars),
+                1.e-10);
+    ASSERT_NEAR(1000.0 * 4.3e-4,
+                rho->getdValue(vars, Fluid::PropertyVariable::T), 1.e-10);
+}
+
+TEST(Material, checkLiquidDensity)
+{
+    const char xml[] =
+        "<density>"
+        "   <type>LiquidDensity</type>"
+        "   <temperature0> 273.15 </temperature0> "
+        "   <p0> 1.e+5 </p0> "
+        "   <bulk_modulus> 2.15e+9 </bulk_modulus> "
+        "   <beta> 2.0e-4 </beta> "
+        "   <rho0>999.8</rho0>"
+        "</density>";
+    const auto rho = createTestFluidDensityModel(xml);
+
+    const double vars[] = {273.15 + 60.0, 1.e+6};
+    const double T0 = 273.15;
+    const double p0 = 1.e+5;
+    const double rho0 = 999.8;
+    const double K = 2.15e+9;
+    const double beta = 2.e-4;
+    const double T = vars[0];
+    const double p = vars[1];
+
+    const double fac_T = 1. + beta * (T - T0);
+    ASSERT_NEAR(rho0 / fac_T / (1. - (p - p0) / K), rho->getValue(vars),
+                1.e-10);
+
+    // Test the derivative with respect to temperature.
+    ASSERT_NEAR(-beta * rho0 / (fac_T * fac_T) / (1. - (p - p0) / K),
+                rho->getdValue(vars, Fluid::PropertyVariable::T), 1.e-10);
+
+    // Test the derivative with respect to pressure.
+    const double fac_p = 1. - (p - p0) / K;
+    ASSERT_NEAR(rho0 / (1. + beta * (T - T0)) / (fac_p * fac_p * K),
+                rho->getdValue(vars, Fluid::PropertyVariable::pl), 1.e-10);
+}
+
+}