diff --git a/dumux/linear/istlsolverfactorybackend.hh b/dumux/linear/istlsolverfactorybackend.hh index 40e6709e10ab259b0d45db448c6629d108a544ca..8bac63d6edaea41fa7be519390f71544fb1ade90 100644 --- a/dumux/linear/istlsolverfactorybackend.hh +++ b/dumux/linear/istlsolverfactorybackend.hh @@ -90,173 +90,6 @@ void initSolverFactories() * to choose the solver and preconditioner at runtime. * \note the solvers are configured via the input file */ -namespace Detail { - -template <class LinearSolverTraits> -class OldIstlSolverFactoryBackend : public LinearSolver -{ -public: - - /*! - * \brief Update the solver after grid adaption - * - * \param gridView the grid view on which we are performing the multi-grid - * \param dofMapper an index mapper for dof entities - */ - void updateAfterGridAdaption(const typename LinearSolverTraits::GridView& gridView, - const typename LinearSolverTraits::DofMapper& dofMapper) - { -#if HAVE_MPI - if (isParallel_) - parallelHelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper); -#endif - } - - /*! - * \brief Solve a linear system. - * - * \param A the matrix - * \param x the seeked solution vector, containing the initial solution upon entry - * \param b the right hand side vector - */ - template<class Matrix, class Vector> - bool solve(Matrix& A, Vector& x, Vector& b) - { -#if HAVE_MPI - solveSequentialOrParallel_(A, x, b); -#else - solveSequential_(A, x, b); -#endif - firstCall_ = false; - return result_.converged; - } - - const Dune::InverseOperatorResult& result() const - { - return result_; - } - - const std::string& name() const - { - return name_; - } - -private: - - void initializeParameters_() - { - params_ = Dumux::LinearSolverParameters<LinearSolverTraits>::createParameterTree(paramGroup_); - checkMandatoryParameters_(); - name_ = params_.get<std::string>("preconditioner.type") + "-preconditioned " + params_.get<std::string>("type"); - if (params_.get<int>("verbose", 0) > 0) - std::cout << "Initialized linear solver of type: " << name_ << std::endl; - } - - void checkMandatoryParameters_() - { - if (!params_.hasKey("type")) - DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Type\" parameter to select the solver"); - - if (!params_.hasKey("preconditioner.type")) - DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Preconditioner.Type\" parameter to select the preconditioner"); - } - -#if HAVE_MPI - template<class Matrix, class Vector> - void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b) - { - // Dune::MultiTypeBlockMatrix does not provide a ColIterator which is needed by Dune::NonoverlappingSchwarzOperator. - // We therefore can only solve these types of systems sequentially. - // TODO: This can be adapted once the situation in Dune ISTL changes. - if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockMatrix<Matrix>::value) - { - if (isParallel_) - { - if (LinearSolverTraits::isNonOverlapping(parallelHelper_->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); - } - else - { - solveSequential_(A, x, b); - } - } - - template<class ParallelTraits, class Matrix, class Vector> - void solveParallel_(Matrix& A, Vector& x, Vector& b) - { - using Comm = typename ParallelTraits::Comm; - using LinearOperator = typename ParallelTraits::LinearOperator; - using ScalarProduct = typename ParallelTraits::ScalarProduct; - - if (firstCall_) - initSolverFactories<Matrix, LinearOperator>(); - - std::shared_ptr<Comm> comm; - std::shared_ptr<LinearOperator> linearOperator; - std::shared_ptr<ScalarProduct> scalarProduct; - prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, comm, linearOperator, scalarProduct, *parallelHelper_); - - // construct solver - auto solver = getSolverFromFactory_(linearOperator); - - // solve linear system - solver->apply(x, b, result_); - } -#endif // HAVE_MPI - - template<class Matrix, class Vector> - void solveSequential_(Matrix& A, Vector& x, Vector& b) - { - // construct linear operator - using Traits = typename LinearSolverTraits::template Sequential<Matrix, Vector>; - using LinearOperator = typename Traits::LinearOperator; - auto linearOperator = std::make_shared<LinearOperator>(A); - - if (firstCall_) - initSolverFactories<Matrix, LinearOperator>(); - - // construct solver - auto solver = getSolverFromFactory_(linearOperator); - - // solve linear system - solver->apply(x, b, result_); - } - - template<class LinearOperator> - auto getSolverFromFactory_(std::shared_ptr<LinearOperator>& fop) - { - try { return Dune::getSolverFromFactory(fop, params_); } - catch(Dune::Exception& e) - { - std::cerr << "Could not create solver with factory" << std::endl; - std::cerr << e.what() << std::endl; - throw e; - } - } - - const std::string paramGroup_; -#if HAVE_MPI - std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> parallelHelper_; -#endif - bool isParallel_ = false; - bool firstCall_; - - Dune::InverseOperatorResult result_; - Dune::ParameterTree params_; - std::string name_; -}; - template <class LinearSolverTraits, class LinearAlgebraTraits> class IstlSolverFactoryBackend : public LinearSolver { @@ -480,32 +313,6 @@ private: std::string name_; }; -template<int i> -struct IstlSolverFactoryBackendImplHelper -{ - template<class ...T> - using type = IstlSolverFactoryBackend<T...>; -}; - -template<> -struct IstlSolverFactoryBackendImplHelper<1> -{ - template<class T> - using type = OldIstlSolverFactoryBackend<T>; -}; - -} // end namespace Detail - -/*! - * \ingroup Linear - * \brief A linear solver using the dune-istl solver factory - * to choose the solver and preconditioner at runtime. - * \note the solvers are configured via the input file - * \note The template arguments are LinearSolverTraits, LinearAlgebraTraits - */ -template<class ...T> -using IstlSolverFactoryBackend = typename Detail::IstlSolverFactoryBackendImplHelper<sizeof...(T)>::template type<T...>; - } // end namespace Dumux #endif