Commit 1265ec85 authored by Santiago Ospina De Los Ríos's avatar Santiago Ospina De Los Ríos
Browse files

Merge branch '62-gradientfluxadapter-problem' into...

Merge branch '62-gradientfluxadapter-problem' into 101-couple-the-transportsimulation-with-richardssimulation
parents f6e01397 e8b7f3be
......@@ -41,7 +41,7 @@ dune_enable_all_packages()
dune_require_cxx_standard(MODULE "dorie" VERSION 14)
# add subdirectories
add_subdirectory(plugins/vendor)
add_subdirectory("plugins/vendor")
add_subdirectory("m4")
add_subdirectory("cmake/modules")
add_subdirectory("python")
......
......@@ -232,10 +232,22 @@ adding an empty line, make text **bold** or ``monospaced``.
<parameter name="fluxReconstruction">
<definition> Apply the flux reconstruction method to the solved matric
head and obtain conservative gradients. </definition>
head and obtain conservative gradients. It always computes (internally)
the local lifting independent whether the ``lifting`` keyword is active
or not.
</definition>
<suggestion> true </suggestion>
<values> true, false </values>
</parameter>
<parameter name="lifting">
<definition> Compute the local lifting for the discrete gradient. It serves
as a measure of the influence of symmetry term on the interior degrees
of freedom.
</definition>
<suggestion> false </suggestion>
<values> true, false </values>
</parameter>
</category>
<category name="NewtonParameters">
......
......@@ -21,4 +21,4 @@ USE_MATHJAX = YES
# EXCLUDE += @top_srcdir@/...
# Include documents with other extesions
FILE_PATTERNS += *.dox
FILE_PATTERNS += *.dox
\ No newline at end of file
......@@ -111,8 +111,29 @@ From the description above one can infer that one has to distinguish between *wa
Flux reconstruction
-------------------
The flux reconstruction is as projection of the fluxes used in the Discontinuous Galerkin method into a vector field function. Using correct elements, it can ensure that the fluxes in normal direction to the entity are *equivalent* to those computed by the Discontinuous Galerkin method, and most imprtantly, it can also ensure the continuity of them. This procedure makes possible to use the results of DORiE useful compute other problems that rely on the fluxes of the water (i.e. solute transport).
The flux reconstruction is a projection of the fluxes used in the Discontinuous Galerkin method into a vector field function. Using correct elements, this procedure can ensure that fluxes in normal direction to the element are *equivalent* to those computed by the Discontinuous Galerkin method, and most importantly, it can also ensure the continuity of them. Hence, the resulting vector field is useful to compute other problems that rely on the fluxes of the water (i.e. solute transport).
The flux reconstruction technique always use Raviar Thomas finite elements of one degree less than the one set for the Richards model. It can be identified in the vtk file by the name ``flux_RT{k-1}``, where ``k`` is the finite element order set for the Richards model. Flux reconstruction is not available for non-conforming grids (i.e. Cube-adaptive grids).
+---------------------------+---+---+---+
| Richards FEorder | 1 | 2 | 3 |
+============+====+=========+===+===+===+
| | 2D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | ✓ | ✓ | ✓ |
| Non-adapt. +----+---------+---+---+---+
| | 3D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | ✓ | ✓ | |
+------------+----+---------+---+---+---+
| | 2D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | | | |
| Adapt. +----+---------+---+---+---+
| | 3D | Simplex | ✓ | ✓ | ✓ |
| | +---------+---+---+---+
| | | Cube | | | |
+------------+----+---------+---+---+---+
Examples
========
......
......@@ -11,6 +11,8 @@
@{
@ingroup Common
@brief Collection of finite elements for DORiE
This group collects implementations of the different finite elements and
finite element maps not provided by dune-localfunctions and dune-pdelab
respectively.
......
......@@ -43,7 +43,9 @@ class RaviartThomasVolumeLocalFiniteElementMap
RaviartThomasVolumeLocalFiniteElement<DF,RF,k,d,gt>,d>
{
public:
/*-----------------------------------------------------------------------*//**
* @brief Constructor of RaviartThomasVolumeLocalFiniteElementMap
*/
RaviartThomasVolumeLocalFiniteElementMap() {}
#ifdef DOXYGEN /* Documentation of the code */
......
......@@ -318,7 +318,7 @@ struct SkeletonFiniteElementMap<DF,RF,dim,k,Dune::GeometryType::simplex>
template<class DF, class RF, int dim, int k>
struct SkeletonFiniteElementMap<DF,RF,dim,k,Dune::GeometryType::cube>
: public SkeletonPkFiniteElementMap<DF,RF,dim,k>
: public SkeletonQkFiniteElementMap<DF,RF,dim,k>
{};
} // namespace Dorie
......
#ifndef DUNE_DORIE_NORM_NORMAL_JUMP_SKELETON_HH
#define DUNE_DORIE_NORM_NORMAL_JUMP_SKELETON_HH
namespace Dune {
namespace Dorie {
/**
* @brief Norm of the jump in normal direction of every face in the
* skeleton.
*
* @param[in] gf A grid function of a vector space.
* @param[in] int_order The integration order
*
* @tparam GF Grid function type
*
* @return The \f$L^2\f$ norm of the jumps.
*/
template<class GF>
double norm_normal_jump_skeleton(const GF& gf, int int_order)
{
using Range = typename GF::Traits::RangeType;
constexpr int dimDomain = GF::Traits::dimDomain;
constexpr int dimRange = GF::Traits::dimRange;
static_assert(dimDomain==dimRange,
"Function 'norm_normal_jump_skeleton' is only defined for vector fields "
"with same dimension as the domain");
auto& gv = gf.getGridView();
auto& is = gv.indexSet();
double norm2;
// Traverse grid view
for (const auto& element : elements(gv))
{
// Traverse intersections
for(const auto& intersection : intersections(gv,element))
{
// Only for skeleton
if (intersection.neighbor())
{
// get entities
const auto& entity_f = intersection;
const auto& entity_i = intersection.inside();
const auto& entity_o = intersection.outside();
auto id_i = is.index(entity_i);
auto id_o = is.index(entity_o);
// Compute face only once
if (id_o > id_i)
{
// Get element geometry
auto geo_f = entity_f.geometry();
auto geo_in_i = entity_f.geometryInInside();
auto geo_in_o = entity_f.geometryInOutside();
// Loop over quadrature points and integrate normal flux
for (const auto& it : quadratureRule(geo_f,int_order))
{
const auto& position_f = it.position();
// position of quadrature point in local coordinates of elements
const auto position_i = geo_in_i.global(position_f);
const auto position_o = geo_in_o.global(position_f);
// face normal vector
const auto normal_f = entity_f.unitOuterNormal(position_f);
Range v_i, v_o;
gf.evaluate(entity_i,position_i,v_i);
gf.evaluate(entity_o,position_o,v_o);
// integration factor
auto factor = it.weight() * geo_f.integrationElement(position_f);
auto jump = normal_f*(v_i-v_o);
norm2 += jump*jump*factor;
}
}
}
}
}
return sqrt(norm2);
}
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_NORM_NORMAL_JUMP_SKELETON_HH
\ No newline at end of file
......@@ -11,6 +11,8 @@
@{
@ingroup Common
@brief Classes and engines to perform flux reconstruction
### Flux reconstruction with Raviart Thomas elements.
Following the notation of Di Petro 2012 the discrete problem for discontinuous
galerkin problems can be written as:
......
......@@ -4,19 +4,22 @@
#include <vector>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/assembler.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/local_engeine.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas/local_engine.hh>
#include <dune/dorie/common/finite_element/skeleton_fem.hh>
#include <dune/dorie/common/finite_element/raviart_thomas_volume_fem.hh>
#include <dune/dorie/common/finite_element/raviart_thomas_fem.hh>
#include <dune/dorie/common/flux_reconstruction/check_normal_jump_skeleton.hh>
namespace Dune{
namespace Dorie{
/**
* @brief Helper for the raviart thomas flux reconstrucion.
* @brief Helper for the raviart thomas flux reconstruction.
* @details This helper exports a bool that says whether the flux
* reconstruction mtehod is available for the convination of the
* reconstruction method is available for the combination of the
* used template parameters.
*
* @tparam GO Grid operator type.
......@@ -32,14 +35,14 @@ private:
static_assert(dim==2||dim==3);
static_assert(gt==GeometryType::simplex||gt==GeometryType::cube);
// hardcoded bool helpers
static constexpr bool value_2d_simplex = true;
static constexpr bool value_3d_simplex = true;
static constexpr bool value_2d_cube = (order <= 2);
static constexpr bool value_3d_cube = (order <= 1);
static constexpr bool value_simplex = (dim==2) ? value_2d_simplex : value_3d_simplex;
static constexpr bool value_cube = (dim==2) ? value_2d_cube : value_3d_cube;
static constexpr bool val_2d_simplex = true;
static constexpr bool val_3d_simplex = true;
static constexpr bool val_2d_cube = (order <= 2);
static constexpr bool val_3d_cube = (order <= 1);
static constexpr bool val_simplex = (dim==2) ? val_2d_simplex : val_3d_simplex;
static constexpr bool val_cube = (dim==2) ? val_2d_cube : val_3d_cube;
public:
static constexpr bool enable = (gt==GeometryType::simplex) ? value_simplex : value_cube;
static constexpr bool enable = (gt==GeometryType::simplex) ? val_simplex : val_cube;
};
/*-------------------------------------------------------------------------*//**
......@@ -63,10 +66,12 @@ class RaviartThomasFluxReconstruction
: public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<
typename GO::Traits::TrialGridFunctionSpace::Traits::GridViewType,
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType,
GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension,
Dune::FieldVector<
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType,
GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension
> >,
RaviartThomasFluxReconstruction<GO,gt,order>
......@@ -84,7 +89,8 @@ class RaviartThomasFluxReconstruction
enum { dim = GV::dimension };
using FEMU = Dune::Dorie::RaviartThomasLocalFiniteElementMap<GV,DF,RF,order,gt>;
using FEMU = Dune::Dorie::
RaviartThomasLocalFiniteElementMap<GV,DF,RF,order,gt>;
using ES = typename GFSW::Traits::EntitySet;
using GFSU = Dune::PDELab::GridFunctionSpace<ES,FEMU>;
using LOP = typename GO::Traits::LocalAssembler::LocalOperator;
......@@ -93,23 +99,33 @@ class RaviartThomasFluxReconstruction
using Range = Dune::PDELab::Backend::Vector<GFSU, RF>;
using LA = typename GOP::Traits::LocalAssembler;
using VolumeFEM = Dune::Dorie::RaviartThomasVolumeLocalFiniteElementMap<DF,RF,order,dim,gt>;
using VolumeFEM = Dune::Dorie::
RaviartThomasVolumeLocalFiniteElementMap<DF,RF,order,dim,gt>;
using GFSVVolume = Dune::PDELab::GridFunctionSpace<ES,VolumeFEM>;
using SkeletonFEM = typename Dune::Dorie::SkeletonFiniteElementMap<DF,RF,dim,order,gt>;
using SkeletonFEM = typename Dune::Dorie::
SkeletonFiniteElementMap<DF,RF,dim,order,gt>;
using GFSVSkeleton = Dune::PDELab::GridFunctionSpace<ES,SkeletonFEM>;
using LocalRaviartThomasEngine = Dune::Dorie::LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton>;
using LocalRaviartThomasEngine = Dune::Dorie::
LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton,true>;
// using Assember = typename GOP::Assembler;
using Assember = Dune::Dorie::DefaultAssembler<GFSU,GFSW,Dune::PDELab::EmptyTransformation,Dune::PDELab::EmptyTransformation>;
using Assember = Dune::Dorie::DefaultAssembler<GFSU,GFSW,
Dune::PDELab::EmptyTransformation,
Dune::PDELab::EmptyTransformation>;
public:
using Traits = Dune::PDELab::GridFunctionTraits<
typename GO::Traits::TrialGridFunctionSpace::Traits::GridViewType,
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
typename GO::Traits::TrialGridFunctionSpace::Traits::
FiniteElementType::Traits::LocalBasisType::
Traits::RangeFieldType,
GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension,
Dune::FieldVector<
typename GO::Traits::TrialGridFunctionSpace::Traits::FiniteElementType::Traits::LocalBasisType::Traits::RangeFieldType,
GO::Traits::TrialGridFunctionSpace::Traits::GridViewType::dimension
typename GO::Traits::TrialGridFunctionSpace::Traits::
FiniteElementType::Traits::
LocalBasisType::Traits::RangeFieldType,
GO::Traits::TrialGridFunctionSpace::Traits::
GridViewType::dimension
> >;
/*---------------------------------------------------------------------*//**
......@@ -148,7 +164,7 @@ public:
MBE mbe(0);
GOP gop(_gfsu,gfsw,local_operator,mbe);
// TODO: send our DefaultAssamber to PDELab so that it can be retrived
// TODO: send our DefaultAssamber to PDELab so that it can be retrieved
// from the grid operator:
// const auto& global_assembler_gop = gop.assembler();
......@@ -157,7 +173,8 @@ public:
const auto& local_assembler_gop = gop.localAssembler();
LocalRaviartThomasEngine local_raviart_thomas_engine(local_assembler_gop,gfsv_volume,gfsv_skeleton,_int_order_add);
LocalRaviartThomasEngine local_raviart_thomas_engine(
local_assembler_gop,gfsv_volume,gfsv_skeleton,_int_order_add);
local_raviart_thomas_engine.setPrescription(p);
local_raviart_thomas_engine.setSolution(_x);
......@@ -176,12 +193,33 @@ public:
void setTime (Time time_)
{}
//! get a reference to the GridView
/*-----------------------------------------------------------------------*//**
* @brief get a reference to the GridView.
*
* @return The grid view.
*/
inline const typename Traits::GridViewType& getGridView () const
{
return _gv;
}
/*-----------------------------------------------------------------------*//**
* @brief Check that jumps in skeleton are below some value.
* @details Tolerance is compared against the
* @f$ L^2(\Omega)\f$ norm of the normal jumps in the skeleton.
*
* @param[in] tol The tolerance.
*
* @return True if jumps are conforming w.r.t the tolerance.
*/
bool check_jumps(double tol)
{
int quadrature_factor = 2; // FIXME!
int int_order = _int_order_add+quadrature_factor*order;
auto jumps_norm = norm_normal_jump_skeleton(*this,int_order);
return Dune::FloatCmp::le(jumps_norm,tol);
}
private:
GO& _go;
GV _gv;
......@@ -190,6 +228,7 @@ private:
Range _x;
Dune::PDELab::DiscreteGridFunctionPiola<GFSU,Range> _dgfp;
const int _int_order_add;
double jumps_norm;
};
} // namespace Dorie
......
......@@ -31,15 +31,15 @@ namespace Dorie{
* operators set up properly. Following these, ideas the local
* operator must be able to handle volume local functions spaces in
* its volume methods and the equivalent for the skeleton methods.
*
*
* @warning When the order of test function is 0, the residuals are entered
* directly to the global solution vector.
* @warning This engeine always enfoce two side visit to the global
* @warning This engine always enforce two side visit to the global
* assembler.
*
*
* @remark This assembler was written following the
* Dune::PDELab::LocalAssemblerEngine interface.
*
*
* @author Santiago Ospina De Los Ríos
* @ingroup FluxReconstruction
*
......@@ -51,15 +51,19 @@ namespace Dorie{
* @todo Check parallel setup.
* @todo Bind the outside entity for skeleton gfsv instead of the inside.
* @todo Speed-up with GPUs calculating all matrices in the postAssembly step.
* @todo Check how would this work with constraits.
* @todo Check how would this work with constraints.
*
* @tparam LA The local assembler
* @tparam GFSVVolume Grid function space for the volume part of the
* ansatz space.
* @tparam GFSVSkeleton Grid function space for the skeleton part of the
* ansatz space.
* @tparam extended Decides whether to perform an extended assembling,
* that is, using the volume test function into the
* skeleton and boundary terms (useful for discrete
* gradient lifting).
*/
template<typename LA, typename GFSVVolume, typename GFSVSkeleton>
template<typename LA, typename GFSVVolume, typename GFSVSkeleton, bool extended = false>
class LocalRaviartThomasAssemblerEngine
: public Dune::PDELab::LocalAssemblerEngineBase
{
......@@ -349,7 +353,7 @@ public:
std::cout << "WARNING: Flux reconstruction for non conforming entities "
<< "produce non-continuos fluxes in normal direction of the "
<< "entity!"
<< std::endl;
<< std::endl; // TODO Change to logger
}
......@@ -382,7 +386,7 @@ public:
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_volume.size(); ++i)
r_vec[i] = rl(lfsv_volume,i);
r_vec[i] += rl(lfsv_volume,i);
auto& lfsu = lfsu_cache.localFunctionSpace();
......@@ -475,7 +479,22 @@ public:
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_skeleton.size(); ++i){
r_vec[i+offset] += weight*rl(lfsv_skeleton,i);
r_vec[i+offset] += rl(lfsv_skeleton,i);
}
if constexpr (extended) {
rl.assign(lfsw_s_cache.size(),0.0);
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaSkeleton>::
alpha_skeleton(lop,ig,
lfsw_s_cache.localFunctionSpace(),pl,lfsv_volume,
lfsw_n_cache.localFunctionSpace(),pn,lfsv_volume, //_n, //TODO add volume outside
rl_view,rn_view);
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_volume.size(); ++i){
r_vec[i] += rl(lfsv_volume,i);
}
}
auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
......@@ -559,6 +578,17 @@ public:
r_vec[i+offset] += rl(lfsv_skeleton,i);
}
if constexpr (extended) {
rl.assign(lfsw_s_cache.size(),0.0);
Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
alpha_boundary(lop,ig,lfsw_s_cache.localFunctionSpace(),pl,lfsv_volume,rl_view);
// Load residual view into the local residual vector
for (unsigned int i = 0; i < lfsv_volume.size(); ++i){
r_vec[i] += rl(lfsv_volume,i);
}
}
auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
// get polynomial degree
......
This diff is collapsed.
......@@ -22,7 +22,9 @@ RichardsSimulation<Traits>::RichardsSimulation (
Traits::dim,Traits::GridGeometryType)),
mbe_tlop(1),
time_before(0.0),
dt_before(0.0)
dt_before(0.0),
enable_fluxrc(_inifile.get<bool>("numerics.fluxReconstruction")),
enable_lifting(_inifile.get<bool>("numerics.lifting"))
{
// --- Grid Function Space ---
this->_log->trace("Setting up GridFunctionSpace");
......@@ -63,11 +65,23 @@ RichardsSimulation<Traits>::RichardsSimulation (
u = std::make_shared<U>(*gfs,.0);
Dune::PDELab::interpolate(*finitial,*gfs,*u);
// --- Operator Setup --- //
// --- Operator Setup ---
operator_setup();
update_adapters();
// --- Update Flux Reconstruction ---
if constexpr (enable_rt_engine)
if (enable_fluxrc)
try {
flux_rec_gf->update(*u);
} catch (Dune::Exception &e) {
this->_log->debug(" Flux Reconstruction failed");
this->_log->error(" Unexpected error in flux reconstruction: {}",
e.what());
DUNE_THROW(Dune::Exception, "Critical error in flux reconstruction");
}
// --- Utility Class Setup --- //
if (output_policy() != OutputPolicy::None)
{
......@@ -93,6 +107,17 @@ RichardsSimulation<Traits>::RichardsSimulation (
}
this->_log->info("Setup complete");
if ((enable_fluxrc or enable_lifting) and not enable_rt_engine) {
if (enable_fluxrc)
this->_log->error(
" Flux reconstruction is not available for this configuration.");
if (enable_lifting)
this->_log->error(
" Local lifting is not available for this configuration.");
DUNE_THROW(Dune::NotImplemented,
"Flux reconstruction engine not implemented!");
}
}
template<typename Traits>
......@@ -131,19 +156,12 @@ void RichardsSimulation<Traits>::operator_setup()
osm->setVerbosityLevel(0);
gfs->update();
if constexpr (enable_fluxrc)
{
if (inifile.get<bool>("numerics.fluxReconstruction"))
{
waterfrgf = std::make_shared<GFWaterFluxReconstruction>(*go0);
}
} else {
if (inifile.get<bool>("numerics.fluxReconstruction"))
{
DUNE_THROW(Dune::NotImplemented,
"Flux reconstruction is not implemented for the selected configuration");
}
if constexpr (enable_rt_engine) {
if (enable_fluxrc)
flux_rec_gf = std::make_shared<GFWaterFluxReconstruction>(*go0);
if (enable_lifting)
lifting_gf = std::make_shared<GFWaterFluxReconstruction>(*go0);
}
}
......@@ -216,6 +234,36 @@ void RichardsSimulation<Traits>::step()
DUNE_THROW(Dune::Exception, "Critical error in Newton solver");
}
// --- Update Flux Reconstruction ---
if constexpr (enable_rt_engine) {
if (enable_fluxrc)
try {
// compute flux reconstruction
flux_rec_gf->update(*u);
} catch (Dune::Exception &e) {
this->_log->debug(" Flux Reconstruction failed");
this->_log->error(" Unexpected error in flux reconstruction: {}",
e.what());
DUNE_THROW(Dune::Exception, "Critical error in flux reconstruction");
}
if (enable_lifting)
try {
// compute lifting
slop->penalty_term(false);
slop->symmetry_term(false);
slop->diffusion_term(false);
lifting_gf->update(*u);
slop->penalty_term(true);
slop->symmetry_term(true);
slop->diffusion_term(true);
} catch (Dune::Exception &e) {
this->_log->debug(" Local lifting failed");
this->_log->error(" Unexpected error in local lifting: {}",
e.what());
DUNE_THROW(Dune::Exception, "Critical error in local lifting");
}
}