TestMPLExponential.cpp
/**
* \file
*
* \copyright
* Copyright (c) 2012-2023, 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 <autocheck/autocheck.hpp>
#include <cmath>
#include "MaterialLib/MPL/Properties/Exponential.h"
#include "Tests/AutoCheckTools.h"
namespace ac = autocheck;
namespace MPL = MaterialPropertyLib;
struct MaterialPropertyLibExponentialProperty : public ::testing::Test
{
void SetUp() override
{
double const y_offset = 1.;
double const y_ref = 1.;
double const T_ref = 1.;
double const factor = 2.;
p = std::make_unique<MPL::Exponential>(
"exponential", y_offset, y_ref,
MPL::ExponentData{MPL::Variable::temperature, T_ref, factor});
}
std::unique_ptr<MPL::Property> p;
ac::gtest_reporter gtest_reporter;
};
// First order derivative approximated with second order central differences.
template <typename F>
double dydx_C2(double const x, F&& f, double const h)
{
return (f(x + h) - f(x - h)) / (2 * h);
}
// Second order derivative approximated with second order central differences.
template <typename F>
double d2ydx2_C2(double const x, F&& f, double const h)
{
return (f(x + h) - 2 * f(x) + f(x - h)) / (h * h);
}
TEST_F(MaterialPropertyLibExponentialProperty, TestNumericalDerivatives)
{
double const eps = 1e-5;
ParameterLib::SpatialPosition const pos;
double const time = std::numeric_limits<double>::quiet_NaN();
double const dt = std::numeric_limits<double>::quiet_NaN();
auto y = [&](double const T)
{
MPL::VariableArray variable_array;
variable_array.temperature = T;
return p->template value<double>(variable_array, pos, time, dt);
};
auto f = [&](double const T)
{
MPL::VariableArray variable_array;
variable_array.temperature = T;
double const v =
p->template value<double>(variable_array, pos, time, dt);
double const dv = p->template dValue<double>(
variable_array, MPL::Variable::temperature, pos, time, dt);
double const dv2 = p->template d2Value<double>(
variable_array, MPL::Variable::temperature,
MPL::Variable::temperature, pos, time, dt);
double const Dv = dydx_C2(T, y, eps);
double const Dv2 = d2ydx2_C2(T, y, eps);
if ((std::abs(dv - Dv) > 1e-9 * v) ||
(std::abs(dv2 - Dv2) > 1.5e-4 * v))
INFO("{} {} {}", T, std::abs(dv - Dv) / v, std::abs(dv2 - Dv2) / v);
return (std::abs(dv - Dv) <= 1e-9 * v) &&
(std::abs(dv2 - Dv2) <= 1.5e-4 * v);
};
// Limit values to avoid +-inf.
auto gen = ac::IntervalGenerator(-100., 100.);
ac::check<double>(f, 10000, ac::make_arbitrary(gen), gtest_reporter);
}
TEST(MaterialPropertyLib, Exponential)
{
double const y_offset = -1e-3;
double const y_ref = 1e-3;
double const reference_condition = 20.0;
double const factor = 1 / 75.0;
MPL::ExponentData const exp_data{MPL::Variable::temperature,
reference_condition, factor};
MPL::Property const& p =
MPL::Exponential{"exponential", y_offset, y_ref, exp_data};
double const T = 20.;
MPL::VariableArray variable_array;
variable_array.temperature = T;
ParameterLib::SpatialPosition const pos;
double const time = std::numeric_limits<double>::quiet_NaN();
double const dt = std::numeric_limits<double>::quiet_NaN();
ASSERT_NEAR(
p.template value<double>(variable_array, pos, time, dt),
y_offset + y_ref * (std::exp(factor * (T - reference_condition))),
1.e-10);
ASSERT_EQ(p.template dValue<double>(
variable_array, MPL::Variable::phase_pressure, pos, time, dt),
0.0);
ASSERT_NEAR(p.template dValue<double>(
variable_array, MPL::Variable::temperature, pos, time, dt),
y_ref * factor * std::exp(factor * (T - reference_condition)),
1.e-16);
ASSERT_NEAR(
p.template d2Value<double>(variable_array, MPL::Variable::temperature,
MPL::Variable::temperature, pos, time, dt),
y_ref * std::pow(factor, 2) *
std::exp(factor * (T - reference_condition)),
1.e-16);
}