...
  View open merge request
Commits (70)
......@@ -9,16 +9,16 @@ set(DORIE_MAX_RORDER_2 6)
set(DORIE_MAX_RORDER_3 6)
# Maximum polynomial orders of transport model for available targets
set(DORIE_MAX_TORDER_2 0)
set(DORIE_MAX_TORDER_3 0)
set(DORIE_MAX_TORDER_2 6)
set(DORIE_MAX_TORDER_3 6)
# Maximum polynomial orders of Richards model for default targets
set(DORIE_MAX_DEFAULT_RORDER_2 3)
set(DORIE_MAX_DEFAULT_RORDER_3 1)
# Maximum polynomial orders of transport model for default targets
set(DORIE_MAX_DEFAULT_TORDER_2 0)
set(DORIE_MAX_DEFAULT_TORDER_3 0)
set(DORIE_MAX_DEFAULT_TORDER_2 1)
set(DORIE_MAX_DEFAULT_TORDER_3 1)
#
# .. cmake_function:: dorie_compile_instance
......
......@@ -229,6 +229,13 @@ adding an empty line, make text **bold** or ``monospaced``.
</category>
<category name="numerics">
<parameter name="penaltyFactor">
<definition> Penalty factor to be used in the Discontinuous Galerkin scheme
</definition>
<values> float </values>
<suggestion> 10 </suggestion>
</parameter>
<parameter name="timestepMethod">
<definition> Numerical scheme to perform time steps in the simulation.
``alex2`` and ``implicit_euler`` are implicit methods.
......@@ -302,4 +309,39 @@ adding an empty line, make text **bold** or ``monospaced``.
<values> float </values>
</parameter>
</category>
<category name="numerics.experimental">
<parameter name="method">
<definition> DG discretization method for skeleton terms.
**SIPG:** Symmetric Interior Penalty
**NIPG:** Non-Symmetric Interior Penalty
**OOB:** Oden, Babuska, Baumann: no penalty term
**IIP:** Incomplete Interior Penalty: no symmetry term
</definition>
<values> SIPG, NIPG, OOB, IIP </values>
<suggestion> SIPG </suggestion>
<comment> Experimental settings are enabled by the appropriate CMake flag.
</comment>
</parameter>
<parameter name="upwinding">
<definition> Upwinding method for skeleton terms.
**fullUpwind:** Apply upwinding on numeric flux.
</definition>
<values> none, semiUpwind, fullUpwind </values>
<suggestion> none </suggestion>
</parameter>
<parameter name="weights">
<definition> Apply harmonic weighting to skeleton term contributions.
</definition>
<values> true, false </values>
<suggestion> true </suggestion>
</parameter>
</category>
</dorie>
......@@ -9,8 +9,11 @@
// Checks for custom setup
#ifndef DORIE_TORDER
#error DORIE_TORDER macro must be defined for Transport factory builds
#elif DORIE_TORDER != 0
#error "Transport FV solver only supports DORIE_TORDER=0"
#else
static_assert(DORIE_TORDER >= 0,
"Negative polynomial orders do not make sense");
static_assert(DORIE_TORDER <= 6,
"DORiE only supports polynomial orders up to 6");
#endif
#ifdef DORIE_RT_ORDER
......
......@@ -15,7 +15,6 @@
@defgroup BuildProcess Model Build Process
@{
@ingroup Models
@brief How DORiE libraries and executables are built and linked
## Overview
......
#ifndef DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#define DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#ifndef DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
#define DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
#include <iostream>
#include <string>
......@@ -22,24 +22,25 @@
namespace Dune{
namespace Dorie{
/// Adaptivity base class. Does nothing.
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U, int order>
class AdaptivityHandlerBase
class RichardsAdaptivityHandler
{
private:
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
using RF = typename Traits::RF;
public:
AdaptivityHandlerBase () = default;
virtual ~AdaptivityHandlerBase () = default;
RichardsAdaptivityHandler () = default;
virtual ~RichardsAdaptivityHandler () = default;
/// Do nothing.
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& u) { }
/// Do nothing.
virtual void mark_grid (Grid& grid, GV& gv, GFS& gfs, const Param& param, const Boundary& boundary, const RF time, U& u) { }
/// Do nothing.
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
/// Do nothing.
virtual void adapt_grid (Grid& grid, GFS& gfs, U& u) { }
};
/// Specialized AdaptivityHandlerBase class providing functions and members for adaptive grid refinement
......@@ -49,7 +50,7 @@ template<typename Traits,
typename Boundary,
typename U,
int order>
class AdaptivityHandler : public AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>
class WaterFluxAdaptivityHandler : public RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>
{
private:
using RF = typename Traits::RF;
......@@ -61,7 +62,7 @@ private:
/// Adaptivity GFS
using AGFS = typename AGFSHelper::Type;
/// Error estimator local operator
using ESTLOP = Dune::Dorie::Operator::FluxErrorEstimator<Traits,Param,Boundary>;
using ESTLOP = Dune::Dorie::Operator::WaterFluxErrorEstimator<Traits,Param,Boundary>;
/// Empty constraints container
using NoTrafo = Dune::PDELab::EmptyTransformation;
/// Solution vector type
......@@ -90,8 +91,8 @@ public:
* \param _inifile Parameter file parser
* \param grid Grid to adapt (reference is not saved)
*/
AdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>(),
WaterFluxAdaptivityHandler (Dune::ParameterTree& _inifile, Grid& grid)
: RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>(),
inifile(_inifile),
maxLevel(inifile.get<int>("adaptivity.maxLevel")),
minLevel(inifile.get<int>("adaptivity.minLevel")),
......@@ -220,7 +221,7 @@ public:
* define which create() function is enabled at compile time.
*/
template<typename Traits, typename GFS, typename Param, typename Boundary, typename U, int order>
struct AdaptivityHandlerFactory
struct RichardsAdaptivityHandlerFactory
{
private:
......@@ -246,7 +247,7 @@ private:
using Grid = typename Traits::Grid;
using AHB = AdaptivityHandlerBase<Traits,GFS,Param,Boundary,U,order>;
using AHB = RichardsAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>;
Dune::ParameterTree& inifile;
Grid& grid;
......@@ -257,17 +258,17 @@ public:
/*---------------------------------------------------------------------*//**
* @brief Create the factory, taking the parameters for building an
* AdaptivityHandler
* RichardsAdaptivityHandler
*
* @param _inifile Parameter file parser
* @param _grid Grid to adapt (reference is not saved)
*/
AdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
RichardsAdaptivityHandlerFactory (Dune::ParameterTree& _inifile, Grid& _grid) :
inifile(_inifile),
grid(_grid)
{ }
/// Create a dummy AdaptivityHandler
/// Create a dummy RichardsAdaptivityHandler
/** Assert that adaptive grid (UG) is used.
*/
template<bool active = enabled>
......@@ -276,15 +277,15 @@ public:
return std::make_unique<AHB>();
}
/// Create a true AdaptivityHandler
/// Create a true RichardsAdaptivityHandler
template<bool active = enabled>
std::enable_if_t<active,std::unique_ptr<AHB>> create ()
{
std::unique_ptr<AHB> p;
if (inifile.hasSub("adaptivity") &&
inifile.get<std::string>("adaptivity.policy") != "none")
p = std::make_unique<AdaptivityHandler<Traits,GFS,Param,Boundary,U,order>>(inifile,grid);
else
inifile.get<std::string>("adaptivity.policy") == "waterFlux")
p = std::make_unique<WaterFluxAdaptivityHandler<Traits,GFS,Param,Boundary,U,order>>(inifile,grid);
else
p = std::make_unique<AHB>();
return p;
}
......@@ -293,5 +294,4 @@ public:
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_ADAPTIVITY_HANDLER_HH
#endif // DUNE_DORIE_RICHARDS_ADAPTIVITY_HANDLER_HH
......@@ -27,9 +27,12 @@ namespace Dune {
* A call to residual() of a grid operator space will assemble
* the quantity \eta_T for each cell.
* Skeleton integral is mostly the same as in the DG operator.
*
* @ingroup LocalOperators
* @ingroup TransportModel
*/
template<typename Traits, typename Parameter, typename Boundary>
class FluxErrorEstimator
class WaterFluxErrorEstimator
: public Dune::PDELab::LocalOperatorDefaultFlags
{
......@@ -62,7 +65,7 @@ namespace Dune {
RF penalty_factor;
RF time;
FluxErrorEstimator (const Dune::ParameterTree& config,
WaterFluxErrorEstimator (const Dune::ParameterTree& config,
const Parameter& param_, const Boundary& boundary_,
RichardsDGMethod method_ = RichardsDGMethod::SIPG,
int intorderadd_ = 2, int quadrature_factor_ = 2)
......
......@@ -47,7 +47,7 @@ enum class RichardsUpwinding
/* \return Tuple of RichardsUpwinding, RichardsDGMethod, and
* RichardsDGWeights to be inserted into the local operator constructor.
*/
inline auto read_operator_settings(
inline auto read_richards_operator_settings(
const Dune::ParameterTree& inifile)
-> std::tuple<Operator::RichardsUpwinding,
Operator::RichardsDGMethod,
......
......@@ -62,7 +62,7 @@ ModelRichards<Traits>::ModelRichards (
// Read local operator settings
namespace OP = Dune::Dorie::Operator;
const auto settings = OP::read_operator_settings(inifile);
const auto settings = OP::read_richards_operator_settings(inifile);
// --- Local Operators ---
if constexpr (order>0)
......
......@@ -185,9 +185,9 @@ struct ModelRichardsTraits : public BaseTraits
/// Custom VTK output writer
using Writer = Dune::Dorie::GridFunctionVTKSequenceWriter<GV>;
/// Grid Adaptivity Handler base class
using AdaptivityHandler = Dune::Dorie::AdaptivityHandlerBase<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
using AdaptivityHandler = Dune::Dorie::RichardsAdaptivityHandler<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
/// Factory for creating the AdaptivityHandler
using AdaptivityHandlerFactory = Dune::Dorie::AdaptivityHandlerFactory<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
using AdaptivityHandlerFactory = Dune::Dorie::RichardsAdaptivityHandlerFactory<BaseTraits,GFS,FlowParameters,FlowBoundary,U,order>;
};
/*-------------------------------------------------------------------------*//**
......
This diff is collapsed.
This diff is collapsed.
......@@ -103,10 +103,6 @@ private:
using WaterFlux = typename GFWaterFlux::Traits::RangeType;
using WaterContent = typename GFWaterContent::Traits::RangeType;
// Extract world dimension
enum {dim = GFWaterContent::Traits::dimDomain};
using Tensor = Dune::FieldMatrix<RF,dim,dim>;
public:
/**
......@@ -116,7 +112,6 @@ public:
* @param[in] boundary The boundary terms
* @param[in] gf_water_flux The grid function of water flux
* @param[in] gf_water_content The grid function of water content
* @param[in] dirichlet_mode The dirichlet mode
*/
TransportFVSpatialOperator(
std::shared_ptr<const Parameter> param,
......@@ -522,9 +517,9 @@ public:
public:
TransportFVTemporalOperator()
: _time(0.)
{ }
TransportFVTemporalOperator()
: _time(0.)
{ }
/*---------------------------------------------------------------------*//**
* @brief Volume integral depending on test and ansatz functions
......
#ifndef DUNE_DORIE_TRANSPORT_LOP_CFG_HH
#define DUNE_DORIE_TRANSPORT_LOP_CFG_HH
#include <tuple>
#include <string>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/dorie/common/logging.hh>
namespace Dune {
namespace Dorie {
namespace Operator {
/**
* \addtogroup TransportModel
* \{
*/
//! Interior penalty type
enum class TransportDGMethod
{
NIPG, //!< Non-symmetric interior penalty
SIPG, //!< Symmetric interior penalty (default)
OBB, //!< Oden, Babuska, Baumann (no penalty term)
IIP //!< Incomplete interior penalty (no symmetry term)
};
//! Gradient weighting
enum class TransportDGWeights
{
weightsOn, //!< harmonic weighting of flux contributions
weightsOff //!< no weighting of flux contributions
};
//! Upwinding type
enum class TransportUpwinding
{
fullUpwind, //!< Upwind conductivity
none //!< No upwinding
};
/// Read operator settings from a parameter file
/* \return Tuple of TransportUpwinding, TransportDGMethod, and
* TransportDGWeights to be inserted into the local operator constructor.
*/
inline auto read_transport_operator_settings(
const Dune::ParameterTree& inifile)
-> std::tuple<Operator::TransportUpwinding,
Operator::TransportDGMethod,
Operator::TransportDGWeights>
{
const auto log = Dorie::get_logger(log_transport);
log->debug("Reading local operator settings:");
TransportUpwinding upwinding;
TransportDGMethod method = TransportDGMethod::SIPG;
TransportDGWeights weights = TransportDGWeights::weightsOn;
// Upwinding
const auto upwinding_str
= inifile.get<std::string>("numerics.upwinding");
log->debug(" Upwinding scheme: {}", upwinding_str);
if (upwinding_str == "none")
upwinding = TransportUpwinding::none;
else if (upwinding_str == "fullUpwind")
upwinding = TransportUpwinding::fullUpwind;
else {
log->error("Unknown upwinding scheme: {}", upwinding_str);
DUNE_THROW(Dune::IOError, "Unknown upwinding scheme");
}
// Return here if the local operator is FV only
if (inifile.get<int>("numerics.FEorder") == 0)
{
log->debug(" Ignoring settings 'numerics.DGMethod' and "
"'numerics.DGWeights' for finite volume solver.");
return std::make_tuple(upwinding, method, weights);
}
// DG Method
const auto method_str = inifile.get<std::string>("numerics.DGMethod");
log->debug(" DG method: {}", method_str);
if (method_str == "SIPG")
method = TransportDGMethod::SIPG;
else if (method_str == "NIPG")
method = TransportDGMethod::NIPG;
else if (method_str == "OBB")
method = TransportDGMethod::OBB;
else if (method_str == "IIP")
method = TransportDGMethod::IIP;
else {
log->error("Unknown DG method: {}", method_str);
DUNE_THROW(Dune::IOError, "Unknown DG method");
}
// Harmonic weighting of fluxes
const auto weights_str = inifile.get<bool>("numerics.DGWeights");
log->debug(" Harmonic weights for interface fluxes (DG): {}",
weights_str);
if (weights_str)
weights = TransportDGWeights::weightsOn;
else
weights = TransportDGWeights::weightsOff;
return std::make_tuple(upwinding, method, weights);
}
/**
* \}
* // group TransportModel
*/
} // namespace Operator
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_TRANSPORT_LOP_CFG_HH
......@@ -59,6 +59,31 @@ ModelTransport<Traits>::ModelTransport(
this->_log);
sinitial = SoluteInitialFactory::create(inifile, gv, this->_log);
// Read local operator settings
namespace OP = Dune::Dorie::Operator;
const auto settings = OP::read_transport_operator_settings(inifile);
// --- Local Operators ---
if constexpr (order>0)
{
this->_log->debug("Setting up local grid operators: DG method");
namespace OP = Dune::Dorie::Operator;
const auto method = std::get<OP::TransportDGMethod>(settings);
const auto upwinding = std::get<OP::TransportUpwinding>(settings);
const auto weights = std::get<OP::TransportDGWeights>(settings);
slop = std::make_unique<SLOP>(inifile, sparam, sboundary,
method, upwinding,
weights);
}
else {
this->_log->debug("Setting up local grid operators: FV method");
slop = std::make_unique<SLOP>(sparam, sboundary/*, upwinding*/);
}
tlop = std::make_unique<TLOP>();
// --- Solution Vectors and Initial Condition ---
const RF ini_value = inifile.get<RF>("initial.value",0.);
u = std::make_shared<U>(*gfs,ini_value);
......@@ -114,27 +139,31 @@ void ModelTransport<Traits>::operator_setup()
dof_count = gv.comm().sum(dof_count);
this->_log->debug(" Total number of DOF: {}", dof_count);
// --- Local Operators ---
slop = std::make_unique<SLOP>(sparam, sboundary);
tlop = std::make_unique<TLOP>();
Dune::PDELab::constraints(*this->gfs,*this->cc,false);
// --- Grid Operators ---
go0 = std::make_shared<GO0>(*gfs,*cc,*gfs,*cc,*slop,mbe_slop);
go1 = std::make_unique<GO1>(*gfs,*cc,*gfs,*cc,*tlop,mbe_tlop);
igo = std::make_unique<IGO>(*go0,*go1);
if (ts_param->implicit())
implicit_igo = std::make_unique<ImplicitIGO>(*go0,*go1);
else
explicit_igo = std::make_unique<ExplicitIGO>(*go0,*go1);
// --- Solvers ---
ls = std::make_unique<LS>(0);
double reduction = inifile.get<double>("solverParameters.reduction");
double min_defect = inifile.get<double>("solverParameters.minDefect");
pdesolver = std::make_unique<PDESOLVER>(*igo,*ls,reduction,min_defect,0);
lsgfs = std::make_unique<LSGFS>(LSGFSHelper::create(gv));
lscc = std::make_unique<LSCC>();
if (ts_param->implicit())
implicit_ls = std::make_unique<ImplicitLS>(*implicit_igo,*cc,*lsgfs,*lscc,1000,0,true,true);
else
explicit_ls = std::make_unique<ExplicitLS>(*explicit_igo,*cc,*lsgfs,*lscc,1000,0,true,true);
// --- Time Step Operators ---
osm = std::make_unique<OSM>(*ts_param,*igo,*pdesolver);
osm->setVerbosityLevel(0);
if (ts_param->implicit()){
implicit_slps = std::make_unique<ImplicitSLPS>(*implicit_igo,*implicit_ls,1e-10);
implicit_osm = std::make_unique<ImplicitOSM>(*ts_param,*implicit_igo,*implicit_slps);
implicit_osm->setVerbosityLevel(0);
} else {
explicit_osm = std::make_unique<ExplicitOSM>(*ts_param,*explicit_igo,*explicit_ls);
explicit_osm->setVerbosityLevel(0);
}
gfs->update();
}
......@@ -189,7 +218,10 @@ void ModelTransport<Traits>::step()
// create a DOF vector for next step
dwarn.push(false);
osm->apply(time, dt, *u, *unext);
if (ts_param->implicit())
implicit_osm->apply(time, dt, *u, *unext);
else
explicit_osm->apply(time, dt, *u, *unext);
dwarn.pop();
// short time step log for level 'info' (after success)
......@@ -201,6 +233,17 @@ void ModelTransport<Traits>::step()
time+dt);
}
if (ts_param->implicit()) {
auto result = implicit_osm->result();
this->_log->trace(" Matrix assembly: {:.2e}s (total), {:.2e}s (success)",
result.total.assembler_time,
result.successful.assembler_time);
this->_log->trace(" Linear solver: {:.2e}s (total), {:.2e}s (success)",
result.total.linear_solver_time,
result.successful.linear_solver_time);
}
u = unext;
// update time and time step
......@@ -296,5 +339,38 @@ void ModelTransport<Traits>::write_data () const
}
}
template<typename Traits>
void ModelTransport<Traits>::mark_grid()
{
if constexpr (AdaptivityHandlerFactory::enabled)
{
if (adaptivity_policy() != AdaptivityPolicy::SoluteFlux)
DUNE_THROW(Dune::NotImplemented,
"Not known adaptivity policy for transport solver");
adaptivity->mark_grid(*grid, gv, *gfs, *sparam, *sboundary,
water_flux_gf, water_content_gf, time_before+dt_before, *u);
} else
DUNE_THROW(NotImplemented,
"Dorie does not implement adaptivity for this simulation object!");
}
template<typename Traits>
void ModelTransport<Traits>::adapt_grid()
{
if constexpr (AdaptivityHandlerFactory::enabled)
adaptivity->adapt_grid(*grid, *gfs, *u);
else
DUNE_THROW(NotImplemented,
"Dorie does not implement adaptivity for this simulation object!");
}
template<typename Traits>
void ModelTransport<Traits>::post_adapt_grid()
{
operator_setup();
}
} // namespace Dorie
} // namespace Dune
This diff is collapsed.
......@@ -53,3 +53,5 @@ time.maxTimestep = 1E6
time.startTimestep = 1E4
numerics.timestepMethod = implicit_euler
numerics.FEorder = 0
\ No newline at end of file