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

Merge branch 'feature/cleanup-parallel-solvers' into 'master'

Feature/cleanup parallel solvers

See merge request !1848
parents 47163f51 0fb06009
......@@ -31,7 +31,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/discretization/method.hh>
#include <dumux/parallel/vertexhandles.hh>
#include <dumux/parallel/vectorcommdatahandle.hh>
#include "jacobianpattern.hh"
#include "diffmethod.hh"
......@@ -191,8 +191,8 @@ public:
if (isBox && gridView().comm().size() > 1)
{
using VertexMapper = typename GridGeometry::VertexMapper;
VertexHandleSum<typename SolutionVector::block_type, SolutionVector, VertexMapper>
sumResidualHandle(residual, gridGeometry_->vertexMapper());
VectorCommDataHandleSum<VertexMapper, SolutionVector, GridGeometry::GridView::dimension>
sumResidualHandle(gridGeometry_->vertexMapper(), residual);
gridView().communicate(sumResidualHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
......
......@@ -32,7 +32,7 @@
#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/discretization/method.hh>
#include <dumux/parallel/vertexhandles.hh>
#include <dumux/parallel/vectorcommdatahandle.hh>
#include "entitycolor.hh"
......@@ -181,8 +181,8 @@ public:
// at this point we communicate the yellow vertices to the
// neighboring processes because a neigbor process may not see
// the red vertex for yellow border vertices
VertexHandleMin<EntityColor, std::vector<EntityColor>, VertexMapper>
minHandle(vertexColor_, vertexMapper);
VectorCommDataHandleMin<VertexMapper, std::vector<EntityColor>, dim>
minHandle(vertexMapper, vertexColor_);
gridView.communicate(minHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
......@@ -235,8 +235,8 @@ public:
}
// demote the border orange vertices
VertexHandleMax<EntityColor, std::vector<EntityColor>, VertexMapper>
maxHandle(vertexColor_, vertexMapper);
VectorCommDataHandleMax<VertexMapper, std::vector<EntityColor>, dim>
maxHandle(vertexMapper, vertexColor_);
gridView.communicate(maxHandle,
Dune::InteriorBorder_InteriorBorder_Interface,
Dune::ForwardCommunication);
......
......@@ -2,8 +2,10 @@ install(FILES
amgbackend.hh
amgparallelhelpers.hh
amgtraits.hh
linearsolvertraits.hh
linearsolveracceptsmultitypematrix.hh
matrixconverter.hh
parallelhelpers.hh
pdesolver.hh
scotchbackend.hh
seqsolverbackend.hh
......
......@@ -43,7 +43,7 @@
#include <dune/istl/solvers.hh>
#include <dumux/linear/solver.hh>
#include <dumux/linear/amgparallelhelpers.hh>
#include <dumux/linear/parallelhelpers.hh>
namespace Dumux {
......@@ -134,12 +134,18 @@ public:
template<class Matrix, class Vector>
bool solve(Matrix& A, Vector& x, Vector& b)
{
int rank = 0;
std::shared_ptr<Comm> comm;
std::shared_ptr<LinearOperator> fop;
std::shared_ptr<ScalarProduct> sp;
static const bool isParallel = AmgTraits::isParallel;
prepareLinearAlgebra_<Matrix, Vector, isParallel>(A, b, rank, comm, fop, sp);
#if HAVE_MPI
if constexpr (AmgTraits::isParallel)
prepareLinearAlgebraParallel<AmgTraits>(A, b, comm, fop, sp, *phelper_, firstCall_);
else
prepareLinearAlgebraSequential<AmgTraits>(A, comm, fop, sp);
#else
prepareLinearAlgebraSequential<AmgTraits>(A, comm, fop, sp);
#endif
using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat, Dune::Amg::FirstDiagonal>>;
......@@ -156,7 +162,7 @@ public:
AMGType amg(*fop, criterion, smootherArgs, *comm);
Dune::BiCGSTABSolver<VType> solver(*fop, *sp, amg, this->residReduction(), this->maxIter(),
rank == 0 ? this->verbosity() : 0);
comm->communicator().rank() == 0 ? this->verbosity() : 0);
solver.apply(x, b, result_);
firstCall_ = false;
......@@ -181,34 +187,6 @@ public:
private:
/*!
* \brief Prepare the linear algebra member variables.
*
* At compile time, correct constructor calls have to be chosen,
* depending on whether the setting is parallel or sequential.
* Since several template parameters are present, this cannot be solved
* by a full function template specialization. Instead, class template
* specialization has to be used.
* This adapts example 4 from http://www.gotw.ca/publications/mill17.htm.
* The function is called from the solve function. The call is
* forwarded to the corresponding function of a class template.
*
* \tparam Matrix the matrix type
* \tparam Vector the vector type
* \tparam isParallel decides if the setting is parallel or sequential
*/
template<class Matrix, class Vector, bool isParallel>
void prepareLinearAlgebra_(Matrix& A, Vector& b, int& rank,
std::shared_ptr<Comm>& comm,
std::shared_ptr<LinearOperator>& fop,
std::shared_ptr<ScalarProduct>& sp)
{
LinearAlgebraPreparator<GridView, AmgTraits, isParallel>
::prepareLinearAlgebra(A, b, rank, comm, fop, sp,
*phelper_, firstCall_);
}
std::shared_ptr<ParallelISTLHelper<GridView, AmgTraits>> phelper_;
Dune::InverseOperatorResult result_;
bool firstCall_;
......@@ -217,7 +195,7 @@ private:
} // namespace Dumux
#include <dumux/common/properties.hh>
#include <dumux/linear/amgtraits.hh>
#include <dumux/linear/linearsolvertraits.hh>
namespace Dumux {
......@@ -228,9 +206,9 @@ namespace Dumux {
* \note This is an adaptor using a TypeTag
*/
template<class TypeTag>
using AMGBackend = ParallelAMGBackend<GetPropType<TypeTag, Properties::GridView>, AmgTraits<GetPropType<TypeTag, Properties::JacobianMatrix>,
GetPropType<TypeTag, Properties::SolutionVector>,
GetPropType<TypeTag, Properties::GridGeometry>>>;
using AMGBackend = ParallelAMGBackend<GetPropType<TypeTag, Properties::GridView>, LinearSolverTraits<GetPropType<TypeTag, Properties::JacobianMatrix>,
GetPropType<TypeTag, Properties::SolutionVector>,
GetPropType<TypeTag, Properties::GridGeometry>>>;
} // namespace Dumux
......
This diff is collapsed.
......@@ -24,151 +24,16 @@
#ifndef DUMUX_AMG_TRAITS_HH
#define DUMUX_AMG_TRAITS_HH
#include <dune/istl/schwarz.hh>
#include <dune/istl/novlpschwarz.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/grid/common/capabilities.hh>
#warning "This header is deprecated and will be removed after release 3.2. Use linearsolver/linearsolvertraits.hh"
#include <dumux/discretization/method.hh>
#include "linearsolvertraits.hh"
// TODO: The following is a temporary solution to make the parallel AMG work
// for UGGrid. Once it is resolved upstream
// (https://gitlab.dune-project.org/core/dune-grid/issues/78),
// it should be guarded by a DUNE_VERSION macro and removed later.
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif // HAVE_UG
namespace Dumux {
namespace Temp {
namespace Capabilities {
template<class Grid, int codim>
struct canCommunicate
{
static const bool v = false;
};
#if HAVE_UG
template<int dim, int codim>
struct canCommunicate<Dune::UGGrid<dim>, codim>
{
static const bool v = true;
};
#endif // HAVE_UG
} // namespace Capabilities
} // namespace Temp
} // namespace Dumux
// end workaround
namespace Dumux {
//! The implementation is specialized for the different discretizations
template<class MType, class VType, class GridGeometry, DiscretizationMethod discMethod> struct AmgTraitsImpl;
//! The type traits required for using the AMG backend
template<class MType, class VType, class GridGeometry>
using AmgTraits = AmgTraitsImpl<MType, VType, GridGeometry, GridGeometry::discMethod>;
//! NonoverlappingSolverTraits used by discretization with non-overlapping parallel model
template <class MType, class VType, bool isParallel>
class NonoverlappingSolverTraits
{
public:
using Comm = Dune::Amg::SequentialInformation;
using LinearOperator = Dune::MatrixAdapter<MType,VType,VType>;
using ScalarProduct = Dune::SeqScalarProduct<VType>;
using Smoother = Dune::SeqSSOR<MType,VType, VType>;
};
#if HAVE_MPI
template <class MType, class VType>
class NonoverlappingSolverTraits<MType, VType, true>
{
public:
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>,int>;
using LinearOperator = Dune::NonoverlappingSchwarzOperator<MType,VType, VType,Comm>;
using ScalarProduct = Dune::NonoverlappingSchwarzScalarProduct<VType,Comm>;
using Smoother = Dune::NonoverlappingBlockPreconditioner<Comm,Dune::SeqSSOR<MType,VType, VType> >;
};
#endif
//! Box: use the non-overlapping AMG
template<class Matrix, class Vector, class GridGeometry>
struct AmgTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::box>
{
using Grid = typename GridGeometry::GridView::Traits::Grid;
enum {
dofCodim = Grid::dimension,
isNonOverlapping = true,
// TODO: see above for description of this workaround, remove second line if fixed upstream
isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
|| Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v
};
using MType = Matrix;
using VType = Dune::BlockVector<Dune::FieldVector<typename Vector::block_type::value_type, Vector::block_type::dimension>>;
using SolverTraits = NonoverlappingSolverTraits<MType, VType, isParallel>;
using Comm = typename SolverTraits::Comm;
using LinearOperator = typename SolverTraits::LinearOperator;
using ScalarProduct = typename SolverTraits::ScalarProduct;
using Smoother = typename SolverTraits::Smoother;
using DofMapper = typename GridGeometry::VertexMapper;
};
//! OverlappingSolverTraits used by discretization with overlapping parallel model
template <class MType, class VType, bool isParallel>
class OverlappingSolverTraits
{
public:
using Comm = Dune::Amg::SequentialInformation;
using LinearOperator = Dune::MatrixAdapter<MType,VType,VType>;
using ScalarProduct = Dune::SeqScalarProduct<VType>;
using Smoother = Dune::SeqSSOR<MType,VType, VType>;
};
#if HAVE_MPI
template <class MType, class VType>
class OverlappingSolverTraits<MType, VType, true>
{
public:
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>,int>;
using LinearOperator = Dune::OverlappingSchwarzOperator<MType,VType, VType,Comm>;
using ScalarProduct = Dune::OverlappingSchwarzScalarProduct<VType,Comm>;
using Smoother = Dune::BlockPreconditioner<VType,VType,Comm,Dune::SeqSSOR<MType,VType, VType> >;
};
#endif
//! Cell-centered tpfa: use the overlapping AMG
template<class Matrix, class Vector, class GridGeometry>
struct AmgTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::cctpfa>
{
using Grid = typename GridGeometry::GridView::Traits::Grid;
enum {
dofCodim = 0,
isNonOverlapping = false,
// TODO: see above for description of this workaround, remove second line if fixed upstream
isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
|| Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v
};
using MType = Matrix;
using VType = Dune::BlockVector<Dune::FieldVector<typename Vector::block_type::value_type, Vector::block_type::dimension>>;
using SolverTraits = OverlappingSolverTraits<MType, VType, isParallel>;
using Comm = typename SolverTraits::Comm;
using LinearOperator = typename SolverTraits::LinearOperator;
using ScalarProduct = typename SolverTraits::ScalarProduct;
using Smoother = typename SolverTraits::Smoother;
using DofMapper = typename GridGeometry::ElementMapper;
};
template<class Matrix, class Vector, class GridGeometry>
struct AmgTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::ccmpfa>
: public AmgTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::cctpfa> {};
using AmgTraits [[deprecated("Use LinearSolverTraits<MType, VType, GridGeometry> instead. AmgTraits will be removed after 3.2!")]] = LinearSolverTraits<MType, VType, GridGeometry>;
} // end namespace Dumux
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup Linear
* \brief Define traits for linear solvers.
*/
#ifndef DUMUX_LINEAR_SOLVER_TRAITS_HH
#define DUMUX_LINEAR_SOLVER_TRAITS_HH
#include <dune/istl/schwarz.hh>
#include <dune/istl/novlpschwarz.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/grid/common/capabilities.hh>
#include <dumux/discretization/method.hh>
// TODO: The following is a temporary solution to make the parallel AMG work
// for UGGrid. Once it is resolved upstream
// (https://gitlab.dune-project.org/core/dune-grid/issues/78),
// it should be guarded by a DUNE_VERSION macro and removed later.
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif // HAVE_UG
namespace Dumux::Temp::Capabilities {
template<class Grid, int codim>
struct canCommunicate
{
static const bool v = false;
};
#if HAVE_UG
template<int dim, int codim>
struct canCommunicate<Dune::UGGrid<dim>, codim>
{
static const bool v = true;
};
#endif // HAVE_UG
} // namespace Dumux::Temp::Capabilities
// end workaround
namespace Dumux {
//! The implementation is specialized for the different discretizations
template<class MType, class VType, class GridGeometry, DiscretizationMethod discMethod> struct LinearSolverTraitsImpl;
//! The type traits required for using the IstlFactoryBackend
template<class MType, class VType, class GridGeometry>
using LinearSolverTraits = LinearSolverTraitsImpl<MType, VType, GridGeometry, GridGeometry::discMethod>;
//! NonoverlappingSolverTraits used by discretization with non-overlapping parallel model
template <class MType, class VType, bool isParallel>
class NonoverlappingSolverTraits
{
public:
using Comm = Dune::Amg::SequentialInformation;
using LinearOperator = Dune::MatrixAdapter<MType,VType,VType>;
using ScalarProduct = Dune::SeqScalarProduct<VType>;
using Smoother = Dune::SeqSSOR<MType,VType, VType>;
};
#if HAVE_MPI
template <class MType, class VType>
class NonoverlappingSolverTraits<MType, VType, true>
{
public:
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>,int>;
using LinearOperator = Dune::NonoverlappingSchwarzOperator<MType,VType, VType,Comm>;
using ScalarProduct = Dune::NonoverlappingSchwarzScalarProduct<VType,Comm>;
using Smoother = Dune::NonoverlappingBlockPreconditioner<Comm,Dune::SeqSSOR<MType,VType, VType> >;
};
#endif
//! Box: use non-overlapping model
template<class Matrix, class Vector, class GridGeometry>
struct LinearSolverTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::box>
{
using Grid = typename GridGeometry::GridView::Traits::Grid;
enum {
dofCodim = Grid::dimension,
isNonOverlapping = true,
// TODO: see above for description of this workaround, remove second line if fixed upstream
isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
|| Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v
};
using MType = Matrix;
using VType = Dune::BlockVector<Dune::FieldVector<typename Vector::block_type::value_type, Vector::block_type::dimension>>;
using SolverTraits = NonoverlappingSolverTraits<MType, VType, isParallel>;
using Comm = typename SolverTraits::Comm;
using LinearOperator = typename SolverTraits::LinearOperator;
using ScalarProduct = typename SolverTraits::ScalarProduct;
using Smoother = typename SolverTraits::Smoother;
using DofMapper = typename GridGeometry::VertexMapper;
};
//! OverlappingSolverTraits used by discretization with overlapping parallel model
template <class MType, class VType, bool isParallel>
class OverlappingSolverTraits
{
public:
using Comm = Dune::Amg::SequentialInformation;
using LinearOperator = Dune::MatrixAdapter<MType,VType,VType>;
using ScalarProduct = Dune::SeqScalarProduct<VType>;
using Smoother = Dune::SeqSSOR<MType,VType, VType>;
};
#if HAVE_MPI
template <class MType, class VType>
class OverlappingSolverTraits<MType, VType, true>
{
public:
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>,int>;
using LinearOperator = Dune::OverlappingSchwarzOperator<MType,VType, VType,Comm>;
using ScalarProduct = Dune::OverlappingSchwarzScalarProduct<VType,Comm>;
using Smoother = Dune::BlockPreconditioner<VType,VType,Comm,Dune::SeqSSOR<MType,VType, VType> >;
};
#endif
//! Cell-centered tpfa: use overlapping model
template<class Matrix, class Vector, class GridGeometry>
struct LinearSolverTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::cctpfa>
{
using Grid = typename GridGeometry::GridView::Traits::Grid;
enum {
dofCodim = 0,
isNonOverlapping = false,
// TODO: see above for description of this workaround, remove second line if fixed upstream
isParallel = Dune::Capabilities::canCommunicate<Grid, dofCodim>::v
|| Dumux::Temp::Capabilities::canCommunicate<Grid, dofCodim>::v
};
using MType = Matrix;
using VType = Dune::BlockVector<Dune::FieldVector<typename Vector::block_type::value_type, Vector::block_type::dimension>>;
using SolverTraits = OverlappingSolverTraits<MType, VType, isParallel>;
using Comm = typename SolverTraits::Comm;
using LinearOperator = typename SolverTraits::LinearOperator;
using ScalarProduct = typename SolverTraits::ScalarProduct;
using Smoother = typename SolverTraits::Smoother;
using DofMapper = typename GridGeometry::ElementMapper;
};
template<class Matrix, class Vector, class GridGeometry>
struct LinearSolverTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::ccmpfa>
: public LinearSolverTraitsImpl<Matrix, Vector, GridGeometry, DiscretizationMethod::cctpfa> {};
} // end namespace Dumux
#endif // DUMUX_LINEAR_SOLVER_TRAITS_HH
This diff is collapsed.
......@@ -24,75 +24,13 @@
#ifndef DUMUX_VECTOR_EXCHANGE_HH
#define DUMUX_VECTOR_EXCHANGE_HH
#include <dune/grid/common/datahandleif.hh>
#warning "This header is deprecated and will be removed after release 3.2. Use linearsolver/linearsolvertraits.hh"
namespace Dumux {
/*!
* \ingroup Linear
* \todo why is this needed? is parallel/vectorexhange.hh not the same?
* \brief A data handle class to exchange entries of a vector
*/
template<class Mapper, class Vector> // mapper type and vector type
class VectorExchange
: public Dune::CommDataHandleIF<VectorExchange<Mapper,Vector>,
typename Vector::value_type>
{
public:
//! export type of data for message buffer
using DataType = typename Vector::value_type;
//! returns true if data for this codim should be communicated
bool contains (int dim, int codim) const
{
return (codim == 0);
}
//! returns true if size per entity of given dim and codim is a constant
bool fixedsize (int dim, int codim) const
{
return true;
}
/*!
* \brief how many objects of type DataType have to be sent for a given entity
* \note Only the sender side needs to know this size.
*/
template<class Entity>
size_t size (Entity& entity) const
{
return 1;
}
//! pack data from user to message buffer
template<class MessageBuffer, class Entity>
void gather (MessageBuffer& buff, const Entity& entity) const
{
buff.write(dataVector_[mapper_.index(entity)]);
}
/*!
* \brief unpack data from message buffer to user
* \note n is the number of objects sent by the sender
*/
template<class MessageBuffer, class Entity>
void scatter (MessageBuffer& buff, const Entity& entity, size_t n)
{
DataType x;
buff.read(x);
dataVector_[mapper_.index(entity)] = x;
}
//! constructor
VectorExchange (const Mapper& mapper, Vector& dataVector)
: mapper_(mapper), dataVector_(dataVector)
{}
private:
const Mapper& mapper_;
Vector& dataVector_;
};
#include <dumux/parallel/vectorcommdatahandle.hh>
namespace Dumux {
template<class Mapper, class Vector>
using VectorExchange [[deprecated("Use VectorCommDataHandleEqual<Mapper, Vector, 0> instead. Will be removed after 3.2!")]] = VectorCommDataHandleEqual<Mapper, Vector, 0/*elementCodim*/>;
} // end namespace Dumux
#endif
install(FILES
vectorcommdatahandle.hh
vertexhandles.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/parallel)
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup Linear
* \brief Contains a class to exchange entries of a vector
*/
#ifndef DUMUX_VECTOR_COMM_DATA_HANDLE_HH
#define DUMUX_VECTOR_COMM_DATA_HANDLE_HH