Commit f35fd900 authored by Santiago Ospina De Los Ríos's avatar Santiago Ospina De Los Ríos

Merge branch...

Merge branch '142-add-parameterization-for-the-effective-hydrodynamic-dispersion-tensor' into 111-add-dg-local-operator-to-the-simulationtransport-class
parents e7b70443 61a3ecf6
volumes:
sand:
index: 0
water_flow:
richards:
type: MvG
parameters:
alpha: -2.3
n: 4.17
k0: 2.2E-5
theta_r: 0.03
theta_s: 0.31
tau: -1.1
parameters:
alpha: -2.3
n: 4.17
k0: 2.2E-5
theta_r: 0.03
theta_s: 0.31
tau: -1.1
silt:
index: 1
water_flow:
richards:
type: MvG
parameters:
alpha: -0.7
n: 1.3
k0: 1.0E-5
theta_r: 0.01
theta_s: 0.41
tau: 0.0
parameters:
alpha: -0.7
n: 1.3
k0: 1.0E-5
theta_r: 0.01
theta_s: 0.41
tau: 0.0
scaling:
type: None
......@@ -67,6 +67,16 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> endOfTransportStep, endOfRichardsStep, none </values>
</parameter>
<parameter name="writeDispersionTensor">
<definition> Defines whether VTK files should include the hydrodynamic
dispersion tensor. Tensors are written in 3D and have 9 componentents
independently of the world dimension. This can be easily be visualizated
in Paraview with the ``Tensor Glyph`` filter.
</definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="subsamplingLevel">
<definition> Plot VTK files with virtually refined grids. VTK only
supports bilinear triangulations and displays higher-order solutions
......
volumes:
sand:
index: 0
solute:
transport:
type: superposition
eff_diff:
type: milligton_quirk_1
......@@ -11,16 +11,16 @@ volumes:
type: const
trans_disp:
type: const
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
silt:
index: 1
solute:
transport:
type: superposition
eff_diff:
type: milligton_quirk_1
......@@ -30,9 +30,9 @@ volumes:
type: const
trans_disp:
type: const
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
\ No newline at end of file
parameters:
char_length: 1.5E-11
mol_diff: 2.E-9
phi: 0.34
lambda_t: 0.0025
lambda_l: 0.025
\ No newline at end of file
......@@ -16,6 +16,70 @@
namespace Dune{
namespace Dorie{
/**
* @brief Converts PDELab grid function to match a dune-function.
* @note This adapter differs from the one provided in PDELab in that this
* one has the current interface of dune-grid to write vtk files.
* The main difference using this adapter is that it allows to write
* more than 3 components.
*
* @tparam GF PDELab grid function type.
*/
template<class GF>
class VTKGridFunctionAdapter
{
using GridFunction = typename std::decay<GF>::type;
using GridView = typename GF::Traits::GridViewType;
using Entity = typename GridView::template Codim<0>::Entity;
using Range = typename GF::Traits::RangeType;
using Domain = typename GF::Traits::DomainType;
public:
/**
* @brief Constructs the object
*
* @param[in] gf pointer to a constant grid function
*/
VTKGridFunctionAdapter(std::shared_ptr<const GF> gf)
: _gf(gf)
{}
/**
* @brief Binds an entity to the grid function
*
* @param[in] e Entity
*/
void bind(const Entity& e)
{
_entity = &e;
}
/**
* @brief Unbinds the previously binded entity
*/
void unbind()
{
_entity = nullptr;
}
/**
* @brief Evaluates the grid function at a given coordinate
*
* @param[in] x Local coordinate with respect to the last binded entity
*/
Range operator()(const Domain& x) const
{
Range y;
_gf->evaluate(*_entity,x,y);
return y;
}
private:
const Entity* _entity;
std::shared_ptr<const GF> _gf;
};
/*-------------------------------------------------------------------------*//**
* @brief Class for grid fucntion vtk sequence writer. This class works
* exacly the same way as a VTKSequence writer but it receives
......@@ -39,13 +103,20 @@ public:
* @tparam GF Type of the grid function
*/
template<class GF>
void addVertexData(const std::shared_ptr<const GF> &gf_ptr, std::string name)
void addVertexData(std::shared_ptr<const GF> gf_ptr, std::string name)
{
static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value,
"GridFunctionVTKSequenceWriter only works with only one GridView type.");
auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
Dune::VTKSequenceWriter<GV>::addVertexData(vtk_gf_ptr);
Dune::Dorie::VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
const int vtk_ncomps = GF::Traits::dimRange;
const int dim = GF::Traits::dimDomain;
auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector
: VTK::FieldInfo::Type::scalar;
auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps);
this->vtkWriter()->addVertexData(vtk_gf,vtk_info);
}
/*-----------------------------------------------------------------------*//**
......@@ -62,14 +133,21 @@ public:
static_assert(std::is_same<typename GF::Traits::GridViewType,GV>::value,
"GridFunctionVTKSequenceWriter only works with only one GridView type.");
auto vtk_gf_ptr = std::make_shared<Dune::PDELab::VTKGridFunctionAdapter<GF> >(gf_ptr, name);
Dune::VTKSequenceWriter<GV>::addCellData(vtk_gf_ptr);
Dune::Dorie::VTKGridFunctionAdapter<GF> vtk_gf(gf_ptr);
const int vtk_ncomps = GF::Traits::dimRange;
const int dim = GF::Traits::dimDomain;
auto vtk_type = (vtk_ncomps == dim) ? VTK::FieldInfo::Type::vector
: VTK::FieldInfo::Type::scalar;
auto vtk_info = VTK::FieldInfo(name,vtk_type,vtk_ncomps);
this->vtkWriter()->addCellData(vtk_gf,vtk_info);
}
/*-----------------------------------------------------------------------*//**
* @brief Clear all the attached VTKFunctions. (Necessary for each
* iteration if VTKFunctions are managed internaly with pointers
* instead of references)
* @brief Clear all the attached VTKFunctions.
* @details This is necessary for each iteration if VTKFunctions are
* managed internaly with pointers instead of references
*/
void clear()
{
......
#ifndef DUNE_DORIE_PARAMETERIZATION_READER_HH
#define DUNE_DORIE_PARAMETERIZATION_READER_HH
#include <unordered_map>
#include <memory>
#include <string>
#include <set>
#include <yaml-cpp/yaml.h>
#include <dune/common/exceptions.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune{
namespace Dorie{
/*-------------------------------------------------------------------------*//**
* @brief Interface and default reader for parameterization factories
*
* @tparam Param The parameter type to build
*/
template<class Param>
struct ParameterizationFactory
{
virtual ~ParameterizationFactory() = default;
/**
* @brief Parameterization selector
* @details Given a type_node in yaml, it selects a -polymorphic-
* parameter.
*
* @param[in] type_node YAML node describing the type of parameterization to
* choose from.
* @param[in] name Name that describes the parametrized material.
*
* @return Shared pointer to the selected parameter.
*/
virtual std::shared_ptr<Param> selector(
const YAML::Node& type_node,
const std::string name) const = 0;
/**
* @brief Parameterization reader
* @details It reads parameterizations and its parameters. For each volume,
* it assings a index to a parameterization (choosen by the
* parameterization selector) and later filled with the paramaters
* in the yaml file.
*
* @param[in] param_file YAML node with information about type and
* parameters
* @param[in] param_name Name of the parameterization.
* @param[in] log Logger
*
* @return (Unsorted) Map with parameterizations for each material.
*/
std::unordered_map<int, std::shared_ptr<Param>>
reader( const YAML::Node& param_file,
const std::string param_name,
const std::shared_ptr<spdlog::logger> log) const
{
// mapping: index on grid to parameterization object (will be returned)
std::unordered_map<int, std::shared_ptr<Param>> ret;
// set to check for duplicate indices
std::set<int> param_index_set;
// iterate over all volume nodes
for (auto&& node : param_file["volumes"])
{
const auto name = node.first.as<std::string>();
const YAML::Node param_info = node.second;
// read relevant data
const auto index = param_info["index"].as<int>();
if (param_index_set.count(index) != 0) {
log->error("Index '{}' in material '{}' of parameterization '{}'"
"already exists.",
index, name, param_name);
DUNE_THROW(IOError, "Parameter file indices must be unique");
}
param_index_set.emplace(index);
// get parameterization object
const YAML::Node& type_node = param_info[param_name];
auto parameterization = this->selector(type_node, name);
auto parameter_values = parameterization->parameters();
// fill parameterization with parameter entries
log->trace(" Reading parameters. Volume: {}, Index: {}",
name, index);
log->trace(" Required parameters: {}", parameter_values.size());
for (auto p : parameter_values)
log->trace(" {}", p.first);
auto yaml_parameters = yaml_sequence_to_map(type_node["parameters"]);
log->trace(" Available parameters: {}",
yaml_parameters.size());
unsigned param_used = 0;
for (auto&& yaml_param : yaml_parameters)
{
const auto p_name = yaml_param.first.as<std::string>();
const auto p_value = yaml_param.second.as<double>();
// try to insert the value into the parameter map
auto [begin,end] = parameter_values.equal_range(p_name);
if (begin==end) { // mark unused parameter
log->trace(" - {}: {}", p_name, p_value);
} else { // mark and set parameter
for (auto p = begin; p != end; ++p) {
p->second = p_value;
log->trace(" + {}: {}", p->first, p->second);
param_used++;
}
}
}
if (param_used!=parameter_values.size()) {
log->error(" Parameters used: {} of {}. "
"Volume: {}, Index: {}.",
param_used,parameter_values.size(),name, index);
DUNE_THROW(IOError, "Missing parameters!");
}
ret.emplace(index, parameterization);
}
return ret;
}
};
}
}
#endif // DUNE_DORIE_PARAMETERIZATION_READER_HH
\ No newline at end of file
......@@ -8,6 +8,9 @@
#include <memory>
#include <utility>
#include <yaml-cpp/yaml.h>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/gridview.hh>
......@@ -113,6 +116,61 @@ inline bool is_none (const std::string& in)
return false;
}
/*-------------------------------------------------------------------------*//**
* @brief Transform a Yaml node with sequences to map
* @details This function check all the nodes within 'in' and transform every
* sequence into a map. In particular it transform vectors and
* tensors into single key-value map pair by adding a suffix
* indicating the component.
*
* @param[in] in Map Yaml node.
*
* @return Map Yaml node with only key-double pairs.
*/
inline YAML::Node yaml_sequence_to_map(const YAML::Node& in)
{
YAML::Node out;
for (auto&& node : in)
{
const auto key = node.first.as<std::string>();
const YAML::Node& value = node.second;
if (value.IsScalar()) {
out[key] = value;
} else if (value.IsSequence()) {
if (value.size() == 1) { // scalar
out[key] = value[0].as<double>();
} else if (value.size() == 2) { // assume vector in 2D
out[key+"_x"] = value[0].as<double>();
out[key+"_y"] = value[1].as<double>();
} else if (value.size() == 3) { // assume vector in 3D
out[key+"_x"] = value[0].as<double>();
out[key+"_y"] = value[1].as<double>();
out[key+"_z"] = value[2].as<double>();
} else if (value.size() == 4) { // assume tensor in 2D
out[key+"_xx"] = value[0].as<double>();
out[key+"_xy"] = value[1].as<double>();
out[key+"_yx"] = value[2].as<double>();
out[key+"_yy"] = value[3].as<double>();
} else if (value.size() == 9) { // assume tensor in 3D
out[key+"_xx"] = value[0].as<double>();
out[key+"_xy"] = value[1].as<double>();
out[key+"_xz"] = value[2].as<double>();
out[key+"_yx"] = value[3].as<double>();
out[key+"_yy"] = value[4].as<double>();
out[key+"_yz"] = value[5].as<double>();
out[key+"_zx"] = value[6].as<double>();
out[key+"_zy"] = value[7].as<double>();
out[key+"_zz"] = value[8].as<double>();
}
} else {
DUNE_THROW(Dune::IOError,
"Yaml sequence '" << key << "' is not convertible to map pair!");
}
}
return out;
}
} // namespace Dorie
} // namespace Dune
......
......@@ -14,8 +14,7 @@
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/dorie/common/logging.hh>
#include <dune/dorie/model/richards/parameterization/interface.hh>
#include <dune/dorie/model/richards/parameterization/mualem_van_genuchten.hh>
#include <dune/dorie/model/richards/parameterization/factory.hh>
#include <dune/dorie/model/richards/parameterization/scaling.hh>
namespace Dune {
......@@ -64,21 +63,23 @@ private:
/// Index type used by the element mapper
using Index = typename Mapper::Index;
/// Base class of the parameterization
using Parameterization = Dorie::RichardsParameterization<Traits>;
using ParameterizationType = Dorie::Parameterization::Richards<Traits>;
/// Parameterization factory
using ParameterizationFactory = Dorie::Parameterization::RichardsFactory<Traits>;
/// Struct for storing scaling factors
using Scaling = Dorie::ScalingFactors<RF>;
using Scaling = Dorie::Parameterization::ScalingFactors<RF>;
/// Base class for scaling adapters
using ScaleAdapter = Dorie::ScalingAdapter<Traits>;
using ScaleAdapter = Dorie::Parameterization::ScalingAdapter<Traits>;
/// Storage for parameters and skaling factors
using ParameterStorage = std::unordered_map<Index,
std::tuple<std::shared_ptr<Parameterization>, Scaling>
std::tuple<std::shared_ptr<ParameterizationType>, Scaling>
>;
/// Need this gruesome typedef because 'map::value_type' has const key
using Cache = std::tuple<std::shared_ptr<Parameterization>, Scaling>;
using Cache = std::tuple<std::shared_ptr<ParameterizationType>, Scaling>;
using ConstCache = std::tuple<std::shared_ptr<const Parameterization>, Scaling>;
using ConstCache = std::tuple<std::shared_ptr<const ParameterizationType>, Scaling>;
/// Configuration file tree
const Dune::ParameterTree& _config;
/// Grid view of the coarsest grid configuration (level 0)
......@@ -99,7 +100,7 @@ private:
/// Check if an entity has been cached
void verify_cache () const
{
if (not std::get<std::shared_ptr<Parameterization>>(_cache)) {
if (not std::get<std::shared_ptr<ParameterizationType>>(_cache)) {
_log->error("Parameterization cache is empty. Call 'bind' before "
"accessing the cache or the parameterization");
DUNE_THROW(Dune::InvalidStateException,
......@@ -126,7 +127,7 @@ public:
for(const auto& [index, tuple] : flow_param._param) {
const auto& [p, sk] = tuple;
// make a hard copy of parameterization
std::shared_ptr<RichardsParameterization<Traits>> _p = p->clone();
std::shared_ptr<Parameterization::Richards<Traits>> _p = p->clone();
this->_param.emplace(index,std::make_tuple(_p,sk));
}
}
......@@ -197,9 +198,9 @@ public:
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
std::get<std::shared_ptr<ParameterizationType>>(_cache);
const auto cond_f = par->conductivity_f();
using Saturation = typename Parameterization::Saturation;
using Saturation = typename ParameterizationType::SaturationType;
const auto& xi_cond = std::get<Scaling>(_cache).scale_cond;
......@@ -217,9 +218,9 @@ public:
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
std::get<std::shared_ptr<ParameterizationType>>(_cache);
const auto sat_f = par->saturation_f();
using MatricHead = typename Parameterization::MatricHead;
using MatricHead = typename ParameterizationType::MatricHeadType;
const auto& scale_head = std::get<Scaling>(_cache).scale_head;
......@@ -237,13 +238,13 @@ public:
{
verify_cache();
const auto& par =
std::get<std::shared_ptr<Parameterization>>(_cache);
std::get<std::shared_ptr<ParameterizationType>>(_cache);
// get water content function and apply the scaling
const auto& scale_por = std::get<Scaling>(_cache).scale_por;
const auto wc_f = par->water_content_f(scale_por);
using Saturation = typename Parameterization::Saturation;
using Saturation = typename ParameterizationType::SaturationType;
return [wc_f](const RF saturation) {
return wc_f(Saturation{saturation}).value;
};
......@@ -259,10 +260,11 @@ private:
_log->info("Loading parameter data file: {}", param_file_name);
YAML::Node param_file = YAML::LoadFile(param_file_name);
// create the parameterization data
const auto parameterization_map = read_parameter_file(param_file);
ParameterizationFactory factory;
const auto parameterization_map = factory.reader(param_file,"richards",_log);
// build global scaling
using SAF = Dorie::ScalingAdapterFactory<Traits>;
using SAF = Dorie::Parameterization::ScalingAdapterFactory<Traits>;
auto scale_cfg = param_file["scaling"];
auto scaling_adapter = SAF::create(scale_cfg["type"].as<std::string>(),
scale_cfg["data"],
......@@ -282,105 +284,6 @@ private:
// check that mapper can map to parameterization
assert(_param.size() == _mapper.size());
}
/// Read the YAML parameter file and insert information into a map.
/** \param param_file The top YAML node to read from (file root)
*/
std::map<int, std::shared_ptr<Parameterization>> read_parameter_file (
const YAML::Node& param_file
) const
{
// mapping: index on grid to parameterization object (will be returned)
std::map<int, std::shared_ptr<Parameterization>> ret;
// set to check for duplicate indices
std::set<int> param_index_set;
// iterate over all volume nodes
for (auto&& node : param_file["volumes"])
{
const auto name = node.first.as<std::string>();
const YAML::Node param_info = node.second;
// read relevant data
const auto index = param_info["index"].as<int>();
if (param_index_set.count(index) != 0) {
_log->error("Index '{}' of parameterization '{}'"
"already exists.",
index, name);
DUNE_THROW(IOError, "Parameter file indices must be unique");
}
param_index_set.emplace(index);
// get parameterization object
const YAML::Node& type_node = param_info["water_flow"];
auto parameterization = parameterization_selector(type_node, name);
auto parameter_values = parameterization->parameters();
// fill parameterization with parameter entries
_log->trace(" Reading parameters. Volume: {}, Index: {}",
name, index);
_log->trace(" Required parameters: {}",
parameter_values.size());
for (auto p : parameter_values)
_log->trace(" {}", p.first);
const YAML::Node& yaml_parameters = param_info["parameters"];
_log->trace(" Available parameters: {}",
yaml_parameters.size());
unsigned param_used = 0;
for (auto&& yaml_param : yaml_parameters)
{
const auto p_name = yaml_param.first.as<std::string>();
const auto p_value = yaml_param.second.as<double>();
// try to insert the value into the parameter map
auto [begin,end] = parameter_values.equal_range(p_name);
std::string marker;
if (begin==end) {
// mark unused parameter
_log->trace(" - {}: {}", p_name, p_value);
} else {
for (auto p = begin; p != end; ++p) {
p->second = p_value;
_log->trace(" + {}: {}", p->first, p->second);
param_used++;
}
}
}
if (param_used!=parameter_values.size()) {
_log->error(" Parameters used: {} of {}. "
"Volume: {}, Index: {}.",
param_used,parameter_values.size(),name, index);
DUNE_THROW(IOError, "Missing parameters!");
}
ret.emplace(index, parameterization);
}
return ret;
}
/// Return a default-initialized parameterization object
/** \param type The type indicating the parameterization implementation
* \param name The name of the parameterization object (layer type)
*/
std::shared_ptr<Parameterization> parameterization_selector (
const YAML::Node& type_node,
const std::string name) const