Commit a11f0799 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[linear] make the ISTL solver backend work for multitype

parent 1722fad1
......@@ -30,6 +30,8 @@
#include <dune/istl/bvector.hh>
#include <dune/common/std/type_traits.hh>
#include<dumux/common/typetraits/vector.hh>
namespace Dumux {
namespace Detail {
......@@ -49,8 +51,8 @@ struct LinearAlgebraTraits
using Vector = V;
};
template<class Assembler>
struct LinearAlgebraTraitsFromAssembler
template<class Assembler, bool isMultiTypeBlockVector>
struct LATraitsFromAssemblerImpl
{
private:
using VectorPossiblyWithState = typename Assembler::ResidualType;
......@@ -62,6 +64,17 @@ public:
using Matrix = typename Assembler::JacobianMatrix;
};
template<class Assembler>
struct LATraitsFromAssemblerImpl<Assembler, true>
{
using Vector = typename Assembler::ResidualType;
using Matrix = typename Assembler::JacobianMatrix;
};
template<class Assembler>
using LinearAlgebraTraitsFromAssembler = LATraitsFromAssemblerImpl<Assembler,
isMultiTypeBlockVector<typename Assembler::ResidualType>::value>;
} // end namespace Dumux
#endif
......@@ -38,7 +38,10 @@
#include <dune/istl/scalarproducts.hh>
#include <dumux/common/typetraits/matrix.hh>
#include <dumux/common/typetraits/vector.hh>
#include <dumux/linear/linearsolverparameters.hh>
#include <dumux/linear/linearsolveracceptsmultitypematrix.hh>
#include <dumux/linear/matrixconverter.hh>
#include <dumux/linear/parallelhelpers.hh>
namespace Dumux::Detail {
......@@ -55,8 +58,8 @@ constexpr std::size_t preconditionerBlockLevel() noexcept
return isMultiTypeBlockMatrix<M>::value ? 2 : 1;
}
template<class LinearSolverTraits>
Dune::SolverCategory::Category solverCategory(const typename LinearSolverTraits::GridView& gridView)
template<class LinearSolverTraits, class GridView>
Dune::SolverCategory::Category solverCategory(const GridView& gridView)
{
if constexpr (LinearSolverTraits::canCommunicate)
{
......@@ -85,12 +88,15 @@ template<class LinearSolverTraits, class LinearAlgebraTraits,
class IstlLinearSolver
{
using Matrix = typename LinearAlgebraTraits::Matrix;
using XVector = typename InverseOperator::domain_type;
using BVector = typename InverseOperator::range_type;
using XVector = typename LinearAlgebraTraits::Vector;
using BVector = typename LinearAlgebraTraits::Vector;
using Scalar = typename InverseOperator::real_type;
#if HAVE_MPI
using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>;
using ScalarProduct = Dune::ScalarProduct<XVector>;
using ScalarProduct = Dune::ScalarProduct<typename InverseOperator::domain_type>;
using ParallelHelper = std::conditional_t<isMultiTypeBlockVector<XVector>::value,
double,
ParallelISTLHelper<LinearSolverTraits>>;
#endif
public:
......@@ -110,8 +116,9 @@ public:
/*!
* \brief Constructor for parallel and sequential solvers
*/
IstlLinearSolver(const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
template <class GridView, class DofMapper>
IstlLinearSolver(const GridView& gridView,
const DofMapper& dofMapper,
const std::string& paramGroup = "")
{
initializeParameters_(paramGroup);
......@@ -135,10 +142,11 @@ public:
/*!
* \brief Constructor with custom scalar product and communication
*/
template <class GridView, class DofMapper>
IstlLinearSolver(std::shared_ptr<Comm> communication,
std::shared_ptr<ScalarProduct> scalarProduct,
const typename LinearSolverTraits::GridView& gridView,
const typename LinearSolverTraits::DofMapper& dofMapper,
const GridView& gridView,
const DofMapper& dofMapper,
const std::string& paramGroup = "")
{
initializeParameters_(paramGroup);
......@@ -165,20 +173,27 @@ public:
Scalar norm(const XVector& x) const
{
#if HAVE_MPI
if (solverCategory_ == Dune::SolverCategory::nonoverlapping)
if constexpr (LinearSolverTraits::canCommunicate)
{
auto y(x); // make a copy because the vector needs to be made consistent
using GV = typename LinearSolverTraits::GridView;
using DM = typename LinearSolverTraits::DofMapper;
ParallelVectorHelper<GV, DM, LinearSolverTraits::dofCodim> vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper());
vectorHelper.makeNonOverlappingConsistent(y);
return scalarProduct_->norm(y);
if (solverCategory_ == Dune::SolverCategory::nonoverlapping)
{
auto y(x); // make a copy because the vector needs to be made consistent
using GV = typename LinearSolverTraits::GridView;
using DM = typename LinearSolverTraits::DofMapper;
ParallelVectorHelper<GV, DM, LinearSolverTraits::dofCodim> vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper());
vectorHelper.makeNonOverlappingConsistent(y);
return scalarProduct_->norm(y);
}
}
else
#endif
if constexpr (!preconditionerAcceptsMultiTypeMatrix<Preconditioner>::value
&& isMultiTypeBlockVector<XVector>::value)
{
return scalarProduct_->norm(x);
auto y = VectorConverter<XVector>::multiTypeToBlockVector(x);
return scalarProduct_->norm(y);
}
else
return scalarProduct_->norm(x);
}
const Dune::InverseOperatorResult& result() const
......@@ -203,13 +218,46 @@ private:
bool solveSequential_(Matrix& A, XVector& x, BVector& b)
{
// construct solver from linear operator
using SequentialTraits = typename LinearSolverTraits::template Sequential<Matrix, XVector>;
auto linearOperator = std::make_shared<typename SequentialTraits::LinearOperator>(A);
auto solver = constructPreconditionedSolver_(linearOperator);
if constexpr (preconditionerAcceptsMultiTypeMatrix<Preconditioner>::value
|| !isMultiTypeBlockMatrix<Matrix>::value)
{
// construct solver from linear operator
using SequentialTraits = typename LinearSolverTraits::template Sequential<Matrix, XVector>;
auto linearOperator = std::make_shared<typename SequentialTraits::LinearOperator>(A);
auto solver = constructPreconditionedSolver_(linearOperator);
// solve linear system
solver.apply(x, b, result_);
}
else
{
// create the bcrs matrix the IterativeSolver backend can handle
auto M = MatrixConverter<Matrix>::multiTypeToBCRSMatrix(A);
// get the new matrix sizes
const std::size_t numRows = M.N();
assert(numRows == M.M());
// create the vector the IterativeSolver backend can handle
auto bTmp = VectorConverter<BVector>::multiTypeToBlockVector(b);
assert(bTmp.size() == numRows);
// create a blockvector to which the linear solver writes the solution
using VectorBlock = typename Dune::FieldVector<Scalar, 1>;
using BlockVector = typename Dune::BlockVector<VectorBlock>;
BlockVector y(numRows);
auto linearOperator = std::make_shared<Dune::MatrixAdapter<decltype(M), decltype(y), decltype(bTmp)>>(M);
auto solver = constructPreconditionedSolver_(linearOperator);
// solve linear system
solver.apply(y, bTmp, result_);
// copy back the result y into x
if(result_.converged)
VectorConverter<XVector>::retrieveValues(x, y);
}
// solve linear system
solver.apply(x, b, result_);
return result_.converged;
}
......@@ -221,7 +269,6 @@ private:
// TODO: This can be adapted once the situation in Dune ISTL changes.
if constexpr (isMultiTypeBlockMatrix<Matrix>::value)
return solveSequential_(A, x, b);
else
{
switch (solverCategory_)
......@@ -274,7 +321,7 @@ private:
}
#if HAVE_MPI
std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> parallelHelper_;
std::unique_ptr<ParallelHelper> parallelHelper_;
std::shared_ptr<Comm> communication_;
#endif
Dune::SolverCategory::Category solverCategory_;
......@@ -285,14 +332,34 @@ private:
std::string name_;
};
template<class LSTraits, class LATraits, bool isMultiType>
struct ILUBiCGSTABIstlSolverImpl;
template<class LSTraits, class LATraits>
struct ILUBiCGSTABIstlSolverImpl<LSTraits, LATraits, false>
{
using type = IstlLinearSolver<LSTraits, LATraits,
Dune::BiCGSTABSolver<typename LATraits::Vector>,
Dune::SeqILU<typename LATraits::Matrix,
typename LATraits::Vector,
typename LATraits::Vector>
>;
};
template<class LSTraits, class LATraits>
struct ILUBiCGSTABIstlSolverImpl<LSTraits, LATraits, true>
{
using type = IstlLinearSolver<LSTraits, LATraits,
Dune::BiCGSTABSolver<decltype(VectorConverter<typename LATraits::Vector>::multiTypeToBlockVector(std::declval<typename LATraits::Vector>()))>,
Dune::SeqILU<decltype(MatrixConverter<typename LATraits::Matrix>::multiTypeToBCRSMatrix(std::declval<typename LATraits::Matrix>())),
decltype(VectorConverter<typename LATraits::Vector>::multiTypeToBlockVector(std::declval<typename LATraits::Vector>())),
decltype(VectorConverter<typename LATraits::Vector>::multiTypeToBlockVector(std::declval<typename LATraits::Vector>()))>
>;
};
template<class LSTraits, class LATraits>
using ILUBiCGSTABIstlSolver = IstlLinearSolver<LSTraits, LATraits,
Dune::BiCGSTABSolver<typename LATraits::Vector>,
Dune::SeqILU<typename LATraits::Matrix,
typename LATraits::Vector,
typename LATraits::Vector,
Detail::preconditionerBlockLevel<typename LATraits::Matrix>()>
>;
using ILUBiCGSTABIstlSolver = typename ILUBiCGSTABIstlSolverImpl<LSTraits, LATraits,
isMultiTypeBlockVector<typename LATraits::Vector>::value>::type;
} // end namespace Dumux
......
......@@ -24,7 +24,10 @@
#ifndef DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH
#define DUMUX_LINEAR_SOLVER_ACCEPTS_MULTITYPEMATRIX_HH
#include <dune/istl/preconditioners.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/istlsolvers.hh>
namespace Dumux {
......@@ -64,6 +67,12 @@ template<>
struct linearSolverAcceptsMultiTypeMatrix<UMFPackBackend> : public std::false_type {};
#endif // HAVE_UMFPACK
template<class Prec>
struct preconditionerAcceptsMultiTypeMatrix : public std::true_type {};
template<class M, class X, class Y>
struct preconditionerAcceptsMultiTypeMatrix<Dune::SeqILU<M, X, Y>> : public std::false_type {};
} // end namespace Dumux
#endif
......@@ -34,6 +34,14 @@
namespace Dumux {
namespace Detail {
template <typename T>
using GVDetector = typename T::GridView;
template <typename T>
using hasGridView = Dune::Std::is_detected<GVDetector, T>;
}
template<class LinearSolverTraits>
class LinearSolverParameters
{
......@@ -62,7 +70,10 @@ public:
params["preconditioner.relaxation"] = "1.0";
params["preconditioner.verbosity"] = "0";
params["preconditioner.defaultAggregationSizeMode"] = "isotropic";
params["preconditioner.defaultAggregationDimension"] = std::to_string(LinearSolverTraits::GridView::dimension);
if constexpr (Detail::hasGridView<LinearSolverTraits>::value)
params["preconditioner.defaultAggregationDimension"] = std::to_string(LinearSolverTraits::GridView::dimension);
else
params["preconditioner.defaultAggregationDimension"] = 3;
params["preconditioner.maxLevel"] = "100";
params["ParameterGroup"] = paramGroup;
params["preconditioner.ParameterGroup"] = paramGroup;
......
......@@ -57,6 +57,14 @@ struct SequentialSolverTraits
using Preconditioner = SeqPreconditioner;
};
struct SeqLinearSolverTraits
{
template<class Matrix, class Vector>
using Sequential = SequentialSolverTraits<Matrix, Vector>;
static constexpr bool canCommunicate = false;
};
#if HAVE_MPI
template <class MType, class VType>
struct NonoverlappingSolverTraits
......
......@@ -1223,14 +1223,14 @@ private:
assert(this->checkSizesOfSubMatrices(A) && "Sub-blocks of MultiTypeBlockMatrix have wrong sizes!");
// create the bcrs matrix the IterativeSolver backend can handle
const auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
auto M = MatrixConverter<JacobianMatrix>::multiTypeToBCRSMatrix(A);
// get the new matrix sizes
const std::size_t numRows = M.N();
assert(numRows == M.M());
// create the vector the IterativeSolver backend can handle
const auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
auto bTmp = VectorConverter<SolutionVector>::multiTypeToBlockVector(b);
assert(bTmp.size() == numRows);
// create a blockvector to which the linear solver writes the solution
......
......@@ -31,7 +31,10 @@
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/linear/istlsolvers.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/algebratraits.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/traits.hh>
......@@ -141,7 +144,8 @@ int main(int argc, char** argv)
couplingManager);
// the linear solver
using LinearSolver = ILU0BiCGSTABBackend;
using LinearSolver = ILUBiCGSTABIstlSolver<SeqLinearSolverTraits,
LinearAlgebraTraitsFromAssembler<Assembler>>;
auto linearSolver = std::make_shared<LinearSolver>();
// the non-linear solver
......
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