Commit 3d862731 authored by Santiago Ospina's avatar Santiago Ospina
Browse files

[Flux Reconstruction] Added support for discrete gradient lifting

parent 6ca7b0c7
......@@ -4,7 +4,7 @@
#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>
......@@ -106,7 +106,7 @@ class RaviartThomasFluxReconstruction
SkeletonFiniteElementMap<DF,RF,dim,order,gt>;
using GFSVSkeleton = Dune::PDELab::GridFunctionSpace<ES,SkeletonFEM>;
using LocalRaviartThomasEngine = Dune::Dorie::
LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton>;
LocalRaviartThomasAssemblerEngine<LA,GFSVVolume,GFSVSkeleton,true>;
// using Assember = typename GOP::Assembler;
using Assember = Dune::Dorie::DefaultAssembler<GFSU,GFSW,
Dune::PDELab::EmptyTransformation,
......
......@@ -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
}
......@@ -478,6 +482,19 @@ public:
r_vec[i+offset] += weight*rl(lfsv_skeleton,i);
}
if constexpr (extended) {
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] += weight*rl(lfsv_volume,i);
}
}
auto& lfsu_s = lfsu_s_cache.localFunctionSpace();
// get polynomial degree
......@@ -559,6 +576,16 @@ public:
r_vec[i+offset] += rl(lfsv_skeleton,i);
}
if constexpr (extended) {
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
......
......@@ -253,8 +253,8 @@ public:
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
{
// Switches between local and global interface in test local space
using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType>;
using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
using FETest = typename LFSV::Traits::FiniteElementType;
using BasisTest = typename FETest::Traits::LocalBasisType;
/* If range dimension of the test local function space is 1 we assume that
* it is a P_k/Q_k space (probably conforming), if it is equal to world
......@@ -262,7 +262,7 @@ public:
* gradient (e.g. Raviart Thomas) which is non conforming with respect to
* the ansatz space.*/
constexpr int dimRangeTest = BasisSwitchTest::dimRange;
constexpr int dimRangeTest = BasisTest::Traits::dimRange;
const int order = lfsu.finiteElement().localBasis().order();
const int intorder = intorderadd + quadrature_factor * order;
......@@ -368,15 +368,16 @@ public:
R& r_s, R& r_n) const
{
// Switches between local and global interface in test local space
using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType>;
using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
using FETest = typename LFSV::Traits::FiniteElementType;
using BasisTest = typename FETest::Traits::LocalBasisType;
/* If domain dimension of the test local function space is dim we assume
* that it is a P_k/Q_k space (probably conforming), if it is dim-1 we
* assume that it is a skeleton finite element function space
* which is non conforming with respect to the ansatz space.*/
constexpr int dimDomainLocalTest = BasisSwitchTest::dimDomainLocal;
constexpr int dimDomainLocalTest = BasisTest::Traits::dimDomain;
constexpr int dimRangeTest = BasisTest::Traits::dimRange;
// get polynomial degree
const int order_s = lfsu_s.finiteElement().localBasis().order();
......@@ -421,8 +422,27 @@ public:
std::vector<Scalar> phiv_s(lfsv_s.size());
std::vector<Scalar> phiv_n(lfsv_n.size());
// evaluate gradient of basis function
const auto& jsu_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
const auto& jsu_n = cache[order_n].evaluateJacobian(p_local_n,lfsu_n.finiteElement().localBasis());
// transform gradients to real element
Tensor jac;
std::vector<Vector> gradphiu_s(lfsu_s.size());
std::vector<Vector> gradphiu_n(lfsu_n.size());
jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
jac.mv(jsu_s[i][0],gradphiu_s[i]);
jac = ig.outside().geometry().jacobianInverseTransposed(p_local_n);
for (unsigned int i = 0; i<lfsu_n.size(); i++)
jac.mv(jsu_n[i][0],gradphiu_n[i]);
std::vector<Vector> gradphiv_s(lfsv_s.size());
std::vector<Vector> gradphiv_n(lfsv_n.size());
// discrete gradient sign.
double dg_sign = 1.;
// evaluate basis functions of test function
if constexpr (dimDomainLocalTest == dim-1) {
// Evaluate position in local coordinates of the inside entity of codim 1
......@@ -447,11 +467,32 @@ public:
lfsv_n.finiteElement().localBasis().evaluateFunction(sub_p_local_n,phiv_n);
// discrete gradient sign.
dg_sign = 0.;
} else if constexpr (dimRangeTest == dim) {
// (we assume lfsv = RTk)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
lfsv_n.finiteElement().localBasis().evaluateFunction(p_local_n,gradphiv_n);
// Piola transformation
auto BTransposed = gtface.jacobianTransposed(it.position());
auto detB = BTransposed.determinant();
for (unsigned int i=0; i<lfsv_s.size(); i++) {
Vector y_s(gradphiv_s[i]);
y_s /= detB;
BTransposed.mtv(y_s,gradphiv_s[i]);
}
for (unsigned int i=0; i<lfsv_n.size(); i++) {
Vector y_n(gradphiv_n[i]);
y_n /= detB;
BTransposed.mtv(y_n,gradphiv_n[i]);
}
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
phiv_n = phiu_n;
gradphiv_s = gradphiu_s;
gradphiv_n = gradphiu_n;
}
// evaluate u
......@@ -461,27 +502,12 @@ public:
for (unsigned int i = 0; i<lfsu_n.size(); i++)
u_n += x_n(lfsu_n,i) * phiu_n[i];
// evaluate gradient of basis function
const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
const auto& gradphi_n = cache[order_n].evaluateJacobian(p_local_n,lfsu_n.finiteElement().localBasis());
// transform gradients to real element
Tensor jac;
std::vector<Vector> tgradphi_s(lfsu_s.size());
std::vector<Vector> tgradphi_n(lfsu_n.size());
jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
jac.mv(gradphi_s[i][0],tgradphi_s[i]);
jac = ig.outside().geometry().jacobianInverseTransposed(p_local_n);
for (unsigned int i = 0; i<lfsu_n.size(); i++)
jac.mv(gradphi_n[i][0],tgradphi_n[i]);
// compute gradient of u
Vector gradu_s(0.), gradu_n(0.);
for (unsigned int i = 0; i<lfsu_s.size(); i++)
gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
gradu_s.axpy(x_s(lfsu_s,i),gradphiu_s[i]);
for (unsigned int i = 0; i<lfsu_n.size(); i++)
gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
gradu_n.axpy(x_n(lfsu_n,i),gradphiu_n[i]);
// update with gravity vector
gradu_s -= param.gravity();
......@@ -553,55 +579,61 @@ public:
relCond_upwind = cond_upwind / satCond_upwind;
}
// update residual (flux term)
// diffusion term
// consistency term
// + penalty term
if (upwinding == RichardsDGUpwinding::none)
if constexpr (dimRangeTest != dim)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, - numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, - relCond_upwind * numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, relCond_upwind * numFlux_upwind * phiv_s[i] * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, - relCond_upwind * numFlux_upwind * phiv_n[i] * factor);
// update residual (flux term)
// diffusion term
// consistency term
// + penalty term
if (upwinding == RichardsDGUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, - dg_sign * numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * relCond_upwind * numFlux * phiv_s[i] * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, - dg_sign * relCond_upwind * numFlux * phiv_n[i] * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * relCond_upwind * numFlux_upwind * phiv_s[i] * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, - dg_sign * relCond_upwind * numFlux_upwind * phiv_n[i] * factor);
}
}
if constexpr (dimDomainLocalTest == dim-1)
return;
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
if (upwinding == RichardsDGUpwinding::none)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (tgradphi_s[i] * normal) * relCond_s * satCond_s * jump * factor);
r_s.accumulate(lfsv_s, i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_s * satCond_s * jump * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (tgradphi_n[i] * normal) * relCond_n * satCond_n * jump * factor);
r_n.accumulate(lfsv_n, i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_n * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::semiUpwind)
{
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (tgradphi_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
r_s.accumulate(lfsv_s,i, dg_sign * theta * omega_s * (gradphiv_s[i] * normal) * relCond_upwind * satCond_s * jump * factor);
for (unsigned int i = 0; i<lfsv_n.size(); i++)
r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (tgradphi_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
r_n.accumulate(lfsv_n,i, dg_sign * theta * omega_n * (gradphiv_n[i] * normal) * relCond_upwind * satCond_n * jump * factor);
}
else if (upwinding == RichardsDGUpwinding::fullUpwind)
{
for (unsigned int i = 0; i < lfsv_s.size(); i++)
r_s.accumulate(lfsv_s, i, dg_sign * theta * (tgradphi_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
r_s.accumulate(lfsv_s, i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
for (unsigned int i = 0; i < lfsv_n.size(); i++)
r_n.accumulate(lfsv_n, i, dg_sign * theta * (tgradphi_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
r_n.accumulate(lfsv_n, i, dg_sign * theta * (gradphiv_n[i] * normal) * relCond_upwind * satCond_upwind * jump * factor);
}
}
}
......@@ -616,15 +648,16 @@ public:
R& r_s) const
{
// Switches between local and global interface in test local space
using FESwitchTest = Dune::FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType>;
using BasisSwitchTest = Dune::BasisInterfaceSwitch<typename FESwitchTest::Basis>;
using FETest = typename LFSV::Traits::FiniteElementType;
using BasisTest = typename FETest::Traits::LocalBasisType;
/* If domain dimension of the test local function space is dim we assume
* that it is a P_k/Q_k space (probably conforming), if it is dim-1 we
* assume that it is a skeleton finite element function space
* which is non conforming with respect to the ansatz space.*/
constexpr int dimDomainLocalTest = BasisSwitchTest::dimDomainLocal;
constexpr int dimDomainLocalTest = BasisTest::Traits::dimDomain;
constexpr int dimRangeTest = BasisTest::Traits::dimRange;
// get polynomial degree
const int order_s = lfsu_s.finiteElement().localBasis().order();
......@@ -656,16 +689,41 @@ public:
std::vector<Scalar> phiv_s(lfsv_s.size());
// evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
const auto& jsu_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
// transform gradients to real element
const auto js = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
std::vector<Vector> gradphiu_s(lfsu_s.size());
for (unsigned int i = 0; i<lfsu_s.size(); i++)
js.mv(jsu_s[i][0],gradphiu_s[i]);
std::vector<Vector> gradphiv_s(lfsv_s.size());
// discrete gradient sign.
double dg_sign = 1.;
// evaluate basis functions of test function
if constexpr (dimDomainLocalTest == dim-1) {
// (we assume non-Galerkin method lfsu != lfsv)
lfsv_s.finiteElement().localBasis().evaluateFunction(it.position(),phiv_s);
dg_sign = 0.;
} else if constexpr (dimRangeTest == dim) {
// (we assume lfsv = RTk)
lfsv_s.finiteElement().localBasis().evaluateFunction(p_local_s,gradphiv_s);
// Piola transformation
auto BTransposed = gtface.jacobianTransposed(it.position());
auto detB = BTransposed.determinant();
for (unsigned int i=0; i<lfsv_s.size(); i++) {
Vector y_s(gradphiv_s[i]);
y_s /= detB;
BTransposed.mtv(y_s,gradphiv_s[i]);
}
dg_sign = -1.;
} else {
// (we assume Galerkin method lfsu = lfsv)
phiv_s = phiu_s;
gradphiv_s = gradphiu_s;
}
// integration factor
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
......@@ -706,26 +764,17 @@ public:
// update residual (flux term)
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, normal_flux * phiv_s[i] * factor);
r_s.accumulate(lfsv_s,i, dg_sign * normal_flux * phiv_s[i] * factor);
continue;
}
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());
// 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]);
gradu_s.axpy(x_s(lfsu_s,i),gradphiu_s[i]);
// update with gravity vector
gradu_s -= param.gravity();
......@@ -770,13 +819,13 @@ public:
// consistency term
// + penalty term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, relCond * numFlux * phiv_s[i] * factor);
r_s.accumulate(lfsv_s,i, dg_sign * relCond * numFlux * phiv_s[i] * factor);
// update residual (symmetry term)
// (non-)symmetric IP term
// symmetry term
for (unsigned int i = 0; i<lfsv_s.size(); i++)
r_s.accumulate(lfsv_s,i, dg_sign * theta * (tgradphi_s[i] * normal) * relCond * satCond_s * jump * factor);
r_s.accumulate(lfsv_s,i, dg_sign * theta * (gradphiv_s[i] * normal) * relCond * satCond_s * jump * factor);
continue;
}
......
Supports Markdown
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