Commit a5743f27 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

remove comm from solver

parent 752f128d
......@@ -23,7 +23,7 @@ add_subdirectory("src")
add_subdirectory("test")
target_link_libraries(amdis fmt)
target_compile_options(amdis PUBLIC "-Wfatal-errors" "-Wall" "-Wpedantic")
target_compile_options(amdis PUBLIC "-Wall" "-Wpedantic")
option(ENABLE_ALL_WARNINGS "enable all meaningful warnings" OFF)
if (ENABLE_ALL_WARNINGS)
......
......@@ -12,6 +12,7 @@ struct NavierStokesBasis
using Grid = Dune::YaspGrid<GRIDDIM>;
using GlobalBasis = typename TaylorHoodBasis<Grid>::GlobalBasis;
using CoefficientType = double;
using LinAlgTraits = BackendTraits<GlobalBasis>;
};
int main(int argc, char** argv)
......
......@@ -21,7 +21,7 @@ namespace AMDiS
* \tparam CB Basis of matrix columns
* \tparam T Coefficient type to store in the matrix
**/
template <class RB, class CB, class T = double, class Traits = BackendTraits<RB,T>>
template <class RB, class CB, class T = double, class Traits = BackendTraits<RB>>
class BiLinearForm
: public MatrixFacade<T, Traits::template MatrixImpl>
, private ObserverSequence<event::adapt,2>
......@@ -39,7 +39,7 @@ namespace AMDiS
using ColLocalView = typename ColBasis::LocalView;
/// The type of the elements of the DOFVector
using CoefficientType = typename Traits::CoefficientType;
using CoefficientType = T;
/// The type of the matrix filled on an element with local contributions
using ElementMatrix = FlatMatrix<CoefficientType>;
......
......@@ -35,7 +35,7 @@ namespace AMDiS
* \tparam T Type of the coefficients
* \tparam TraitsType Collection of parameter for the linear-algebra backend
**/
template <class GB, class T = double, class TraitsType = BackendTraits<GB,T>>
template <class GB, class T = double, class TraitsType = BackendTraits<GB>>
class DOFVector
: public VectorFacade<T, TraitsType::template VectorImpl>
, private Observer<event::preAdapt>
......
......@@ -28,7 +28,7 @@
#else // ISTL
#include <amdis/linearalgebra/istl/Constraints.hpp>
#include <amdis/linearalgebra/istl/ISTLRunner.hpp>
#include <amdis/linearalgebra/istl/ISTLSolver.hpp>
#include <amdis/linearalgebra/istl/PreconCreator.hpp>
#include <amdis/linearalgebra/istl/SolverCreator.hpp>
#include <amdis/linearalgebra/istl/Traits.hpp>
......
......@@ -19,7 +19,7 @@ namespace AMDiS
* \tparam T Coefficient type to store in the vector
* \tparam Traits Collection of parameter for the linear-algebra backend
**/
template <class GB, class T = double, class Traits = BackendTraits<GB,T>>
template <class GB, class T = double, class Traits = BackendTraits<GB>>
class LinearForm
: public VectorFacade<T, Traits::template VectorImpl>
{
......@@ -32,7 +32,7 @@ namespace AMDiS
using LocalView = typename GlobalBasis::LocalView;
/// The type of the elements of the DOFVector
using CoefficientType = typename Traits::CoefficientType;
using CoefficientType = T;
/// The type of the vector filled on an element with local contributions
using ElementVector = FlatVector<CoefficientType>;
......
......@@ -25,7 +25,7 @@ namespace AMDiS
using Self = ProblemInstat;
using ProblemType = ProblemStat<Traits>;
using CoefficientVector = typename ProblemType::CoefficientVector;
using SolutionVector = typename ProblemType::SolutionVector;
public:
/// Constructs a ProblemInstat with prob as its stationary problem, stored as reference.
......@@ -54,7 +54,7 @@ namespace AMDiS
ProblemType const& problemStat() const { return *problemStat_; }
/// Returns const-ref of \ref oldSolution.
std::shared_ptr<CoefficientVector const> oldSolutionVector() const
std::shared_ptr<SolutionVector const> oldSolutionVector() const
{
test_exit_dbg(bool(oldSolution_),
"OldSolution need to be created. Call initialize with INIT_UH_OLD.");
......@@ -82,7 +82,7 @@ namespace AMDiS
ProblemType* problemStat_;
/// Solution of the last timestep.
std::shared_ptr<CoefficientVector> oldSolution_;
std::shared_ptr<SolutionVector> oldSolution_;
};
......
......@@ -43,7 +43,7 @@ void ProblemInstat<Traits>::createUhOld()
if (oldSolution_)
warning("oldSolution already created\n");
else // create oldSolution
oldSolution_.reset(new CoefficientVector(*problemStat_->solutionVector()));
oldSolution_.reset(new SolutionVector(*problemStat_->solutionVector()));
}
......
......@@ -76,13 +76,17 @@ namespace AMDiS
/// Dimension of the world
static constexpr int dow = Grid::dimensionworld;
using LinearSolverTraits = BackendTraits<GlobalBasis, typename Traits::CoefficientType>;
using LinearSolverType = LinearSolverInterface<LinearSolverTraits>;
using CommunicationType = typename LinearSolverTraits::Comm;
private:
using LinAlgTraits = typename Traits::LinAlgTraits;
using Mat = typename LinAlgTraits::template MatrixImpl<typename Traits::CoefficientType>;
using Vec = typename LinAlgTraits::template VectorImpl<typename Traits::CoefficientType>;
using PartitionSet = typename LinAlgTraits::PartitionSet;
using SystemMatrix = BiLinearForm<GlobalBasis, GlobalBasis, typename Traits::CoefficientType>;
using SystemVector = LinearForm<GlobalBasis, typename Traits::CoefficientType>;
using CoefficientVector = DOFVector<GlobalBasis, typename Traits::CoefficientType>;
public:
using LinearSolver = LinearSolverInterface<Mat,Vec>;
using SystemMatrix = BiLinearForm<GlobalBasis, GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
using SystemVector = LinearForm<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
using SolutionVector = DOFVector<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
public:
/**
......@@ -334,16 +338,16 @@ namespace AMDiS
std::shared_ptr<GlobalBasis const> globalBasis() const { return globalBasis_; }
/// Return a reference to the linear solver, \ref linearSolver
std::shared_ptr<LinearSolverType> solver() { return linearSolver_; }
std::shared_ptr<LinearSolverType const> solver() const { return linearSolver_; }
std::shared_ptr<LinearSolver> solver() { return linearSolver_; }
std::shared_ptr<LinearSolver const> solver() const { return linearSolver_; }
/// Returns a reference to system-matrix, \ref systemMatrix_
std::shared_ptr<SystemMatrix> systemMatrix() { return systemMatrix_; }
std::shared_ptr<SystemMatrix const> systemMatrix() const { return systemMatrix_; }
/// Returns a reference to the solution vector, \ref solution_
std::shared_ptr<CoefficientVector> solutionVector() { return solution_; }
std::shared_ptr<CoefficientVector const> solutionVector() const { return solution_; }
std::shared_ptr<SolutionVector> solutionVector() { return solution_; }
std::shared_ptr<SolutionVector const> solutionVector() const { return solution_; }
/// Return a reference to the rhs system-vector, \ref rhs
std::shared_ptr<SystemVector> rhsVector() { return rhs_; }
......@@ -485,13 +489,13 @@ namespace AMDiS
// std::vector<Estimator*> estimator;
/// An object of the linearSolver Interface
std::shared_ptr<LinearSolverType> linearSolver_;
std::shared_ptr<LinearSolver> linearSolver_;
/// Matrix that is filled during assembling
std::shared_ptr<SystemMatrix> systemMatrix_;
/// Vector with the solution components
std::shared_ptr<CoefficientVector> solution_;
std::shared_ptr<SolutionVector> solution_;
/// Vector (load-vector) corresponding to the right-hand side
/// of the equation, filled during assembling
......
......@@ -240,7 +240,7 @@ template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
solution_ = std::make_shared<CoefficientVector>(globalBasis_);
solution_ = std::make_shared<SolutionVector>(globalBasis_);
rhs_ = std::make_shared<SystemVector>(globalBasis_);
auto localView = globalBasis_->localView();
......@@ -261,7 +261,7 @@ void ProblemStat<Traits>::createSolver()
Parameters::get(name_ + "->solver", solverName);
auto solverCreator
= named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver"));
= named(CreatorMap<LinearSolver>::getCreator(solverName, name_ + "->solver"));
linearSolver_ = solverCreator->createWithString(name_ + "->solver");
}
......@@ -291,7 +291,7 @@ void ProblemStat<Traits>::createMarker()
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
{
FileWriterCreator<CoefficientVector> creator(solution_, boundaryManager_);
FileWriterCreator<SolutionVector> creator(solution_, boundaryManager_);
filewriter_.clear();
auto localView = globalBasis_->localView();
......@@ -382,8 +382,7 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
solverInfo.setStoreMatrixData(storeMatrixData);
solution_->resize();
linearSolver_->solve(systemMatrix_->impl().matrix(), solution_->impl().vector(),
rhs_->impl().vector(), globalBasis_->communicator(), solverInfo);
linearSolver_->solve(*systemMatrix_, *solution_, *rhs_, solverInfo);
if (solverInfo.info() > 0) {
msg("solution of discrete system needed {} seconds", t.elapsed());
......@@ -525,7 +524,7 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as
msg("{} local DOFs", rhs_->localSize());
// 2. traverse grid and assemble operators on the elements
for (auto const& element : elements(gridView(), typename LinearSolverTraits::PartitionSet{})) {
for (auto const& element : elements(gridView(), PartitionSet{})) {
localView.bind(element);
if (asmMatrix)
......
......@@ -59,16 +59,20 @@ namespace AMDiS
} // end namespace Impl
/// An Exemplar for ProblemStatTraits
template <class GlobalBasisType, class T = double>
/// Wrapper around a global basis providing default traits
template <class GB, class T = double,
template <class> class TraitsImpl = BackendTraits>
struct DefaultProblemTraits
{
using GlobalBasis = GlobalBasisType;
using GlobalBasis = GB;
using CoefficientType = T;
using LinAlgTraits = TraitsImpl<GB>;
};
template <class HostGrid, class PreBasisCreator, class T = double>
/// Generator for a basis and default problem traits
template <class HostGrid, class PreBasisCreator, class T = double,
template <class> class TraitsImpl = BackendTraits>
struct DefaultBasisCreator
{
using Grid = AdaptiveGrid_t<HostGrid>;
......@@ -84,10 +88,12 @@ namespace AMDiS
return makeGlobalBasis(gridView, PreBasisCreator::create());
}
using GlobalBasis = Underlying_t<decltype(create(std::declval<GridView>()))>;
using GlobalBasis = decltype(create(std::declval<GridView>()));
using CoefficientType = T;
using LinAlgTraits = TraitsImpl<GlobalBasis>;
};
/// \brief ProblemStatTraits for a (composite) basis composed of
/// lagrange bases of different degree.
template <class Grid, int... degrees>
......
......@@ -8,13 +8,12 @@ install(FILES
LinearSolverInterface.hpp
MatrixFacade.hpp
ParallelIndexSet.hpp
PreconditionerInterface.hpp
RunnerInterface.hpp
SolverInfo.hpp
SparsityPattern.hpp
SymmetryStructure.hpp
Traits.hpp
VectorDefaultImplementation.hpp
VectorBase.hpp
VectorFacade.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra)
......
......@@ -21,16 +21,12 @@ namespace AMDiS
* solvers where the backend provides an interface, can be assigned
* by different Runner objects.
**/
template <class Traits, class Runner>
template <class Mat, class Vec, class Runner>
class LinearSolver
: public LinearSolverInterface<Traits>
: public LinearSolverInterface<Mat,Vec>
{
using Self = LinearSolver;
using Super = LinearSolverInterface<Traits>;
using Mat = typename Traits::Mat;
using Vec = typename Traits::Vec;
using Comm = typename Traits::Comm;
using Super = LinearSolverInterface<Mat,Vec>;
public:
/// A creator to be used instead of the constructor.
......@@ -49,23 +45,25 @@ namespace AMDiS
: runner_(prefix, FWD(args)...)
{}
Runner& runner() { return runner_; }
Runner const& runner() const { return runner_; }
Runner* runner() override
{
return &runner_;
}
private:
/// Implements \ref LinearSolverInterface::solveSystemImpl()
void solveImpl(Mat const& A, Vec& x, Vec const& b, Comm const& comm, SolverInfo& solverInfo) override
protected:
/// Implements \ref LinearSolverInterface::solveImpl()
void solveImpl(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) override
{
Dune::Timer t;
if (solverInfo.doCreateMatrixData()) {
// init matrix or wrap block-matrix or ...
runner_.init(A, comm);
runner_.init(A.matrix());
}
if (solverInfo.info() > 0)
msg("prepare solver needed {} seconds", t.elapsed());
int error = runner_.solve(A, x, b, solverInfo);
int error = runner_.solve(A.matrix(), x.vector(), b.vector(), solverInfo);
solverInfo.setError(error);
if (!solverInfo.doStoreMatrixData())
......@@ -73,7 +71,7 @@ namespace AMDiS
}
private:
/// redirect the implementation to a runner. Implementing a \ref RunnerInterface
/// redirect the implementation to a runner implementing a \ref RunnerInterface
Runner runner_;
};
......
#pragma once
/** \file LinearSolverInterface.h */
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/RunnerInterface.hpp>
/**
* \defgroup Solver Solver module
......@@ -15,15 +16,16 @@ namespace AMDiS
{
class SolverInfo;
template <class T, template <class> class MatrixImpl>
class MatrixFacade;
template <class T, template <class> class VectorImpl>
class VectorFacade;
/// Abstract base class for linear solvers
template <class Traits>
template <class Mat, class Vec>
class LinearSolverInterface
{
protected:
using Mat = typename Traits::Mat;
using Vec = typename Traits::Vec;
using Comm = typename Traits::Comm;
public:
/// Destructor.
virtual ~LinearSolverInterface() = default;
......@@ -38,14 +40,21 @@ namespace AMDiS
* \p x A [block-]vector for the unknown components.
* \p b A [block-]vector for the right-hand side of the linear system.
**/
void solve(Mat const& A, Vec& x, Vec const& b, Comm const& comm, SolverInfo& solverInfo)
template <class TA, class TX, class TY, template <class> class MI, template <class> class VI>
void solve(MatrixFacade<TA,MI> const& A, VectorFacade<TX,VI>& x, VectorFacade<TY,VI> const& b, SolverInfo& solverInfo)
{
solveImpl(A, x, b, comm, solverInfo);
solveImpl(A.impl(), x.impl(), b.impl(), solverInfo);
}
private:
virtual RunnerInterface<Mat,Vec>* runner()
{
error_exit("Must be implemented by derived class.");
return nullptr;
}
protected:
/// main methods that all solvers must implement
virtual void solveImpl(Mat const& A, Vec& x, Vec const& b, Comm const& comm, SolverInfo& solverInfo) = 0;
virtual void solveImpl(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) = 0;
};
} // end namespace AMDiS
......@@ -7,29 +7,28 @@ namespace AMDiS
class SolverInfo;
/// Interface for Runner / Worker types used in solver classes
template <class Traits>
template <class Mat, class Vec>
class RunnerInterface
{
protected:
using Mat = typename Traits::Mat;
using Vec = typename Traits::Vec;
using Comm = typename Traits::Comm;
using M = typename Mat::BaseMatrix;
using X = typename Vec::BaseVector;
using Y = typename Vec::BaseVector;
public:
/// virtual destructor
virtual ~RunnerInterface() = default;
/// Is called at the beginning of a solution procedure
virtual void init(Mat const& A, Comm const& comm) = 0;
virtual void init(M const& A) = 0;
/// Is called at the end of a solution procedure
virtual void exit() = 0;
/// Solve the system A*x = b
virtual int solve(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) = 0;
virtual int solve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) = 0;
/// Solve the adjoint system A^T*x = b
virtual int adjointSolve(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo)
/// Solve the system A*x = b
virtual int adjointSolve(M const& A, X& x, Y const& b, SolverInfo& solverInfo)
{
error_exit("Must be implemented by derived class.");
return 0;
......
......@@ -22,17 +22,17 @@ namespace AMDiS
* required to provide the typedefs listed here.
*
* \tparam Basis A global basis from dune-functions
* \tparam T The type of the coefficients stored in the matrix and vector
*/
template <class Basis, class T = double>
template <class Basis>
class BackendTraits
{
public:
using Mat = implementation_defined; //< The backend matrix type
using Vec = implementation_defined; //< The backend vector type
using Comm = implementation_defined; //< The communication type
template <class T>
using MatrixImpl = implementation_defined; //< The backend matrix type
template <class T>
using VectorImpl = implementation_defined; //< The backend vector type
using CoefficientType = T; //< The type of the matrix/vector entries
using SparsityPattern = implementation_defined; //< The SparsityPattern for the matrix type
using PartitionSet = Dune::Partitions::All; //< The dune partition set where to assemble operators
};
......
......@@ -14,16 +14,16 @@ namespace AMDiS
* \class AMDiS::DirectRunner
* \brief \implements RunnerInterface for the (external) direct solvers
*/
template <template <class> class Solver, class Traits>
template <class Mat, class Vec, template <class> class Solver>
class DirectRunner
: public RunnerInterface<Traits>
: public RunnerInterface<Mat,Vec>
{
using M = typename Mat::BaseMatrix;
using X = typename Vec::BaseVector;
using Y = typename Vec::BaseVector;
protected:
using Super = RunnerInterface<Traits>;
using Mat = typename Traits::Mat;
using Vec = typename Traits::Vec;
using Comm = typename Traits::Comm;
using EigenSolver = Solver<Mat>;
using EigenSolver = Solver<M>;
public:
/// Constructor.
......@@ -35,10 +35,8 @@ namespace AMDiS
}
/// Implements \ref RunnerInterface::init()
void init(Mat const& A, Comm const& comm) override
void init(M const& A) override
{
DUNE_UNUSED_PARAMETER(comm);
if (!reusePattern_ || !initialized_) {
solver_.analyzePattern(A);
initialized_ = true;
......@@ -56,13 +54,11 @@ namespace AMDiS
}
/// Implements \ref RunnerInterface::solve()
int solve(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) override
int solve(M const& A, X& x, Y const& b, SolverInfo& solverInfo)
{
DUNE_UNUSED_PARAMETER(A);
x = solver_.solve(b);
auto r = Vec(b);
Y r = b;
if (x.norm() != 0)
r -= A * x;
......
......@@ -9,15 +9,16 @@
namespace AMDiS
{
template <class IterativeSolver, class Traits>
template <class Mat, class Vec, class IterativeSolver>
class IterativeRunner
: public RunnerInterface<Traits>
: public RunnerInterface<Mat,Vec>
{
using M = typename Mat::BaseMatrix;
using X = typename Vec::BaseVector;
using Y = typename Vec::BaseVector;
using SolverCfg = SolverConfig<IterativeSolver>;
using PreconCfg = PreconConfig<typename IterativeSolver::Preconditioner>;
using Mat = typename Traits::Mat;
using Vec = typename Traits::Vec;
using Comm = typename Traits::Comm;
public:
IterativeRunner(std::string const& prefix)
......@@ -28,12 +29,9 @@ namespace AMDiS
Parameters::get(prefix + "->reuse pattern", reusePattern_);
}
/// Implements \ref RunnerInterface::init()
void init(Mat const& A, Comm const& comm) override
void init(M const& A) override
{
DUNE_UNUSED_PARAMETER(comm);
if (!reusePattern_ || !initialized_) {
solver_.analyzePattern(A);
initialized_ = true;
......@@ -51,14 +49,12 @@ namespace AMDiS
}
/// Implements \ref RunnerInterface::solve()
int solve(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) override
int solve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) override
{
DUNE_UNUSED_PARAMETER(A);
solver_.setTolerance(solverInfo.relTolerance());
x = solver_.solveWithGuess(b, x);
auto r = Vec(b);
Y r = b;
if (x.norm() != 0)
r -= A * x;
......
......@@ -20,18 +20,17 @@
namespace AMDiS
{
template <class Traits, template <class> class IterativeSolver>
template <class Mat, class Vec, template <class> class IterativeSolver>
struct IterativeSolverCreator
: public CreatorInterfaceName<LinearSolverInterface<Traits>>
: public CreatorInterfaceName<LinearSolverInterface<Mat,Vec>>
{
template <class Precon>
using SolverCreator
= typename LinearSolver<Traits,
IterativeRunner<IterativeSolver<Precon>, Traits>>::Creator;
= typename LinearSolver<Mat,Vec, IterativeRunner<Mat, Vec, IterativeSolver<Precon>>>::Creator;
using SolverBase = LinearSolverInterface<Traits>;
using Mat = typename Traits::Mat;
using Scalar = typename Mat::Scalar;
using SolverBase = LinearSolverInterface<Mat,Vec>;
using M = typename Mat::BaseMatrix;
using Scalar = typename M::Scalar;
std::unique_ptr<SolverBase> createWithString(std::string prefix) override
{
......@@ -71,11 +70,11 @@ namespace AMDiS
Parameters::get(prefix + "->precon->ordering", ordering);
if (ordering == "amd") {
using AMD = Eigen::AMDOrdering<typename Mat::StorageIndex>;
using AMD = Eigen::AMDOrdering<typename M::StorageIndex>;
return IncompleteCholesky<AMD>{}.createWithString(prefix);
}
else {
using NATURAL = Eigen::NaturalOrdering<typename Mat::StorageIndex>;
using NATURAL = Eigen::NaturalOrdering<typename M::StorageIndex>;
return IncompleteCholesky<NATURAL>{}.createWithString(prefix);
}
}
......@@ -92,37 +91,42 @@ namespace AMDiS