Commit 64c9e1cc authored by Norbert Grunwald's avatar Norbert Grunwald
Browse files

[T/MPL] New SaturationBrooksCorey test based on numerical derivatives.

parent a6456e0e
Pipeline #6763 failed with stages
in 76 minutes and 48 seconds
/**
* \file
*
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
......@@ -10,74 +9,62 @@
*/
#include <gtest/gtest.h>
#include <sstream>
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationBrooksCorey.h"
#include "TestMPL.h"
#include "Tests/TestTools.h"
namespace MPL = MaterialPropertyLib;
TEST(MaterialPropertyLib, SaturationBrooksCorey)
{
const double ref_lambda = 2.5;
const double ref_residual_liquid_saturation = 0.12;
const double ref_residual_gas_saturation = 0.06;
const double ref_entry_pressure = 5678.54;
const double max_saturation = 1. - ref_residual_gas_saturation;
double const residual_liquid_saturation = 0.1;
double const residual_gas_saturation = 0.05;
double const exponent = 3.0;
double const p_b = 5000;
MPL::Property const& pressure_saturation =
MPL::SaturationBrooksCorey{"saturation", residual_liquid_saturation,
residual_gas_saturation, exponent, p_b};
std::stringstream m;
m << "<medium>\n";
m << "<phases></phases>\n";
m << "<properties>\n";
m << " <property>\n";
m << " <name>saturation</name>\n";
m << " <type>SaturationBrooksCorey</type>\n";
m << " <residual_liquid_saturation>" << ref_residual_liquid_saturation
<< "</residual_liquid_saturation>\n";
m << " <residual_gas_saturation>" << ref_residual_gas_saturation
<< "</residual_gas_saturation>\n";
m << " <lambda>" << ref_lambda << "</lambda>\n";
m << " <entry_pressure>" << ref_entry_pressure << "</entry_pressure>\n";
m << " </property> \n";
m << "</properties>\n";
m << "</medium>\n";
auto const& medium = Tests::createTestMaterial(m.str());
MaterialPropertyLib::VariableArray variable_array;
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();
for (double p_cap = 1.0; p_cap < 1.0e10; p_cap *= 1.5)
double const p_0 = -1e6;
double const p_max = 4000;
int const n_steps = 10000;
for (int i = 0; i <= n_steps; ++i)
{
variable_array[static_cast<int>(
MaterialPropertyLib::Variable::capillary_pressure)] = p_cap;
double const p_L = p_0 + i * (p_max - p_0) / n_steps;
variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
-p_L;
double const S = pressure_saturation.template value<double>(
variable_array, pos, t, dt);
double const dS = pressure_saturation.template dValue<double>(
variable_array, MPL::Variable::capillary_pressure, pos, t, dt);
double const dS2 = pressure_saturation.template d2Value<double>(
variable_array, MPL::Variable::capillary_pressure,
MPL::Variable::capillary_pressure, pos, t, dt);
auto s_L =
medium->property(MaterialPropertyLib::PropertyType::saturation)
.template value<double>(variable_array, pos, t, dt);
auto ds_L_dp_cap =
medium->property(MaterialPropertyLib::PropertyType::saturation)
.template dValue<double>(
variable_array,
MaterialPropertyLib::Variable::capillary_pressure, pos, t,
dt);
double const eps = 1e-2;
variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
-p_L - eps;
double const S_minus = pressure_saturation.template value<double>(
variable_array, pos, t, dt);
variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
-p_L + eps;
double const S_plus = pressure_saturation.template value<double>(
variable_array, pos, t, dt);
const double s_eff =
std::pow(ref_entry_pressure / std::max(p_cap, ref_entry_pressure),
ref_lambda);
const double s_ref =
s_eff * (max_saturation - ref_residual_liquid_saturation) +
ref_residual_liquid_saturation;
const double ds_eff_dpc =
-ref_lambda / std::max(p_cap, ref_entry_pressure) * s_ref;
const double ds_L_ds_eff =
1. / (max_saturation - ref_residual_liquid_saturation);
const double ds_L_dpc = ds_L_ds_eff * ds_eff_dpc;
double const DS = (S_plus - S_minus) / 2 / eps;
double const DS2 = (S_plus - 2 * S + S_minus) / (eps * eps);
ASSERT_NEAR(s_L, s_ref, 1.e-10);
ASSERT_NEAR(ds_L_dp_cap, ds_L_dpc, 1.e-10);
ASSERT_LE(std::abs(dS - DS), 1e-9)
<< "for capillary pressure " << -p_L << " and saturation " << S;
ASSERT_LE(std::abs(dS2 - DS2), 1e-9)
<< "for capillary pressure " << -p_L << " and saturation " << S;
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment