[Common] Add parser to analytic initial condition

parent 1d254707
......@@ -107,23 +107,16 @@ adding an empty line, make text **bold** or ``monospaced``.
<category name="initial">
<parameter name="type">
<definition> There are currently two possible initial conditions to choose from: A HDF datafile (``data``), or linear potential gradient (``linear``).</definition>
<values> data, linear </values>
<suggestion> linear </suggestion>
<comment> Choose initial condition type: data, linear </comment>
<definition> There are currently two possible initial conditions to choose from: A HDF datafile (``data``), or analytic equations (``analytic``).</definition>
<values> data, analytic </values>
<suggestion> analytic </suggestion>
<comment> Choose initial condition type: data, analytic </comment>
<parameter name="headLower">
<definition> Initial matric head at the lower boundary of the domain </definition>
<values> float </values>
<suggestion> 0 </suggestion>
<parameter name="headGradient">
<definition> Gradient of the matric head towards the upper boundary
(positive y coordinate). </definition>
<values> float </values>
<suggestion> -1 </suggestion>
<parameter name="equation">
<definition> Equation for the initial condition </definition>
<values> equation [x,y,z] </values>
<suggestion> -y </suggestion>
<parameter name="datafile">
#include <string>
#include <plugins//vendor/exprtk/exprtk.hpp>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
......@@ -13,25 +15,27 @@ namespace Dorie {
* @brief Initial conditions for flow equation, realized by determining the
* value of the matric head
* @brief Initial conditions
* @tparam T BaseTraits
* @tparam m dimension of the Range type
template<class T, int m>
class InitialConditionLinear
class InitialConditionAnalytic
: public InitialCondition<T,m>
using RF = typename T::RF;
using Domain = typename T::Domain;
using GV = typename T::GV;
using Base = InitialCondition<T,m>;
static constexpr int dim = T::dim;
using SymbolTable = exprtk::symbol_table<RF>;
using Expression = exprtk::expression<RF>;
using Parser = exprtk::parser<RF>;
const RF _const_lower; //!< Matric head at lower boundary
const RF _const_grad; //!< Gradient of matric head for hydrostatic equilibrium
using Traits = typename Base::Traits;
......@@ -39,30 +43,45 @@ public:
/// Initial Condition Class Constructor. Reads initial condition parameters
/** \throw Dune::IOError BC File Type not supported
InitialConditionLinear (const GV& gv, RF const_lower, RF const_grad)
InitialConditionAnalytic (const GV& gv, const std::string& equation_str)
: Base(gv)
, _const_lower(const_lower)
, _const_grad(const_grad)
, _equation_str(equation_str)
SymbolTable symbol_table;
if constexpr (dim == 3)
Parser parser;
if (!parser.compile(_equation_str,_expression))
DUNE_THROW(IOError,"Unexpected error with math expression!");
/// Evaluate the initial condition at certain position
/** \param y Value of matric head at given position
/** \param y
* \param x Position in local element coordinates
* \param e Element enitity
* \param xGlobal Position in global coordinates
void evaluate ( const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const override
constexpr int dim = GV::dimension;
auto xGlobal = e.geometry().global(x);
y = _const_lower+_const_grad*xGlobal[dim-1];
_x_global = e.geometry().global(x);
y = _expression.value();
mutable Domain _x_global;
Expression _expression;
const std::string _equation_str;
......@@ -2,7 +2,7 @@
#include <dune/dorie/common/initial_condition/linear.hh>
#include <dune/dorie/common/initial_condition/analytic.hh>
namespace Dune{
namespace Dorie{
......@@ -21,7 +21,6 @@ struct InitialConditionFactory
static auto create( const Dune::ParameterTree& ini_file,
const typename T::GV& grid_view)
using RF = typename T::RF;
std::unique_ptr<InitialCondition<T,m>> ic;
auto ic_type = ini_file.get<std::string>("initial.type");
......@@ -39,12 +38,11 @@ struct InitialConditionFactory
"Data file type '." << file_type << "' not supported!");
} else if (ic_type == "linear") {
auto head_low = ini_file.get<RF>("initial.headLower");
auto head_grad = ini_file.get<RF>("initial.headGradient");
} else if (ic_type == "analytic") {
auto equation_str = ini_file.get<std::string>("initial.equation");
using ICLinear = InitialConditionLinear<T,m>;
ic = std::make_unique<ICLinear>(grid_view,head_low,head_grad);
using ICAnalytic = InitialConditionAnalytic<T,m>;
ic = std::make_unique<ICAnalytic>(grid_view,equation_str);
} else {
// log->error("Unknown initial condition type: {}",ic_type);
......@@ -5,4 +5,3 @@ set(DUNE_REENABLE_ADD_TEST True)
add_subdirectory(spdlog EXCLUDE_FROM_ALL)
add_subdirectory(exprtk EXCLUDE_FROM_ALL)
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