Revision 62cbe69c65fa6b7ad56ddc7fde5fe7e98a5895be authored by Norbert Grunwald on 17 March 2021, 13:07:44 UTC, committed by Norbert Grunwald on 01 April 2021, 07:28:49 UTC
1 parent f7ca4da
TestMPLRelPermVanGenuchten.cpp
/**
* \file
* \copyright
* Copyright (c) 2012-2021, 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/Medium.h"
#include "MaterialLib/MPL/Properties/RelativePermeability/RelPermVanGenuchten.h"
#include "TestMPL.h"
#include "Tests/TestTools.h"
namespace MPL = MaterialPropertyLib;
TEST(MaterialPropertyLib, RelPermVanGenuchten)
{
double const residual_liquid_saturation = 0.1;
double const residual_gas_saturation = 0.05;
double const exponent = 0.79;
double const min_relative_permeability_liquid = 1e-12;
MPL::Property const& permeability = MPL::RelPermVanGenuchten{
"relative_permeability", residual_liquid_saturation,
residual_gas_saturation, min_relative_permeability_liquid, exponent};
MPL::VariableArray variable_array;
ParameterLib::SpatialPosition const pos;
double const t = std::numeric_limits<double>::quiet_NaN();
double const dt = std::numeric_limits<double>::quiet_NaN();
double const S_0 = -0.1;
double const S_max = 1.1;
int const n_steps = 1000;
for (int i = 0; i <= n_steps; ++i)
{
double const S_L = S_0 + i * (S_max - S_0) / n_steps;
variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L;
double const k_rel =
permeability.template value<double>(variable_array, pos, t, dt);
double const dk_rel = permeability.template dValue<double>(
variable_array, MPL::Variable::liquid_saturation, pos, t, dt);
double const eps = 1e-8;
variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L - eps;
double const k_rel_minus =
permeability.template value<double>(variable_array, pos, t, dt);
variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
S_L + eps;
double const k_rel_plus =
permeability.template value<double>(variable_array, pos, t, dt);
double const Dk_rel = (k_rel_plus - k_rel_minus) / 2 / eps;
if (S_L < 1 - residual_gas_saturation)
{
ASSERT_LE(std::abs(dk_rel - Dk_rel), 1e-7)
<< "for saturation " << S_L << " and relative permeability "
<< k_rel;
}
else
{
ASSERT_EQ(dk_rel, 0.) << "for saturation " << S_L
<< " and relative permeability " << k_rel;
}
}
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...