Commit a8a230c4 authored by Lukas Riedel's avatar Lukas Riedel

Remove BC options "evaporation" and "limitedInflux"

* remove options from input scheme
* remove options in Richards LOP
* remove grid view from LOP constructor arguments
* remove BC type cache from LOP
parent fbb7f2fa
......@@ -120,7 +120,7 @@ namespace Dorie{
/*-----------------------------------------------------------------------*//**
* @brief Read necessary parameters. Then read out BC data file
* with rectangular boundary decomposition.
*
*
* @throws Dune::IOError BC File path not specified
*
* @param[in] config The configuration
......@@ -643,19 +643,6 @@ namespace Dorie{
bcData[i][j].headValue.push_back(stod(buffer));
bcData[i][j].fluxValue.push_back(0.0);
}
else if(buffer == "limited_influx"){
bcData[i][j].type.push_back(BoundaryCondition::LimitedInFlux);
instream >> buffer;
bcData[i][j].fluxValue.push_back(stod(buffer));
bcData[i][j].headValue.push_back(0.0);
}
else if(buffer == "evaporation"){
bcData[i][j].type.push_back(BoundaryCondition::Evaporation);
instream >> buffer;
bcData[i][j].fluxValue.push_back(stod(buffer));
instream >> buffer;
bcData[i][j].headValue.push_back(stod(buffer));
}
else
DUNE_THROW(IOError,"In BC file line " + std::to_string(counter) + ": BC side "
+ std::to_string(i) + ": BC declaration " + std::to_string(j+1) + ": BC type unknown!");
......@@ -821,11 +808,6 @@ namespace Dorie{
std::cout << "neumann : ";
else if(bcData[i][j].type[k]==BoundaryCondition::Dirichlet)
std::cout << "dirichlet : ";
else if(bcData[i][j].type[k]==BoundaryCondition::LimitedInFlux)
std::cout << "limited_influx : ";
else if(bcData[i][j].type[k]==BoundaryCondition::Evaporation)
std::cout << "evaporation : ";
std::cout << bcData[i][j].fluxValue[k] << " , " << bcData[i][j].headValue[k] << std::endl;
}
index++;
}
......
......@@ -13,8 +13,6 @@ namespace Dune{
enum Type {
Neumann, //!< Fixed flux at boundary
Dirichlet, //!< Fixed matric head at boundary
Evaporation, //!< Special BC which switches between Neumann and Dirichlet depending on the water content
LimitedInFlux, //!< Special Neumann BC which stops infiltration into saturated medium
Other //!< UNUSED
};
......@@ -28,16 +26,6 @@ namespace Dune{
{
return (i == Neumann);
}
//! Test for Limited Flux boundary condition
static bool isLimitedInFlux (Type i)
{
return (i == LimitedInFlux);
}
//! Test for Evaporation boundary condition
static bool isEvaporation (Type i)
{
return (i == Evaporation);
}
//! UNUSED: Test for other(?) boundary condition
static bool isOther (Type i)
{
......
......@@ -132,7 +132,7 @@ inline auto read_experimental_operator_settings(
* to allow for different orders of the shape functions, which should be
* enough to support p-adaptivity. (Another likely candidate would be
* differing geometry types, i.e. hybrid meshes.)
*
*
* @ingroup LocalOperators
*/
template<typename Traits, typename Parameter, typename Boundary, typename SourceTerm, typename FEM, bool adjoint>
......@@ -148,7 +148,6 @@ template<typename Traits, typename Parameter, typename Boundary, typename Source
protected:
enum { dim = Traits::dim };
typedef typename Traits::GridView GridView;
typedef typename Traits::RangeField RF;
typedef typename Traits::Scalar Scalar;
typedef typename Traits::Vector Vector;
......@@ -161,7 +160,6 @@ protected:
typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
typedef Dune::PDELab::LocalBasisCache<LocalBasisType> Cache;
using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
const Parameter& param; //!< class providing equation parameters
const Boundary& boundary; //!< class providing boundary conditions
......@@ -169,15 +167,12 @@ protected:
const RichardsDGMethod::Type method; //!< Interior penalty type
const RichardsDGUpwinding::Type upwinding; //!< DG method weights switch
const RichardsDGWeights::Type weights; //!< Gradient weighting
const Mapper mapper; //!< maps grid elements (cells) to natural numbers
const int intorderadd; //!< integration order addend
const int quadrature_factor; //!< factor to FEM integration order
RF penalty_factor; //!< Interior penalty factor
RF theta; //!< Symmetry factor
RF time; //!< operator internal time for BC, parameter, and source term queries
//! Cache for fixing the BC type over one entire time step
mutable std::map<typename std::size_t,BoundaryCondition::Type> bc_type_cache;
std::vector<Cache> cache; //!< Local basis cache
......@@ -196,7 +191,6 @@ public:
/// Create operator.
/**
* \param param_ Object providing equation parameters
* \param view_ Grid view
* \param boundary_ Object prividing boundary conditions
* \param sourceTerm_ Object providing source term information
* \param method_ DG method type \see RichardsDGMethod::Type
......@@ -205,18 +199,28 @@ public:
* \param intorderadd_ Addend to integration order
* \param quadrature_factor_ Factor to FEM integration order
*/
RichardsDGSpatialOperator (const Dune::ParameterTree& config, const Parameter& param_, const GridView view_,
Boundary& boundary_, const SourceTerm& sourceTerm_,
RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::none,
RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
int intorderadd_ = 2, int quadrature_factor_ = 2)
: Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
RichardsDGSpatialOperator (
const Dune::ParameterTree& config,
const Parameter& param_,
const Boundary& boundary_,
const SourceTerm& sourceTerm_,
const RichardsDGMethod::Type method_ = RichardsDGMethod::SIPG,
const RichardsDGUpwinding::Type upwinding_ = RichardsDGUpwinding::none,
const RichardsDGWeights::Type weights_ = RichardsDGWeights::weightsOn,
int intorderadd_ = 2,
int quadrature_factor_ = 2
) :
Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), upwinding(upwinding_), weights(weights_),
mapper(view_, Dune::mcmgElementLayout()),
intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
param(param_),
boundary(boundary_),
sourceTerm(sourceTerm_),
method(method_),
upwinding(upwinding_),
weights(weights_),
intorderadd(intorderadd_),
quadrature_factor(quadrature_factor_),
penalty_factor(config.get<RF>("numerics.penaltyFactor")),
time(0.0),
cache(20)
......@@ -226,7 +230,7 @@ public:
theta = 1.;
else if (method == RichardsDGMethod::SIPG)
theta = -1.;
if (method == RichardsDGMethod::OOB)
penalty_factor = 0.;
}
......@@ -547,13 +551,8 @@ public:
const RF satCond_s = conductivity_f_s(1.0);
// loop over quadrature points and integrate normal flux
std::size_t i = 0;
for (const auto& it : quadratureRule(gtface,intorder))
{
// calculate unique id
const auto it_id = create_unique_id(ig.intersection(),i);
i++;
// position of quadrature point in local coordinates of elements
const auto p_local_s = ig.geometryInInside().global(it.position());
......@@ -568,15 +567,9 @@ public:
for (unsigned int i = 0; i<lfsu_s.size(); i++)
u_s += x_s(lfsu_s,i) * phi_s[i];
// check if a BC type is cached
BoundaryCondition::Type bc;
const auto bc_cached = bc_type_cache.find(it_id);
if(bc_cached!=bc_type_cache.end())
bc = bc_cached->second;
else // query type of boundary condition
bc = boundary.bc(ig.intersection(),it.position(),time);
// query the type of boundary condition
// set the internal bcType flag
const auto bc = boundary.bc(ig.intersection(),it.position(),time);
auto bcType = BCType::none;
RF normal_flux;
......@@ -593,83 +586,6 @@ public:
bcType = BCType::dirichlet;
}
else if (BoundaryCondition::isLimitedInFlux(bc))
{
// modified Neumann boundary condition, set flux to 0 where potential is positive
normal_flux = .0;
if(u_s < .0) normal_flux = boundary.j(ig.intersection(),it.position(),time);
bcType = BCType::Neumann;
}
else if (BoundaryCondition::isEvaporation(bc))
{
// estimate the dirichlet flux
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
// transform gradients to real element
const auto jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
std::vector<Vector> tgradphi_s(lfsu_s.size());
for (unsigned int i = 0; i<lfsu_s.size(); i++)
jac.mv(gradphi_s[i][0],tgradphi_s[i]);
// compute gradient of u
Vector gradu_s(0.);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
// update with gravity vector
gradu_s -= param.gravity();
// face normal vector
const Domain normal = ig.unitOuterNormal(it.position());
// jump relative to Dirichlet value
const RF g = boundary.g(ig.intersection(),it.position(),time);
const RF jump = u_s - g;
// conductivity factors
const RF relCond_s =
conductivity_f_s(saturation_f_s(u_s)) / satCond_s;
const RF relCond_n =
conductivity_f_s(saturation_f_s(g)) / satCond_s;
// penalty
const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
// compute numerical flux
const RF numFlux = satCond_s * (-(gradu_s * normal) + penalty * jump);
// compute conductivity factor
RF relCond;
if (upwinding == RichardsDGUpwinding::semiUpwind
|| upwinding == RichardsDGUpwinding::fullUpwind)
{
if (numFlux > 0)
relCond = relCond_s;
else
relCond = relCond_n;
}
else // (upwinding == RichardsDGUpwinding::none)
{
// harmonic average of conductivity factors
relCond = 2 * relCond_s * relCond_n / (relCond_s + relCond_n);
}
// switches between Neumann and Dirichlet bc
// this choice is kept constant for one time step
normal_flux = boundary.j(ig.intersection(),it.position(),time);
if(relCond * numFlux < normal_flux) {
bcType = BCType::dirichlet;
bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
}
else {
bcType = BCType::Neumann;
bc_type_cache.emplace(it_id,BoundaryCondition::Neumann);
}
}
else if (BoundaryCondition::isOther(bc))
bcType = BCType::none;
......@@ -838,39 +754,12 @@ public:
void preStep (RF t, RF dt, int stages)
{
time = t;
bc_type_cache.clear();
}
private:
//! create a unique ID for caching the BC type at single quadrature points
/** Use Cantor's pairing function
* \param it Intersection object
* \param qp_id Index of the quadrature point
*/
template<typename Intersection>
size_t create_unique_id (const Intersection& it, const std::size_t qp_id) const
{
const auto id_entity = mapper.index(it.inside());
const auto id_inter = it.indexInInside();
const std::size_t id = cantor_pairing(id_entity,id_inter);
return (std::size_t)cantor_pairing(id,qp_id);
}
//! Cantor's Pairing Function: Map two natural numbers to one
template<typename T, typename U>
decltype(auto) cantor_pairing (const T val1, const U val2) const
{
static_assert(std::is_integral<T>::value,"argument must be integral");
static_assert(std::is_integral<U>::value,"argument must be integral");
return (val1+val2)*(val1+val2+1)/2 +val2;
}
};
/**
* @brief DG local operator for the temporal derivative of the Richards equation
*
*
* @ingroup LocalOperators
*/
template<typename Traits, typename Parameter, typename FEM, bool adjoint>
......@@ -911,11 +800,17 @@ template<typename Traits, typename Parameter, typename FEM, bool adjoint>
/**
* @brief Constructor
*/
RichardsDGTemporalOperator (const Dune::ParameterTree& config, const Parameter& param_, const typename Traits::GridView view,
int intorderadd_ = 0, int quadrature_factor_ = 2)
: Dune::PDELab::NumericalJacobianVolume<RichardsDGTemporalOperator<Traits, Parameter, FEM, adjoint> >(1.e-7),
param(param_), intorderadd(intorderadd_), quadrature_factor(quadrature_factor_)
{}
RichardsDGTemporalOperator (
const Dune::ParameterTree& config,
const Parameter& param_,
int intorderadd_ = 0,
int quadrature_factor_ = 2
) :
Dune::PDELab::NumericalJacobianVolume<RichardsDGTemporalOperator<Traits, Parameter, FEM, adjoint> >(1.e-7),
param(param_),
intorderadd(intorderadd_),
quadrature_factor(quadrature_factor_)
{ }
/**
* @brief Volume integral depending on test and ansatz functions
......
......@@ -49,13 +49,13 @@ RichardsSimulation<Traits>::RichardsSimulation (
const auto upwinding = std::get<OP::RichardsDGUpwinding::Type>(settings);
const auto weights = std::get<OP::RichardsDGWeights::Type>(settings);
slop = std::make_unique<SLOP>(inifile, *fparam, gv, *fboundary, *fsource,
slop = std::make_unique<SLOP>(inifile, *fparam, *fboundary, *fsource,
method, upwinding, weights);
#else
slop = std::make_unique<SLOP>(inifile, *fparam, gv, *fboundary, *fsource);
slop = std::make_unique<SLOP>(inifile, *fparam, *fboundary, *fsource);
#endif // EXPERIMENTAL_DG_FEATURES
tlop = std::make_unique<TLOP>(inifile, *fparam, gv);
tlop = std::make_unique<TLOP>(inifile, *fparam);
controller = std::make_unique<CalculationController>(
inifile, *fboundary, this->_log);
......
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