...
  View open merge request
Commits (17)
......@@ -82,7 +82,7 @@ private:
using FEMU = Dune::Dorie::
RaviartThomasLocalFiniteElementMap<GV,DF,RF,order,gt>;
using ES = Dune::PDELab::OverlappingEntitySet<GV>;
using ES = Dune::PDELab::AllEntitySet<GV>;
using GFSU = Dune::PDELab::GridFunctionSpace<ES,FEMU>;
using X = Dune::PDELab::Backend::Vector<GFSU, RF>;
using GFDGPiola = Dune::PDELab::DiscreteGridFunctionPiola<GFSU,X>;
......
......@@ -71,9 +71,9 @@ public:
_config(config),
_log(Dorie::get_logger(log_base))
{
check_parallel_allowed();
const auto grid_type = _config.get<std::string>("grid.gridType");
check_parallel_allowed(grid_type);
if (grid_type == "rectangular")
{
_grid_mode = GridMode::rectangular;
......@@ -128,13 +128,15 @@ public:
/** Currently, we only allow parallel execution on YaspGrid due to
* constraints imposed by the linear solver.
*/
bool check_parallel_allowed () const
bool check_parallel_allowed (const std::string& grid_type) const
{
if (_helper.size() > 1
&& not std::is_same_v<YaspGrid<Grid::dimension>, Grid>) {
&& not std::is_same_v<YaspGrid<Grid::dimension>, Grid>
&& grid_type == "rectangular") {
_log->warn("DORiE's solver does not support parallel "
"execution on unstructured grids. This will likely "
"result in a fatal exception (later on)");
"execution on non-conforming grids. This will likely "
"result in a fatal exception or a MPI deadlock "
"(later on)");
return false;
}
......
......@@ -10,8 +10,11 @@
#include <dune/geometry/type.hh>
#include <dune/grid/uggrid.hh>
#include <dune/pdelab/gridfunctionspace/gridfunctionspace.hh>
#include <dune/pdelab/constraints/p0.hh>
#include <dune/pdelab/constraints/p0ghost.hh>
#include <dune/pdelab/finiteelementmap/p0fem.hh>
#include <dune/pdelab/finiteelementmap/qkdg.hh>
#include <dune/dorie/common/finite_element/pk_dg_fem.hh>
......@@ -53,6 +56,16 @@ constexpr R estimate_mbe_entries (const int dim, const Dune::GeometryType::Basic
"for given geometry type!");
}
/// Provides correct constraints for parallel overlapping methods in DG or FV
template<typename Grid>
struct P0ParallelConstraints
: public Dune::PDELab::P0ParallelConstraints {};
/// Provides correct constraints for parallel overlapping methods in DG or FV
template<int dim>
struct P0ParallelConstraints<Dune::UGGrid<dim>>
: public Dune::PDELab::P0ParallelGhostConstraints {};
/// Provide types and construction of the GFS for the linear solver
template<typename GridView, typename RF, Dune::GeometryType::BasicType GeometryType>
struct LinearSolverGridFunctionSpaceHelper
......@@ -68,11 +81,11 @@ private:
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
using ES = Dune::PDELab::AllEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
using CON = Dune::Dorie::P0ParallelConstraints<typename GridView::Grid>;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::ISTL::VectorBackend<>>;
......@@ -99,11 +112,11 @@ private:
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
using ES = Dune::PDELab::AllEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
using CON = Dune::Dorie::P0ParallelConstraints<typename GridView::Grid>;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::ISTL::VectorBackend<>>;
......@@ -137,11 +150,11 @@ private:
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
using ES = Dune::PDELab::AllEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::Dorie::PkDGLocalFiniteElementMap<DF, RF, dim, order>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
using CON = Dune::Dorie::P0ParallelConstraints<typename GridView::Grid>;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::ISTL::VectorBackend<
......@@ -173,11 +186,11 @@ private:
public:
/// Entity set of the GFS
using ES = Dune::PDELab::OverlappingEntitySet<GridView>;
using ES = Dune::PDELab::AllEntitySet<GridView>;
/// FiniteElementMap type of GFS
using FEM = typename Dune::PDELab::QkDGLocalFiniteElementMap<DF,RF,order,dim>;
/// Constraints type of the GFS
using CON = Dune::PDELab::P0ParallelConstraints;
using CON = Dune::Dorie::P0ParallelConstraints<typename GridView::Grid>;
/// GFS type
using Type = typename Dune::PDELab::GridFunctionSpace<ES,FEM,CON,
Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::fixed>
......
......@@ -121,6 +121,7 @@ void RichardsSimulation<Traits>::operator_setup()
this->_log->debug(" Total number of DOF: {}", dof_count);
Dune::PDELab::constraints(*this->gfs,*this->cc,false);
// --- Grid Operators ---
go0 = std::make_unique<GO0>(*gfs,*cc,*gfs,*cc,*slop,mbe_slop);
go1 = std::make_unique<GO1>(*gfs,*cc,*gfs,*cc,*tlop,mbe_tlop);
......@@ -145,6 +146,13 @@ void RichardsSimulation<Traits>::operator_setup()
osm->setVerbosityLevel(0);
gfs->update();
// Comunicate DOF in ghost cells
if (gfs->gridView().comm().size()>1) {
Dune::PDELab::CopyDataHandle<GFS,U> handler(*gfs,*u);
auto interface = Dune::InteriorBorder_All_Interface;
gfs->gridView().communicate(handler,interface,Dune::ForwardCommunication);
}
}
template<typename Traits>
......
......@@ -128,6 +128,13 @@ void TransportSimulation<Traits>::operator_setup()
osm->setVerbosityLevel(0);
gfs->update();
// Comunicate DOF in ghost cells
if (gfs->gridView().comm().size()>1) {
Dune::PDELab::CopyDataHandle<GFS,U> handler(*gfs,*u);
auto interface = Dune::InteriorBorder_All_Interface;
gfs->gridView().communicate(handler,interface,Dune::ForwardCommunication);
}
}
template<typename Traits>
......