Commit 78ccc14e authored by Lukas Riedel's avatar Lukas Riedel

Optimize new MvG parameterization wrt speed

* Remove static lambdas
* Remove scale operation if no scaling is applied
parent 9d0a7179
......@@ -105,16 +105,19 @@ public:
{
verify_cache();
auto& par = std::get<std::shared_ptr<Parameterization>>(_cache.second);
auto& type = std::get<ScalingType>(_cache.second);
auto cond_f = par->conductivity_f();
using Saturation = typename Parameterization::Saturation;
RF scale = 1.0;
if (std::get<ScalingType>(_cache.second) == ScalingType::Miller) {
scale = std::get<ScalingFactor>(_cache.second).value;
if (type == ScalingType::Miller) {
auto& scale = std::get<ScalingFactor>(_cache.second).value;
return [cond_f, &scale](const RF saturation){
return cond_f(Saturation{saturation}).value * scale * scale;
};
}
using Saturation = typename Parameterization::Saturation;
return [&par, scale](const RF saturation){
return par->conductivity_f()(Saturation{saturation}).value
* scale * scale;
return [cond_f](const RF saturation){
return cond_f(Saturation{saturation}).value;
};
}
......@@ -123,15 +126,19 @@ public:
{
verify_cache();
auto& par = std::get<std::shared_ptr<Parameterization>>(_cache.second);
auto& type = std::get<ScalingType>(_cache.second);
auto sat_f = par->saturation_f();
using MatricHead = typename Parameterization::MatricHead;
RF scale = 1.0;
if (std::get<ScalingType>(_cache.second) == ScalingType::Miller) {
scale = std::get<ScalingFactor>(_cache.second).value;
if (type == ScalingType::Miller) {
auto& scale = std::get<ScalingFactor>(_cache.second).value;
return [sat_f, &scale](const RF matric_head){
return sat_f(MatricHead{matric_head * scale}).value;
};
}
using MatricHead = typename Parameterization::MatricHead;
return [&par, scale](const RF matric_head){
return par->saturation_f()(MatricHead{matric_head * scale}).value;
return [sat_f](const RF matric_head){
return sat_f(MatricHead{matric_head}).value;
};
}
......
......@@ -54,40 +54,6 @@ public:
Tau _tau;
N _n;
protected:
inline static std::function<Conductivity(const Saturation,
const SaturatedConductivity,
const N,
const Tau)>
_conductivity_f = [](const Saturation sat,
const SaturatedConductivity k0,
const N n,
const Tau tau) {
const auto m = 1.0 - 1.0 / n.value;
return Conductivity {
k0.value
* std::pow(sat.value, tau.value)
* std::pow(1 - std::pow(1 - std::pow(sat.value, 1.0/m), m), 2)
};
};
inline static std::function<Saturation(const MatricHead,
const Alpha,
const N)>
_saturation_f = [](const MatricHead head,
const Alpha alpha,
const N n) {
if (head.value >= 0.0) {
return Saturation{1.0};
}
const auto m = 1.0 - 1.0 / n.value;
return Saturation {
std::pow(1 + std::pow(alpha.value * head.value, n.value), -m)
};
};
public:
MualemVanGenuchtenParameterization (
const ResidualWaterContent theta_r,
......@@ -109,19 +75,37 @@ public:
std::function<Conductivity(const Saturation)> conductivity_f () const
override
{
return std::bind(_conductivity_f,
std::placeholders::_1,
std::cref(this->_k0), std::cref(_n), std::cref(_tau)
);
const auto m = 1.0 - 1.0 / _n.value;
return [this, m](const Saturation sat) {
return Conductivity {
this->_k0.value
* std::pow(sat.value, _tau.value)
* std::pow(1 - std::pow(1 - std::pow(sat.value, 1.0/m), m), 2)
};
};
}
std::function<Saturation(const MatricHead)> saturation_f () const
override
{
return std::bind(_saturation_f,
std::placeholders::_1,
std::cref(_alpha), std::cref(_n)
);
// return std::bind(_saturation_f,
// std::placeholders::_1,
// std::cref(_alpha), std::cref(_n)
// );
// return [this](const MatricHead head){
// return _saturation_f(head, _alpha, _n);
// };
const auto m = 1.0 - 1.0 / _n.value;
return [this, m](const MatricHead head) {
if (head.value >= 0.0) {
return Saturation{1.0};
}
return Saturation {
std::pow(1 + std::pow(_alpha.value * head.value, _n.value), -m)
};
};
}
};
......
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