Commit 09174215 authored by Timo Koch's avatar Timo Koch
Browse files

[parallel] Make box work for overlapping grids

parent 11f0721b
......@@ -35,7 +35,6 @@
#include <dumux/common/indextraits.hh>
#include <dumux/common/defaultmappertraits.hh>
#include <dumux/discretization/basegridgeometry.hh>
#include <dumux/discretization/checkoverlapsize.hh>
#include <dumux/discretization/box/boxgeometryhelper.hh>
#include <dumux/discretization/box/fvelementgeometry.hh>
#include <dumux/discretization/box/subcontrolvolume.hh>
......@@ -118,12 +117,7 @@ public:
//! Constructor
BoxFVGridGeometry(const GridView gridView)
: ParentType(gridView)
{
// Check if the overlap size is what we expect
if (!CheckOverlapSize<DiscretizationMethod::box>::isValid(gridView))
DUNE_THROW(Dune::InvalidStateException, "The box discretization method only works with zero overlap for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
}
{}
//! the vertex mapper is the dofMapper
//! this is convenience to have better chance to have the same main files for box/tpfa/mpfa...
......@@ -391,12 +385,7 @@ public:
//! Constructor
BoxFVGridGeometry(const GridView gridView)
: ParentType(gridView)
{
// Check if the overlap size is what we expect
if (!CheckOverlapSize<DiscretizationMethod::box>::isValid(gridView))
DUNE_THROW(Dune::InvalidStateException, "The box discretization method only works with zero overlap for parallel computations. "
<< " Set the parameter \"Grid.Overlap\" in the input file.");
}
{}
//! the vertex mapper is the dofMapper
//! this is convenience to have better chance to have the same main files for box/tpfa/mpfa...
......
......@@ -130,7 +130,18 @@ private:
if constexpr (LinearSolverTraits::canCommunicate)
{
if (isParallel_)
solveParallel_(A, x, b);
{
if (LinearSolverTraits::isNonOverlapping(phelper_->gridView()))
{
using PTraits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>;
solveParallel_<PTraits>(A, x, b);
}
else
{
using PTraits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>;
solveParallel_<PTraits>(A, x, b);
}
}
else
solveSequential_(A, x, b);
}
......@@ -140,21 +151,20 @@ private:
}
}
template<class Matrix, class Vector>
template<class ParallelTraits, class Matrix, class Vector>
void solveParallel_(Matrix& A, Vector& x, Vector& b)
{
using Traits = typename LinearSolverTraits::template Parallel<Matrix, Vector>;
using Comm = typename Traits::Comm;
using LinearOperator = typename Traits::LinearOperator;
using ScalarProduct = typename Traits::ScalarProduct;
using Comm = typename ParallelTraits::Comm;
using LinearOperator = typename ParallelTraits::LinearOperator;
using ScalarProduct = typename ParallelTraits::ScalarProduct;
std::shared_ptr<Comm> comm;
std::shared_ptr<LinearOperator> linearOperator;
std::shared_ptr<ScalarProduct> scalarProduct;
prepareLinearAlgebraParallel<LinearSolverTraits>(A, b, comm, linearOperator, scalarProduct, *phelper_, firstCall_);
prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, comm, linearOperator, scalarProduct, *phelper_, firstCall_);
using SeqSmoother = Dune::SeqSSOR<Matrix, Vector, Vector>;
using Smoother = typename Traits::template Preconditioner<SeqSmoother>;
using Smoother = typename ParallelTraits::template Preconditioner<SeqSmoother>;
solveWithAmg_<Smoother>(A, x, b, linearOperator, comm, scalarProduct);
}
#endif // HAVE_MPI
......
......@@ -94,6 +94,7 @@ public:
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>;
using LinearOperator = Dune::NonoverlappingSchwarzOperator<MType, VType, VType, Comm>;
using ScalarProduct = Dune::NonoverlappingSchwarzScalarProduct<VType, Comm>;
static constexpr bool isNonOverlapping = true;
template<class SeqPreconditioner>
using Preconditioner = Dune::NonoverlappingBlockPreconditioner<Comm, SeqPreconditioner>;
......@@ -108,58 +109,65 @@ public:
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>;
using LinearOperator = Dune::OverlappingSchwarzOperator<MType, VType, VType, Comm>;
using ScalarProduct = Dune::OverlappingSchwarzScalarProduct<VType, Comm>;
static constexpr bool isNonOverlapping = false;
template<class SeqPreconditioner>
using Preconditioner = Dune::BlockPreconditioner<VType, VType, Comm, SeqPreconditioner>;
};
#endif
//! Box: use non-overlapping model
template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethod::box>
struct LinearSolverTraitsBase
{
using GridView = typename GridGeometry::GridView;
using Grid = typename GridGeometry::GridView::Traits::Grid;
using DofMapper = typename GridGeometry::VertexMapper;
template<class Matrix, class Vector>
using Sequential = SequentialSolverTraits<Matrix, Vector>;
template<class Matrix, class Vector>
using ParallelOverlapping = OverlappingSolverTraits<Matrix, Vector>;
template<class Matrix, class Vector>
using ParallelNonoverlapping = NonoverlappingSolverTraits<Matrix, Vector>;
};
//! Box: use overlapping or non-overlapping model depending on the grid
template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethod::box>
: public LinearSolverTraitsBase<GridGeometry>
{
using DofMapper = typename GridGeometry::VertexMapper;
using Grid = typename GridGeometry::GridView::Traits::Grid;
static constexpr int dofCodim = Grid::dimension;
static constexpr bool isNonOverlapping = true;
// TODO: see above for description of this workaround, remove second line if fixed upstream
static constexpr bool canCommunicate =
Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
|| Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v;
template<class Matrix, class Vector>
using Sequential = SequentialSolverTraits<Matrix, Vector>;
template<class Matrix, class Vector>
using Parallel = std::conditional_t<isNonOverlapping,
NonoverlappingSolverTraits<Matrix, Vector>,
OverlappingSolverTraits<Matrix, Vector>>;
template<class GridView>
static bool isNonOverlapping(const GridView& gridView)
{ return gridView.overlapSize(0) == 0; }
};
//! Cell-centered tpfa: use overlapping model
template<class GridGeometry>
struct LinearSolverTraitsImpl<GridGeometry, DiscretizationMethod::cctpfa>
: public LinearSolverTraitsBase<GridGeometry>
{
using GridView = typename GridGeometry::GridView;
using DofMapper = typename GridGeometry::ElementMapper;
using Grid = typename GridGeometry::GridView::Traits::Grid;
using DofMapper = typename GridGeometry::VertexMapper;
static constexpr int dofCodim = 0;
static constexpr bool isNonOverlapping = false;
// TODO: see above for description of this workaround, remove second line if fixed upstream
static constexpr bool canCommunicate =
Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
|| Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v;
template<class Matrix, class Vector>
using Sequential = SequentialSolverTraits<Matrix, Vector>;
template<class Matrix, class Vector>
using Parallel = std::conditional_t<isNonOverlapping,
NonoverlappingSolverTraits<Matrix, Vector>,
OverlappingSolverTraits<Matrix, Vector>>;
template<class GridView>
static bool isNonOverlapping(const GridView& gridView)
{ return false; }
};
//! Cell-centered mpfa: use overlapping model
......
......@@ -786,21 +786,22 @@ private:
/*!
* \brief Prepare linear algebra variables for parallel solvers
*/
template<class LinearSolverTraits, class Matrix, class Vector, class ParallelHelper>
template<class LinearSolverTraits, class ParallelTraits,
class Matrix, class Vector, class ParallelHelper>
void prepareLinearAlgebraParallel(Matrix& A, Vector& b,
std::shared_ptr<typename LinearSolverTraits::template Parallel<Matrix, Vector>::Comm>& comm,
std::shared_ptr<typename LinearSolverTraits::template Parallel<Matrix, Vector>::LinearOperator>& fop,
std::shared_ptr<typename LinearSolverTraits::template Parallel<Matrix, Vector>::ScalarProduct>& sp,
ParallelHelper& pHelper,
const bool firstCall)
std::shared_ptr<typename ParallelTraits::Comm>& comm,
std::shared_ptr<typename ParallelTraits::LinearOperator>& fop,
std::shared_ptr<typename ParallelTraits::ScalarProduct>& sp,
ParallelHelper& pHelper,
const bool firstCall)
{
if (firstCall) pHelper.initGhostsAndOwners();
if (LinearSolverTraits::isNonOverlapping)
if constexpr (ParallelTraits::isNonOverlapping)
{
// extend the matrix pattern such that it is usable for a parallel solver
// and make right-hand side consistent
using GridView = std::decay_t<decltype(pHelper.gridView())>;
using GridView = typename LinearSolverTraits::GridView;
using DofMapper = typename LinearSolverTraits::DofMapper;
static constexpr int dofCodim = LinearSolverTraits::dofCodim;
ParallelMatrixHelper<Matrix, GridView, DofMapper, dofCodim> matrixHelper(pHelper.gridView(), pHelper.dofMapper());
......@@ -810,7 +811,7 @@ void prepareLinearAlgebraParallel(Matrix& A, Vector& b,
pHelper.makeNonOverlappingConsistent(b);
// create commicator, operator, scalar product
using Traits = typename LinearSolverTraits::template Parallel<Matrix, Vector>;
using Traits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>;
const auto category = Dune::SolverCategory::nonoverlapping;
comm = std::make_shared<typename Traits::Comm>(pHelper.gridView().comm(), category);
pHelper.createParallelIndexSet(*comm);
......@@ -820,7 +821,7 @@ void prepareLinearAlgebraParallel(Matrix& A, Vector& b,
else
{
// create commicator, operator, scalar product
using Traits = typename LinearSolverTraits::template Parallel<Matrix, Vector>;
using Traits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>;
const auto category = Dune::SolverCategory::overlapping;
comm = std::make_shared<typename Traits::Comm>(pHelper.gridView().comm(), category);
pHelper.createParallelIndexSet(*comm);
......
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