Commit 38c78add authored by Lukas Riedel's avatar Lukas Riedel 🎧
Browse files

Add weighting options for scaling factor manipulation

Add options 'none', 'norm', and 'correlation', which set functions
that calculate the new scaling field based on the chosen kernel(s).

* 'none': Simple superposition.
* 'norm': Weighted superposition.
* 'correlation': Correlated sum of scaling factors. Use kernels as
  local correlation. This involves solving a linear equation system.

The call to `set_scaling_factor` in the KnoFuInterface must now
specify a weighting type.
parent 87b24ba9
......@@ -60,6 +60,7 @@ private:
LevelGridView
>;
using Vector = typename Traits::Vector;
using Domain = typename Traits::Domain;
// Define the storage and cache for parameterizations and scaling
/// Index type used by the element mapper
......@@ -106,7 +107,12 @@ public:
/// Type of a single scaling value
using ScalingValue = RF;
/// Type of a scaling manipulation
using Kernel = std::function<ScalingValue(const typename Traits::Domain&)>;
using Kernel = std::function<ScalingValue(const Domain&)>;
/// Type of scaling manipulation applicator
using KernelWeightApplicator
= std::function<ScalingValue(const std::vector<ScalingValue>&,
const std::vector<Kernel>&,
const Domain&)>;
private:
......@@ -205,38 +211,47 @@ public:
* to evaluate the target factor throughout the entire physical space.
*
* \param target_factor Name of the target scaling factor
* \param weight_f Function for applying kernels and respective weights
* \param values Scaling values to set
* \param kernels Interpolation kernels to use
*/
void manipulate_scaling (const std::string_view target_factor,
void manipulate_scaling (const std::string_view target_factor_str,
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))
{
const auto center = cell.geometry().center();
// evaluate scaling kernels and weighting
const auto pos = cell.geometry().center();
const auto s_value = weight_f(values, kernels, pos);
// apply to stored scaling value
const auto index = _mapper.index(cell);
auto& scaling = std::get<Scaling>(_param[index]);
// calculate scaling value
ScalingValue s_value = 0.0;
for (size_t i = 0; i < values.size(); ++i) {
const auto kern = kernels[i](center);
s_value += kern * (values[i] - 1.0);
}
// apply to scaling value
if (target_factor == "head") {
scaling.scale_head = 1.0 + s_value;
}
else if (target_factor == "cond") {
scaling.scale_cond = 1.0 + s_value;
if (target_factor == TargetFactor::head) {
scaling.scale_head = s_value;
}
else {
DUNE_THROW(Exception, "Unknown target scaling factor");
else if (target_factor == TargetFactor::cond) {
scaling.scale_cond = s_value;
}
}
}
......
......@@ -14,6 +14,13 @@
#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solver.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/common/function.hh>
#include <dune/pdelab/constraints/common/constraintstransformation.hh>
......@@ -227,6 +234,8 @@ public:
* \param target_factor Name of the scaling factor to manipulate
* \param kernal_name Name of the interpolation kernel to use.
* Options: `gaussian`.
* \param weighting Type of weighting applied to kernels to receive the
* final scaling field.
* \param values List of scaling values to set.
* \param positions Positions of scaling values on the grid.
* \param length_scales Length scales of scaling values to use in the
......@@ -234,6 +243,7 @@ public:
*/
void set_scaling_factor (const std::string_view target_factor,
const std::string_view kernel_name,
const std::string_view weighting,
const std::vector<ScalingValue>& values,
const std::vector<Domain>& positions,
const std::vector<DF>& length_scales)
......@@ -241,31 +251,116 @@ public:
assert(values.size() == positions.size()
&& values.size() == length_scales.size());
using Kernel = typename Simulation::Traits::FlowParameters::Kernel;
if (kernel_name != "gaussian") {
DUNE_THROW(Exception, "Unknown manipulation kernel name");
// Choose the kernel function
using FullKernel = std::function<ScalingValue(const Domain&,
const Domain&,
const DF)>;
FullKernel kernel;
if (kernel_name == "gaussian") {
kernel = [](const Domain& x, const Domain& mu, const DF sigma){
const auto deviation = x - mu;
const RF exponent = - deviation.two_norm2()
/ ( 2.0 * std::pow(sigma, 2) );
return std::exp(exponent);
};
}
else {
DUNE_THROW(IOError, "Unknown manipulation kernel name");
}
auto kernel = [](const auto x, const auto mu, const auto sigma){
const auto deviation = x - mu;
const RF exponent = - deviation.two_norm2()
/ ( 2.0 * std::pow(sigma, 2) );
return std::exp(exponent);
};
// bind kernels to the positions
using Kernel = typename Simulation::Traits::FlowParameters::Kernel;
std::vector<Kernel> kernels;
for (size_t i = 0; i < values.size(); ++i) {
using namespace std::placeholders;
kernels.push_back(std::bind(kernel,
_1,
std::placeholders::_1,
std::cref(positions[i]),
std::cref(length_scales[i])));
}
// Choose the weighting function
using KernelWeightApplicator = typename Simulation::Traits
::FlowParameters::KernelWeightApplicator;
KernelWeightApplicator weighting_f;
// None: Apply kernels independently
if (weighting == "none") {
weighting_f = [](const std::vector<ScalingValue>& values,
const std::vector<Kernel>& kernels,
const Domain& pos)
{
ScalingValue s_value = 1.0;
for (size_t i = 0; i < values.size(); ++i) {
s_value += kernels[i](pos) * (values[i] - 1.0);
}
return s_value;
};
}
// Norm: Compute a norm from all kernel contributions
else if (weighting == "norm") {
weighting_f = [](const std::vector<ScalingValue>& values,
const std::vector<Kernel>& kernels,
const Domain& pos)
{
ScalingValue s_value = 0.0, s_norm = 0.0;
for (size_t i = 0; i < values.size(); ++i) {
const auto kern = kernels[i](pos);
s_value += kern * values[i];
s_norm += kern;
}
s_value /= s_norm;
return s_value;
};
}
// Correlation: Use kernel contributions as correlation information
else if (weighting == "correlation") {
using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
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)
{
// initialize objects
const auto size = values.size();
Matrix corr_mat(size, size, size, 1.0,
Matrix::BuildMode::implicit);
Vector weights(size), miller(size);
// fill Miller 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;
// compute kernel correlations
for (size_t j = 0; j < i; ++j) {
const auto corr = kernels[i](pos)
* kernels[j](pos);
corr_mat.entry(i, j) = corr;
corr_mat.entry(j, i) = corr;
}
}
// solve the equation system
corr_mat.compress();
SuperLU solver(corr_mat);
Dune::InverseOperatorResult info;
Vector result(size);
solver.apply(result, miller, 1e-6, info);
// compute the dot product
return 1.0 + Dune::dot(weights, result);
};
}
else {
DUNE_THROW(IOError, "Unknown scaling factor kernel weighting");
}
// pass to the parameter interface
this->fparam->manipulate_scaling(target_factor,
weighting_f,
values,
kernels);
}
......
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