Skip to content

GitLab

Commit c3fd7c01 authored by Lukas Riedel's avatar Lukas Riedel
Browse files
Options

Merge branch '124-remove-boundary-conditions-evaporation-and-limitedinflux' into 'master' 

Resolve "Remove boundary conditions "evaporation" and "limitedInflux""

Closes #124

See merge request !120
parents fbb7f2fa dfacfe26
Hide whitespace changes
Inline Side-by-side
Showing with 46 additions and 230 deletions
CHANGELOG.md
View file @ c3fd7c01
......@@ -134,7 +134,7 @@
### Removed
* The class `H5File::AttributeReader` was completely removed.
* Boundary conditions `evaporation` and `limitedInflux` !120
## 1.1.1 (2018-08-21)
......
doc/examples/bc_files/2d_evap.bcdat deleted 100644 → 0
View file @ fbb7f2fa
spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 1
0 evaporation 3.55e-8 -6.0 dirichlet 0
\ No newline at end of file
doc/examples/bc_files/3d_evap.bcdat deleted 100644 → 0
View file @ fbb7f2fa
spatial_resolution_north_we 0
spatial_resolution_north_fb 0
spatial_resolution_south_we 0
spatial_resolution_south_fb 0
spatial_resolution_west_sn -1
spatial_resolution_west_fb -1
spatial_resolution_east_sn -1
spatial_resolution_east_fb -1
spatial_resolution_front_sn -1
spatial_resolution_front_we -1
spatial_resolution_back_sn -1
spatial_resolution_back_we -1
number_BC_change_times 1
0 evaporation 3.55e-8 -6.0 dirichlet 0
\ No newline at end of file
doc/man-bcfile.rst
View file @ c3fd7c01
......@@ -30,20 +30,6 @@ Neumann
Set a fixed volumetric flux [m/s] at the boundary. Positive values imply fluxes out of the domain, negative values imply fluxes into the domain.
Limited Influx
++++++++++++++
.. object:: limited_influx <flux>
Set a volumetric flux [m/s] at the boundary (Neumann BC). Fluxes into the domain are stopped as soon as the matric head at the respective boundary segment reaches :math:`h_m \geq 0.0 \, \text{m}`.
Evaporation
+++++++++++
.. object:: evaporation <flux> <head>
Set a volumetric flux [m/s] (Neumann BC) and a fallback matric head [m] (Dirichlet BC) at the boundary. The Neumann BC is invoked as long as the water content inside the soil can sustain the flux. Before every time step, the numeric flux caused by the given Dirichlet BC is estimated. If it is lower than the given Neumann flux, the Dirichlet BC is applied.
Datafile Structure
==================
......@@ -96,9 +82,8 @@ These lines follow a simple grammar:
.. productionlist:: bc
bc_line: time { group }*
time: `float`
group: ( bc_type1 `float` ) | ( bc_type2 `float` `float` )
bc_type1: "neumann" | "dirichlet" | "limited_influx"
bc_type2: "evaporation"
group: ( bc_type `float` )
bc_type: "neumann" | "dirichlet"
The boundary conditions defined here are parsed in the same order as the boundary segments have been specified. In 3D, the rectangular boundary segments are parsed in a tabular fashion, where columns run faster than rows. Columns are defined along the first direction specified in the `Boundary Segmentation`_, and rows are defined along the second direction.
......
dune/dorie/common/boundary_condition.hh
View file @ c3fd7c01
......@@ -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++;
}
......
dune/dorie/common/typedefs.hh
View file @ c3fd7c01
......@@ -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)
{
......
dune/dorie/model/richards/local_operator.hh
View file @ c3fd7c01
......@@ -60,7 +60,7 @@ struct RichardsDGWeights
struct BCType
{
enum Type { Neumann, //!< Fixed flux at boundary
dirichlet, //!< Fixed matric head at boundary
Dirichlet, //!< Fixed matric head at boundary
none //!< No BC specified
};
};
......@@ -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;
......@@ -590,84 +583,7 @@ public:
else if (BoundaryCondition::isDirichlet(bc))
{
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);
}
bcType = BCType::Dirichlet;
}
else if (BoundaryCondition::isOther(bc))
......@@ -687,7 +603,7 @@ public:
continue;
}
else if (bcType == BCType::dirichlet)
else if (bcType == BCType::Dirichlet)
{
// 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());
......@@ -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
......
dune/dorie/model/richards/richards.cc
View file @ c3fd7c01
......@@ -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);
......
test/bcs/evaporation_2d.dat deleted 100644 → 0
View file @ fbb7f2fa
spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 1
0 evaporation 2e-7 -100 dirichlet 0
test/bcs/infiltration_evaporation_2d.dat deleted 100644 → 0
View file @ fbb7f2fa
spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 3
0 neumann -5.55e-6 dirichlet 0
1e5 neumann 0 dirichlet 0
1e5 evaporation 1e-7 -10 dirichlet 0
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