Commit 36db13b0 authored by Santiago Ospina's avatar Santiago Ospina

more on lifting: Global lifting solved with DG (does not converge) -> TODO:...

more on lifting: Global lifting solved with DG (does not converge) -> TODO: Calculate it locally [ci skip]
parent 8ce31d12
......@@ -211,36 +211,61 @@ void RichardsSimulation<Traits>::update_adapters ()
template<typename Traits>
void RichardsSimulation<Traits>::write_data (RF time)
{
constexpr int l = 0; // FIXME! this can be 0 or 1
if constexpr (Traits::write_output)
{
update_adapters();
std::shared_ptr<Dune::Dorie::FluxReconstruction<typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
1,
Traits::GridGeometryType>> lift = std::make_shared<
Dune::Dorie::FluxReconstruction<typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
1,
Traits::GridGeometryType>>(udgf);
constexpr bool frimp = FluxReconstructionImplemented<typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
Traits::fem_order - l,
Traits::GridGeometryType>::value;
if (inifile.get<bool>("output.vertexData")) {
vtkwriter->addVertexData(udgf,"head");
vtkwriter->addVertexData(fluxdgf,"flux");
vtkwriter->addVertexData(condgf,"K_0");
vtkwriter->addVertexData(waterdgf,"theta_w");
vtkwriter->addVertexData(satdgf,"Theta");
// vtkwriter->addVertexData(lift,"lift");
vtkwriter->template addVertexData<GFMatricHead>(udgf,"head");
vtkwriter->template addVertexData<GFWaterFlux>(fluxdgf,"flux");
vtkwriter->template addVertexData<GFConductivity>(condgf,"K_0");
vtkwriter->template addVertexData<GFWaterContent>(waterdgf,"theta_w");
vtkwriter->template addVertexData<GFSaturation>(satdgf,"Theta");
if constexpr (frimp)
{
using GFFluxReconstruction = Dune::Dorie::FluxReconstruction<
typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
Traits::fem_order - l,
Traits::GridGeometryType>;
std::shared_ptr<GFFluxReconstruction> frgf = std::make_shared<GFFluxReconstruction>(udgf);
vtkwriter->template addVertexData<GFFluxReconstruction>(frgf,"flux (reconstruction)");
} else {
std::cout << "WARNING: Flux reconstruction not available for this simulation!" << std::endl;
}
} else {
vtkwriter->addCellData(udgf,"head");
vtkwriter->addCellData(fluxdgf,"flux");
vtkwriter->addCellData(condgf,"K_0");
vtkwriter->addCellData(waterdgf,"theta_w");
vtkwriter->addCellData(satdgf,"Theta");
// vtkwriter->addCellData(lift,"lift");
vtkwriter->template addCellData<GFMatricHead>(udgf,"head");
vtkwriter->template addCellData<GFWaterFlux>(fluxdgf,"flux");
vtkwriter->template addCellData<GFConductivity>(condgf,"K_0");
vtkwriter->template addCellData<GFWaterContent>(waterdgf,"theta_w");
vtkwriter->template addCellData<GFSaturation>(satdgf,"Theta");
if constexpr (frimp)
{
using GFFluxReconstruction = Dune::Dorie::FluxReconstruction<
typename Traits::GV,
typename Traits::RangeField,
Traits::GV::dimension,
Traits::fem_order - l,
Traits::GridGeometryType>;
std::shared_ptr<GFFluxReconstruction> frgf = std::make_shared<GFFluxReconstruction>(udgf);
vtkwriter->template addCellData<GFFluxReconstruction>(frgf,"flux (reconstruction)");
} else {
std::cout << "WARNING: Flux reconstruction not available for this simulation!" << std::endl;
}
}
try{
vtkwriter->write(time,output_type);
}
......
......@@ -34,15 +34,15 @@ namespace Dune{
* @param gfs The GridFunctionsSpace
* @param x The coefficients vector
*/
DiscreteGridFunctionGradient(std::shared_ptr<GFS> gfs, std::shared_ptr<X> x)
DiscreteGridFunctionGradient(const std::shared_ptr<GFS> gfs, std::shared_ptr<X> x)
: Base(*gfs,*x)
, pgfs(gfs)
, px(x)
{}
private:
std::shared_ptr<GFS> pgfs;
std::shared_ptr<X> px;
const std::shared_ptr<GFS> pgfs;
const std::shared_ptr<X> px;
};
}
......
......@@ -3,6 +3,12 @@
#define DUNE_DORIE_FLUX_RECONSTRUCTION_HH
#include <dune/pdelab/finiteelementmap/raviartthomasfem.hh>
#include <dune/pdelab/stationary/linearproblem.hh>
#include <dune/pdelab/gridfunctionspace/subspace.hh>
// #include <dune/typetree/compositenode.hh>
// #include <dune/typetree/utility.hh>
namespace Dune {
namespace Dorie {
......@@ -12,7 +18,7 @@ namespace Dorie {
*
* @tparam <unnamed> { description }
*/
template<class GF>
template<class FEM, class GF>
class GlobalLifting
: public Dune::PDELab::LocalOperatorDefaultFlags
{
......@@ -21,8 +27,101 @@ public:
: _grid_function(grid_function)
{}
enum {doAlphaSkeleton = true};
enum {doSkeletonTwoSided = true}; // FIXME, my programer was too lazy to do the two sides at once
enum {doPatternVolume = true };
enum {doPatternSkeleton = true };
enum {doAlphaSkeleton = true};
enum {doAlphaVolume = true};
// jacobian of volume term
template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
void jacobian_volume (const EG& eg, const LFSU& lfsu,
const X& x, const LFSV& lfsv,
M& mat) const
{
const unsigned int size = LFSU::CHILDREN;
using namespace Indices;
const auto& lfsv_base = child(lfsv,_0);
const auto& lfsu_base = child(lfsu,_0);
// const int order = lfsu.finiteElement().localBasis().order();
const int intorder = 5; //FIXME
// get element geomentry
auto gt = eg.geometry();
// loop over quadrature points
for (const auto& it : quadratureRule(gt,intorder))
{
// evaluate position in element local and global coordinates
const auto p_local = it.position();
// evaluate basis function (assume power lfs and galerkin method)
auto& phi = cache.evaluateFunction(p_local,lfsu_base.finiteElement().localBasis());
auto& psi = cache.evaluateFunction(p_local,lfsv_base.finiteElement().localBasis());
// integration factor
const double factor = it.weight() * gt.integrationElement(p_local);
for (unsigned int i=0; i<size; i++)
for (unsigned int j=0; j<size; j++)
for (unsigned int k = 0; k<lfsv.child(i).size(); k++)
for (unsigned int l = 0; l<lfsv.child(j).size(); l++)
mat.accumulate(lfsu.child(i),k,lfsv.child(j), l, phi[k] * psi[l] * factor);
}
}
template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x,
const LFSV& lfsv, R& r) const
{
const unsigned int size = LFSU::CHILDREN;
using namespace Indices;
using LFSU_BASE = TypeTree::Child<LFSU,_0>;
const auto& lfsv_base = child(lfsv,_0);
const auto& lfsu_base = child(lfsu,_0);
using RF = typename LFSU_BASE::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
using Range = Dune::FieldVector<RF,size>;
// const int order = lfsu.finiteElement().localBasis().order();
const int intorder = 5; //FIXME
// get element geomentry
auto gt = eg.geometry();
// loop over quadrature points
for (const auto& it : quadratureRule(gt,intorder))
{
// evaluate position in element local and global coordinates
const auto p_local = it.position();
// evaluate basis function (assume power lfs and galerkin method)
auto& phi = cache.evaluateFunction(p_local,lfsu_base.finiteElement().localBasis());
auto& psi = cache.evaluateFunction(p_local,lfsv_base.finiteElement().localBasis());
// evaluate u
Range u;
for(unsigned int d=0; d<size; ++d) {
const auto& lfsu_i = lfsu.child(d);
for(unsigned int i=0; i<size; i++)
u[d] += x(lfsu_i,i) * phi[i];
}
// integration factor
const double factor = it.weight() * gt.integrationElement(p_local);
// update residual
for (unsigned int i = 0; i<size; i++)
for (unsigned int j = 0; j<lfsv.child(i).size(); j++)
r.accumulate(lfsv.child(i),j, u[i] * psi[j] * factor);
}
}
/*-----------------------------------------------------------------------*//**
* @brief Skeleton integral depending on test and ansatz functions. Each
......@@ -51,44 +150,315 @@ public:
const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o,
R& r_i, R& r_o) const
{
const unsigned int size = LFSU::CHILDREN;
using namespace Indices;
using LFSU_BASE = TypeTree::Child<LFSU,_0>;
const auto& lfsu_base_i = child(lfsu_i,_0);
const auto& lfsu_base_o = child(lfsu_o,_0);
using RF = typename LFSU::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
using RF = typename LFSU_BASE::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
// const int order = lfsu.finiteElement().localBasis().order();
const int intorder = 5; //FIXME
const int intorder = 1; //FIXME
auto face_geo = ig.geometry();
for (const auto& it : quadratureRule(face_geo,intorder)) {
assert(lfsu_i.size()==lfsv_i.size());
assert(lfsu_o.size()==lfsv_o.size());
const auto p_local_i = ig.geometryInInside().global(it.position());
const auto p_local_o = ig.geometryInOutside().global(it.position());
auto phi_i = cache.evaluateFunction(p_local_i,lfsu_base_i.finiteElement().localBasis());
auto phi_o = cache.evaluateFunction(p_local_o,lfsu_base_o.finiteElement().localBasis());
int n = lfsu_i.size();
Dune::DynamicMatrix<RF> A(n,n);
Dune::DynamicVector<RF> b(n);
const auto normal = ig.unitOuterNormal(it.position());
// loop over quadrature points and integrate normal flux
for (const auto& it : quadratureRule(face_geo,intorder)) {
for (std::size_t i = 0; i < lfsv_i.size(); ++i) {
for (std::size_t j = 0; j < lfsv_o.size(); ++j) {
if (lfsv_i.localIndex(i) == lfsv_o.localIndex(j)) {
const auto average = 0.5*(x_i(lfsv_i,i) + x_o(lfsv_o,j)); // FIXME: average of a vector field!
const auto normal = ig.unitOuterNormal(it.position());
typename GF::Traits::RangeType flux;
typename GF::Traits::DomainType x_local = ig.geometryInInside().global(it.position());
_grid_function->evaluate(ig.inside(),x_local,flux);
b[i] += average * flux * normal ; // FIXME eval theta = jump in flux
}
}
for (std::size_t j = 0; j < lfsu_i.size(); ++j)
A[j][i] = x_i(lfsu_i,j)*x_i(lfsv_i,i); // Check order of i and j access
typename GF::Traits::RangeType gf_i, gf_o;
_grid_function->evaluate(ig.inside(),p_local_i,gf_i);
_grid_function->evaluate(ig.outside(),p_local_o,gf_o);
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
// update residual
for (unsigned int i = 0; i<size; i++)
{
for (unsigned int j = 0; j<lfsv_i.child(i).size(); j++)
r_i.accumulate(lfsv_i.child(i),j, (gf_i - gf_o) * normal[i] * phi_i[j] * 0.5 * factor);
for (unsigned int j = 0; j<lfsv_o.child(i).size(); j++)
r_o.accumulate(lfsv_o.child(i),j, (gf_i - gf_o) * normal[i] * phi_o[j] * 0.5 * factor);
}
}
Dune::DynamicVector<RF> x(n);
A.solve(b,x);
for (std::size_t i = 0; i < lfsu_i.size(); ++i)
r_i.accumulate(lfsu_i,i,x[i]);
}
template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
void jacobian_skeleton(const IG& ig,
const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o,
M& mat_ii, M& mat_io,
M& mat_oi, M& mat_oo) const
{
const unsigned int size = LFSU::CHILDREN;
using namespace Indices;
using LFSU_BASE = TypeTree::Child<LFSU,_0>;
const auto& lfsu_base_i = child(lfsu_i,_0);
const auto& lfsu_base_o = child(lfsu_o,_0);
using RF = typename LFSU_BASE::Traits::FiniteElementType::
Traits::LocalBasisType::Traits::RangeFieldType;
// const int order = lfsu.finiteElement().localBasis().order();
const int intorder = 5; //FIXME
auto face_geo = ig.geometry();
for (const auto& it : quadratureRule(face_geo,intorder)) {
const auto p_local_i = ig.geometryInInside().global(it.position());
const auto p_local_o = ig.geometryInOutside().global(it.position());
auto phi_i = cache.evaluateFunction(p_local_i,lfsu_base_i.finiteElement().localBasis());
auto phi_o = cache.evaluateFunction(p_local_o,lfsu_base_o.finiteElement().localBasis());
const auto normal = ig.unitOuterNormal(it.position());
typename GF::Traits::RangeType gf_i, gf_o;
_grid_function->evaluate(ig.inside(),p_local_i,gf_i);
_grid_function->evaluate(ig.outside(),p_local_o,gf_o);
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
for (unsigned int i=0; i<size; i++)
for (unsigned int j=0; j<size; j++)
for (unsigned int k = 0; k<lfsv_i.child(i).size(); k++)
for (unsigned int l = 0; l<lfsv_o.child(j).size(); l++)
{
mat_ii.accumulate(lfsu_i.child(i),k,lfsv_i.child(j), l, (gf_i - gf_o) * normal[i] * phi_i[k] * 0.5 * factor);
mat_oo.accumulate(lfsu_o.child(i),k,lfsv_o.child(j), l, (gf_i - gf_o) * normal[j] * phi_o[l] * 0.5 * factor);
}
}
}
// define sparsity pattern of operator representation
template<typename LFSU, typename LFSV, typename LocalPattern>
void pattern_volume (const LFSU& lfsu, const LFSV& lfsv,
LocalPattern& pattern) const
{
const unsigned int size = LFSU::CHILDREN;
for (unsigned int i=0; i<size; i++)
for (unsigned int j=0; j<size; j++)
for (size_t k=0; k<lfsv.child(i).size(); ++k)
for (size_t l=0; l<lfsu.child(j).size(); ++l)
pattern.addLink(lfsv.child(i),k,lfsu.child(j),l);
}
// define sparsity pattern connecting self and neighbor dofs
template<typename LFSU, typename LFSV, typename LocalPattern>
void pattern_skeleton (const LFSU& lfsu_s, const LFSV& lfsv_s, const LFSU& lfsu_n, const LFSV& lfsv_n,
LocalPattern& pattern_sn,
LocalPattern& pattern_ns) const
{
const unsigned int size = LFSU::CHILDREN;
for (unsigned int i=0; i<size; i++)
for (unsigned int j=0; j<size; j++)
for (size_t k=0; k<lfsv_s.child(i).size(); ++k)
for (size_t l=0; l<lfsu_n.child(j).size(); ++l)
pattern_sn.addLink(lfsv_s.child(i),k,lfsu_n.child(j),l);
for (unsigned int i=0; i<size; i++)
for (unsigned int j=0; j<size; j++)
for (size_t k=0; k<lfsv_n.child(i).size(); ++k)
for (size_t l=0; l<lfsu_s.child(j).size(); ++l)
pattern_ns.addLink(lfsv_n.child(i),k,lfsu_s.child(j),l);
}
private:
std::shared_ptr<GF> _grid_function;
typedef typename FEM::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
Dune::PDELab::LocalBasisCache<LocalBasisType> cache;
};
template<class GFS,class X>
class PowerDiscreteGridFunction : public Dune::PDELab::GridFunctionBase<
Dune::PDELab::GridFunctionTraits<
typename Dune::PDELab::GridFunctionSubSpace<GFS,Dune::TypeTree::TreePath<0> >::Traits::GridViewType,
typename BasisInterfaceSwitch<
typename FiniteElementInterfaceSwitch<
typename Dune::PDELab::GridFunctionSubSpace<GFS,Dune::TypeTree::TreePath<0> >::Traits::FiniteElementType
>::Basis
>::RangeField,
GFS::CHILDREN,
Dune::FieldVector<
typename BasisInterfaceSwitch<
typename FiniteElementInterfaceSwitch<
typename Dune::PDELab::GridFunctionSubSpace<GFS,Dune::TypeTree::TreePath<0> >::Traits::FiniteElementType
>::Basis
>::RangeField,
GFS::CHILDREN
>
>,
PowerDiscreteGridFunction<GFS,X>
>
{
template<int I>
using SUBGFS = Dune::PDELab::GridFunctionSubSpace<GFS,Dune::TypeTree::TreePath<I> >;
/**
* @brief { struct_description }
*
* @tparam I Number of Children
* @tparam GFSs { description }
*/
template<int I, class...>
struct ContainerGFSHelper {
using type = typename ContainerGFSHelper<I-1,SUBGFS<I-1>>::type;
};
template<class... GFSs>
struct ContainerGFSHelper<0,GFSs...> {
using type = std::tuple<std::shared_ptr<GFSs>...>;
};
template<int I, class...>
struct ContainerDGFHelper {
using type = typename ContainerDGFHelper<I-1,Dune::PDELab::DiscreteGridFunction<SUBGFS<I-1>,X>>::type;
};
template<class... GFSs>
struct ContainerDGFHelper<0,GFSs...> {
using type = std::tuple<std::shared_ptr<GFSs>...>;
};
static constexpr int size = GFS::CHILDREN;
using ContainerDGF = typename ContainerDGFHelper<size>::type;
using ContainerGFS = typename ContainerGFSHelper<size>::type;
public:
using Traits = Dune::PDELab::GridFunctionTraits<
typename Dune::PDELab::GridFunctionSubSpace<GFS,Dune::TypeTree::TreePath<0> >::Traits::GridViewType,
typename BasisInterfaceSwitch<
typename FiniteElementInterfaceSwitch<
typename Dune::PDELab::GridFunctionSubSpace<GFS,Dune::TypeTree::TreePath<0> >::Traits::FiniteElementType
>::Basis
>::RangeField,
GFS::CHILDREN,
Dune::FieldVector<
typename BasisInterfaceSwitch<
typename FiniteElementInterfaceSwitch<
typename Dune::PDELab::GridFunctionSubSpace<GFS,Dune::TypeTree::TreePath<0> >::Traits::FiniteElementType
>::Basis
>::RangeField,
GFS::CHILDREN
>
>;
/** \brief Construct a PowerDiscreteGridFunction
*
* \param gfs The PowerGridFunctionsSpace
* \param x_ The coefficients vector
*/
PowerDiscreteGridFunction (const GFS& gfs, const X& x)
: px(stackobject_to_shared_ptr(x))
{
Dune::Hybrid::forEach(std::make_index_sequence<std::tuple_size<ContainerGFS>::value>{},
[&](auto i) {
using GFS_ptr = typename std::tuple_element<i,ContainerGFS>::type;
using GFS_i = typename GFS_ptr::element_type;
std::get<i>(_gfs) = std::make_shared<GFS_i>(gfs);
using DGF_ptr = typename std::tuple_element<i,ContainerDGF>::type;
using DGF_i = typename DGF_ptr::element_type;
std::get<i>(_dgf) = std::make_shared<DGF_i>(std::get<i>(_gfs),px);
}
);
}
/** \brief Construct a PowerDiscreteGridFunction
*
* \param gfs shared pointer to the PowerGridFunctionsSpace
* \param x_ shared pointer to the coefficients vector
*/
PowerDiscreteGridFunction (std::shared_ptr<const GFS> gfs, std::shared_ptr<const X> x)
: px(x)
{
Dune::Hybrid::forEach(std::make_index_sequence<std::tuple_size<ContainerGFS>::value>{},
[&](auto i)
{
using GFS_ptr = typename std::tuple_element<i,ContainerGFS>::type;
using GFS_i = typename GFS_ptr::element_type;
std::get<i>(_gfs) = std::make_shared<GFS_i>(gfs);
using DGF_ptr = typename std::tuple_element<i,ContainerDGF>::type;
using DGF_i = typename DGF_ptr::element_type;
std::get<i>(_dgf) = std::make_shared<DGF_i>(std::get<i>(_gfs),px);
}
);
}
void evaluate ( const typename Traits::ElementType& e,
const typename Traits::DomainType& x,
typename Traits::RangeType& y) const
{
Dune::Hybrid::forEach(std::make_index_sequence<std::tuple_size<ContainerGFS>::value>{},
[&](auto i)
{
using DGF_ptr = typename std::tuple_element<i,ContainerDGF>::type;
using DGF_i = typename DGF_ptr::element_type;
using RangeType_i = typename DGF_i::Traits::RangeType;
RangeType_i y_i;
static_assert(RangeType_i::dimension==1,"All subspaces must have a one dimensional range!");
std::get<i>(_dgf)->evaluate(e,x,y_i);
y[i] = y_i[0];
}
);
}
const typename Traits::GridViewType& getGridView() const
{
return std::get<0>(_dgf)->getGridView();
}
ContainerGFS _gfs;
ContainerDGF _dgf;
std::shared_ptr<const X> px;
};
template<class GV, class RF, int d, int k, Dune::GeometryType::BasicType gt>
struct FluxReconstructionImplemented
{
static constexpr bool value = false;
};
template<class GV, class RF, int d, int k>
struct FluxReconstructionImplemented<GV,RF,d,k,Dune::GeometryType::BasicType::simplex>
{
static constexpr bool value = true; //(d==2 && (k==0 || k==1));
};
template<class GV, class RF, int d, int k>
struct FluxReconstructionImplemented<GV,RF,d,k,Dune::GeometryType::BasicType::cube>
{
static constexpr bool value = true; // ((d==2) && (k==0 || k==1 || k==2)) || ((d==3) && (k==0 || k==1));
};
template<class GV, class RF, int d, int k, Dune::GeometryType::BasicType gt>
......@@ -102,40 +472,45 @@ private:
static constexpr int dim = d;
using DF = typename Traits::DomainFieldType;
using FEM = Dune::PDELab::RaviartThomasLocalFiniteElementMap<GV,DF,RF,k,gt>;
using CON = Dune::PDELab::NoConstraints;
using VBE = Dune::PDELab::ISTL::VectorBackend<>;
using GFS = Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE>;