The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

Commit fcd50020 authored by Lukas Riedel's avatar Lukas Riedel

Support setting porosity scaling factors in KnoFu interface

* Use enum to identify target scaling factor.
* Improve error messages.
parent 10c6c9db
......@@ -113,6 +113,8 @@ public:
= std::function<ScalingValue(const std::vector<ScalingValue>&,
const std::vector<Kernel>&,
const Domain&)>;
/// Scaling factors that can be manipulated with manipulate_scaling()
enum class TargetScalingFactor {head, cond, por};
private:
......@@ -215,27 +217,13 @@ public:
* \param values Scaling values to set
* \param kernels Interpolation kernels to use
*/
void manipulate_scaling (const std::string_view target_factor_str,
void manipulate_scaling (const TargetScalingFactor target_factor,
const KernelWeightApplicator weight_f,
const std::vector<ScalingValue>& values,
const std::vector<Kernel>& kernels)
{
assert(values.size() == kernels.size());
// translate string to enum for faster application in loop
enum class TargetFactor {head, cond, por};
TargetFactor target_factor;
if (target_factor_str == "head") {
target_factor = TargetFactor::head;
}
else if (target_factor_str == "cond") {
target_factor = TargetFactor::cond;
}
else {
DUNE_THROW(Exception, "Unknown target scaling factor");
}
// iterate over all cells
for (const auto& cell : elements(_gv))
{
......@@ -247,12 +235,16 @@ public:
const auto index = _mapper.index(cell);
auto& scaling = std::get<Scaling>(_param[index]);
if (target_factor == TargetFactor::head) {
if (target_factor == TargetScalingFactor::head) {
scaling.scale_head = s_value;
}
else if (target_factor == TargetFactor::cond) {
else if (target_factor == TargetScalingFactor::cond) {
scaling.scale_cond = s_value;
}
// Last case: TargetScalingFactor::por
else {
scaling.scale_por = s_value;
}
}
}
......
......@@ -116,6 +116,8 @@ public:
using BackendVector = typename Dune::PDELab::Backend::Native<U>;
/// Type of scaling value to insert
using ScalingValue = typename Simulation::Traits::FlowParameters::ScalingValue;
/// Names of scaling factors to manipulate
using TargetScalingFactor = typename Parameters::TargetScalingFactor;
/// Construct this object
KnoFuInterface (
......@@ -261,7 +263,7 @@ public:
* \param length_scales Length scales of scaling values to use in the
* kernel interpolation.
*/
void set_scaling_factor (const std::string_view target_factor,
void set_scaling_factor (const TargetScalingFactor target_factor,
const std::string_view kernel_name,
const std::string_view weighting,
const std::vector<ScalingValue>& values,
......@@ -286,9 +288,15 @@ public:
};
}
else {
this->_log->error("Unknown manipulation kernel name: {}",
kernel_name);
DUNE_THROW(IOError, "Unknown manipulation kernel name");
}
// Set the default scaling value (no scaling)
const RF s_default =
(target_factor == TargetScalingFactor::por) ? 0.0 : 1.0;
// bind kernels to the positions
using Kernel = typename Simulation::Traits::FlowParameters::Kernel;
std::vector<Kernel> kernels;
......@@ -306,13 +314,13 @@ public:
// None: Apply kernels independently
if (weighting == "none") {
weighting_f = [](const std::vector<ScalingValue>& values,
const std::vector<Kernel>& kernels,
const Domain& pos)
weighting_f = [s_default](const std::vector<ScalingValue>& values,
const std::vector<Kernel>& kernels,
const Domain& pos)
{
ScalingValue s_value = 1.0;
ScalingValue s_value = s_default;
for (size_t i = 0; i < values.size(); ++i) {
s_value += kernels[i](pos) * (values[i] - 1.0);
s_value += kernels[i](pos) * (values[i] - s_default);
}
return s_value;
};
......@@ -339,21 +347,21 @@ public:
using Vector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using SuperLU = Dune::SuperLU<Matrix>;
weighting_f = [](const std::vector<ScalingValue>& values,
const std::vector<Kernel>& kernels,
const Domain& pos)
weighting_f = [s_default](const std::vector<ScalingValue>& values,
const std::vector<Kernel>& kernels,
const Domain& pos)
{
// initialize objects
const auto size = values.size();
Matrix corr_mat(size, size, size, 1.0,
Matrix::BuildMode::implicit);
Vector weights(size), miller(size);
Vector weights(size), scale(size);
// fill Miller scaling vectors and correlation matrix
// fill scaling vectors and correlation matrix
for (size_t i = 0; i < size; ++i) {
corr_mat.entry(i, i) = 1.0;
weights[i] = kernels[i](pos);
miller[i] = values[i] - 1.0;
scale[i] = values[i] - s_default;
// compute kernel correlations
for (size_t j = 0; j < i; ++j) {
const auto corr = kernels[i](pos)
......@@ -368,13 +376,15 @@ public:
SuperLU solver(corr_mat);
Dune::InverseOperatorResult info;
Vector result(size);
solver.apply(result, miller, 1e-6, info);
solver.apply(result, scale, 1e-6, info);
// compute the dot product
return 1.0 + Dune::dot(weights, result);
return s_default + Dune::dot(weights, result);
};
}
else {
this->_log->error("Unknown scaling factor kernel weighting: {}",
weighting);
DUNE_THROW(IOError, "Unknown scaling factor kernel weighting");
}
......
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