Commit ed160bfa authored by Lukas Riedel's avatar Lukas Riedel

Add Dirichlet and Outflow BCs with muparser expressions

parent 233bd1b9
#ifndef DUNE_DORIE_COMMON_BOUNDARY_CONDITION_DIRICHLET_HH
#define DUNE_DORIE_COMMON_BOUNDARY_CONDITION_DIRICHLET_HH
#include <array>
#include <muParser.h>
#include <dune/dorie/common/boundary_condition/boundary_condition.hh>
namespace Dune {
......@@ -81,6 +85,97 @@ public:
}
};
/// Dirichlet boundary condition with muParser expression
/**
* \ingroup Boundary
*/
template<typename RF>
class DirichletExprBoundaryCondition : public DirichletBoundaryCondition<RF>
{
protected:
using Base = DirichletBoundaryCondition<RF>;
using Time = typename Base::Time;
const std::string _expr;
/// The parser instance
mutable mu::Parser _parser;
/// Position cache
mutable std::array<RF, 3> _pos{0.0, 0.0, 0.0};
mutable RF _time;
public:
/// Construct a boundary condition
DirichletExprBoundaryCondition (
const std::string_view name,
const TimeInterval<Time> interval,
const std::string_view expr,
const SoluteConcentration concentration=SoluteConcentration::none)
:
Base(name, interval, 0.0),
_expr(expr)
{
// Initialize parser
_parser.DefineConst("pi",
Dune::StandardMathematicalConstants<double>::pi());
_parser.DefineVar("x", &_pos[0]);
_parser.DefineVar("y", &_pos[1]);
_parser.DefineVar("z", &_pos[2]);
_parser.DefineVar("t", &_time);
_parser.DefineVar("h", &_pos[1]);
// set up parser expression
try
{
_parser.SetExpr(_expr);
// try to evaluate once
_parser.Eval();
}
catch (mu::Parser::exception_type& e)
{
handle_parser_error(e);
}
}
template<typename Vector>
RF evaluate (const Time time, const Vector& pos) const
{
if (not this->in_interval(time)) {
DUNE_THROW(InvalidStateException, "Evaluating boundary condition "
"with invalid time");
}
_time = time;
std::copy(pos.begin(), pos.end(), _pos.begin());
if (pos.size() == 3) {
_parser.DefineVar("h", &_pos[2]);
}
RF ret;
try {
ret = _parser.Eval();
}
catch (mu::Parser::exception_type& e) {
handle_parser_error(e);
}
return ret;
}
/// Destroy this object
virtual ~DirichletExprBoundaryCondition () override = default;
/// Return the type of this boundary condition
virtual Operator::BCType type () const override
{
return Operator::BCType::DirichletExpr;
}
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "dirichletexpr";
};
} // namespace Dorie
} // namespace Dune
......
......@@ -55,10 +55,12 @@ private:
};
/// List of know type strings for boundary conditions
static constexpr inline std::array<std::string_view, 3> known_type_strs = {
static constexpr inline std::array<std::string_view, 5> known_type_strs = {
DirichletBoundaryCondition<RF>::type_str,
NeumannBoundaryCondition<RF>::type_str,
OutflowBoundaryCondition<RF>::type_str
OutflowBoundaryCondition<RF>::type_str,
DirichletExprBoundaryCondition<RF>::type_str,
OutflowExprBoundaryCondition<RF>::type_str
};
/// Return the type strings of permitted boundary conditions
......@@ -287,6 +289,20 @@ public:
// read the time interval from the config
const auto t_int = get_time_interval(config, time_end_default);
// Dirichlet BC with muParser expression
if (type == DirichletExprBoundaryCondition<RF>::type_str) {
const YAML::Node& node = config["value"];
return std::make_shared<DirichletExprBoundaryCondition<RF>>(
name, t_int, node.as<std::string>(), SoluteConcentration::none
);
}
else if (type == OutflowExprBoundaryCondition<RF>::type_str) {
const YAML::Node& node = config["value"];
return std::make_shared<OutflowExprBoundaryCondition<RF>>(
name, t_int, node.as<std::string>(), SoluteConcentration::none
);
}
// read the values from the config
const auto [v_begin, v_end] = get_values(config);
......
......@@ -3,6 +3,8 @@
#include <dune/dorie/common/boundary_condition/boundary_condition.hh>
#include <dune/dorie/common/boundary_condition/dirichlet.hh>
namespace Dune {
namespace Dorie {
......@@ -32,7 +34,36 @@ public:
}
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "outflow";
static constexpr std::string_view type_str = "outflow";
};
/// Class representing an Outflow BC with muparser expression
/** The exact implementation depends on the model.
*
* \ingroup Boundary
*/
template<typename RF>
class OutflowExprBoundaryCondition : public DirichletExprBoundaryCondition<RF>
{
protected:
using Base = DirichletExprBoundaryCondition<RF>;
using Time = typename Base::Time;
public:
/// Re-use base constructor
using Base::Base;
/// Destroy this object
virtual ~OutflowExprBoundaryCondition () override = default;
/// Return the type of this boundary condition
virtual Operator::BCType type () const override
{
return Operator::BCType::OutflowExpr;
}
/// The type string of this boundary condition class
static constexpr std::string_view type_str = "outflowexpr";
};
} // namespace Dorie
......
......@@ -77,7 +77,7 @@ public:
// try to evaluate once
_parser.Eval();
} catch (mu::Parser::exception_type& e) {
handle_parser_error(e);
handle_parser_error(e, _log);
}
}
......@@ -101,7 +101,7 @@ public:
y = _parser.Eval();
}
catch (mu::Parser::exception_type& e) {
handle_parser_error(e);
handle_parser_error(e, _log);
}
}
......@@ -112,22 +112,6 @@ private:
mutable Domain _pos_global;
/// The parser instance
mu::Parser _parser;
/// Output information on the parser error and throw DUNE exception
/**
* \param e Exception thrown by the parser
* \throw IOError (always throws)
*/
void handle_parser_error (const mu::Parser::exception_type& e) const
{
_log->error("Evaluating analytic initial condition failed:");
_log->error(" Parsed expression: {}",e.GetExpr());
_log->error(" Token: {}",e.GetToken());
_log->error(" Error position: {}",e.GetPos());
_log->error(" Error code: {}",e.GetCode());
_log->error(" Error message: {}",e.GetMsg());
DUNE_THROW(IOError, "Error evaluating analytic initial condition");
}
};
}
......
......@@ -17,8 +17,11 @@ namespace Operator {
enum class BCType
{
Neumann, //!< Fixed flux at boundary
NeumannExpr,
Dirichlet, //!< Fixed matric head at boundary
Outflow, //!< Free flow of solute through boundary
DirichletExpr,
Outflow, //!< Free flow of solute through boundary
OutflowExpr,
none //!< No BC specified
};
......
......@@ -16,12 +16,16 @@
#include <yaml-cpp/yaml.h>
#include <muParser.h>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune {
namespace Dorie {
......@@ -226,6 +230,25 @@ inline YAML::Node yaml_sequence_to_map(const YAML::Node& in)
return out;
}
/// Output information on a muparser error and throw DUNE exception
/**
* \param e Exception thrown by the parser
* \throw IOError (always throws)
*/
inline void
handle_parser_error(
const mu::Parser::exception_type& e,
const std::shared_ptr<spdlog::logger> log = get_logger(log_base))
{
log->error("Evaluating analytic initial condition failed:");
log->error(" Parsed expression: {}", e.GetExpr());
log->error(" Token: {}", e.GetToken());
log->error(" Error position: {}", e.GetPos());
log->error(" Error code: {}", e.GetCode());
log->error(" Error message: {}", e.GetMsg());
DUNE_THROW(IOError, "Error evaluating analytic initial condition");
}
} // namespace Dorie
} // namespace Dune
......
......@@ -674,7 +674,10 @@ public:
continue;
}
else if (bc_type == BCType::Dirichlet or bc_type == BCType::Outflow)
else if (bc_type == BCType::Dirichlet
or bc_type == BCType::DirichletExpr
or bc_type == BCType::Outflow
or bc_type == BCType::OutflowExpr)
{
// compute gradient of u
Vector gradu_s(0.);
......@@ -688,7 +691,16 @@ public:
const Domain normal = ig.unitOuterNormal(it.position());
// jump relative to Dirichlet value
const RF g = bc->evaluate(time);
RF g = bc->evaluate(time);
if (bc_type == BCType::DirichletExpr
or bc_type == BCType::OutflowExpr)
{
auto& bc_expr =
static_cast<const DirichletExprBoundaryCondition<RF>&>(*bc);
g = bc_expr.evaluate(time, gtface.global(it.position()));
}
const RF jump = u_s - g;
// conductivity factors
......@@ -720,8 +732,9 @@ public:
}
// Avoid accumulating inflow for outflow BC
if (bc_type == BCType::Outflow and
FloatCmp::lt(double(numFlux), 0.0)) {
if ((bc_type == BCType::Outflow or bc_type == BCType::OutflowExpr)
and FloatCmp::lt(double(numFlux), 0.0))
{
continue;
}
......
......@@ -264,7 +264,8 @@ public:
const auto conductivity_f_i = _param->conductivity_f();
const auto bc_type = bc->type();
if (bc_type == BCType::Dirichlet or bc_type == BCType::Outflow)
if (bc_type == BCType::Dirichlet or bc_type == BCType::DirichletExpr
or bc_type == BCType::Outflow or bc_type == BCType::OutflowExpr)
{
const auto geo_i = entity_i.geometry();
const auto center_position_i_g = geo_i.center();
......@@ -277,7 +278,14 @@ public:
const auto distance = position_diff.two_norm();
// Evaluate Dirichlet condition
const auto g = bc->evaluate(_time);
auto g = bc->evaluate(_time);
if (bc_type == BCType::DirichletExpr or bc_type == BCType::OutflowExpr)
{
auto& bc_expr =
static_cast<const DirichletExprBoundaryCondition<RangeFieldU>&>(*bc);
g = bc_expr.evaluate(_time, center_position_f_g);
}
// Inside unknown value
const RangeU u_i = x_i(lfsu_i,0);
......@@ -295,8 +303,9 @@ public:
const auto water_flux_n = - cond_i*dudn;
// Avoid accumulation of inflow for outflow BC
if (bc_type == BCType::Outflow and
FloatCmp::lt(double(water_flux_n), 0.0)) {
if ((bc_type == BCType::Outflow or bc_type == BCType::OutflowExpr)
and FloatCmp::lt(double(water_flux_n), 0.0))
{
return;
}
......
......@@ -126,6 +126,28 @@ TEST_F(BCFactory, Outflow)
EXPECT_NO_THROW(_bcf_transport->create("name", config, 1.0));
}
TEST_F(BCFactory, DirichletExpr)
{
YAML::Node config = YAML::Load("{type: DirichletExpr,"
"value: x+y,"
"time: 0.0"
"}");
EXPECT_NO_THROW(_bcf_none->create("name", config, 1.0));
EXPECT_NO_THROW(_bcf_richards->create("name", config, 1.0));
EXPECT_NO_THROW(_bcf_transport->create("name", config, 1.0));
}
TEST_F(BCFactory, OutflowExpr)
{
YAML::Node config = YAML::Load("{type: OutflowExpr,"
"value: x+y,"
"time: 0.0"
"}");
EXPECT_NO_THROW(_bcf_none->create("name", config, 1.0));
EXPECT_NO_THROW(_bcf_richards->create("name", config, 1.0));
EXPECT_NO_THROW(_bcf_transport->create("name", config, 1.0));
}
/* ------------------------------------------------------------------------- */
// Test BC manager and loaded BC data.
// This verifies that YAML data is properly interpreted
......
......@@ -6,6 +6,7 @@
#include <type_traits>
#include <string>
#include <string_view>
#include <dune/common/exceptions.hh>
......@@ -163,3 +164,59 @@ TEST(Outflow, StaticInformation)
EXPECT_EQ(bc.type(), Operator::BCType::Outflow);
EXPECT_EQ(bc.type_str, "outflow");
}
/* ------------------------------------------------------------------------- */
// Test Dirichlet BC with muParser expression
template<template<class> class ExprBC>
ExprBC<RF>
create_expr_bc(std::string_view expr)
{
Dune::Dorie::create_base_logger();
return ExprBC<RF>(
"bc", { 0.0, 1.0 }, expr, Dune::Dorie::SoluteConcentration::none);
}
TEST(DirichletExpr, StaticInformation)
{
const auto bc =
create_expr_bc<Dune::Dorie::DirichletExprBoundaryCondition>("0");
EXPECT_EQ(bc.type(), Dune::Dorie::Operator::BCType::DirichletExpr);
EXPECT_EQ(bc.type_str, "dirichletexpr");
}
TEST(DirichletExpr, Evaluate2D)
{
const auto bc = create_expr_bc<Dune::Dorie::DirichletExprBoundaryCondition>("x+2*y-t");
std::array<double, 2> pos{0.0, 0.0};
EXPECT_DOUBLE_EQ(bc.evaluate(0.0, pos), 0.0);
pos = {1.0, 0.0};
EXPECT_DOUBLE_EQ(bc.evaluate(0.0, pos), 1.0);
pos = {1.0, 1.0};
EXPECT_DOUBLE_EQ(bc.evaluate(0.0, pos), 3.0);
EXPECT_DOUBLE_EQ(bc.evaluate(1.0, pos), 2.0);
}
TEST(DirichletExpr, Evaluate3D)
{
const auto bc = create_expr_bc<Dune::Dorie::DirichletExprBoundaryCondition>("x+2*y+3*h-t");
std::array<double, 3> pos{0.0, 0.0, 0.0};
EXPECT_DOUBLE_EQ(bc.evaluate(0.0, pos), 0.0);
pos = {1.0, 0.0, 0.0};
EXPECT_DOUBLE_EQ(bc.evaluate(0.0, pos), 1.0);
pos = {1.0, 1.0, 0.0};
EXPECT_DOUBLE_EQ(bc.evaluate(0.0, pos), 3.0);
pos = {1.0, 1.0, 1.0};
EXPECT_DOUBLE_EQ(bc.evaluate(0.0, pos), 6.0);
EXPECT_DOUBLE_EQ(bc.evaluate(1.0, pos), 5.0);
}
TEST(OutflowExpr, StaticInformation)
{
const auto bc =
create_expr_bc<Dune::Dorie::OutflowExprBoundaryCondition>("0");
EXPECT_EQ(bc.type(), Dune::Dorie::Operator::BCType::OutflowExpr);
EXPECT_EQ(bc.type_str, "outflowexpr");
}
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