From b57c4f41d7058afad39f2870b6701639a4c6a352 Mon Sep 17 00:00:00 2001 From: Kilian Date: Thu, 16 Apr 2020 16:24:38 +0200 Subject: [PATCH 1/8] [linear] Add first draft of SIMPLE preconditioner --- dumux/linear/preconditioners.hh | 273 ++++++++++++++++++ .../freeflow/navierstokes/sincos/params.input | 11 +- 2 files changed, 280 insertions(+), 4 deletions(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index 08b10a7eff..5d22ccd7be 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -29,6 +29,7 @@ #include #include #include +#include #include #if HAVE_UMFPACK @@ -333,6 +334,278 @@ private: DUMUX_REGISTER_PRECONDITIONER("uzawa", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator()); + +//! See https://www.cs.umd.edu/~elman/papers/tax.pdf +template +class SeqSimple : public Dune::Preconditioner +{ + static_assert(Dumux::isMultiTypeBlockMatrix::value && M::M() == 2 && M::N() == 2, "SeqSimple expects a 2x2 MultiTypeBlockMatrix."); + static_assert(l== 1, "SeqSimple expects a block level of 1."); + + using A = std::decay_t()[Dune::Indices::_0][Dune::Indices::_0])>; + using U = std::decay_t()[Dune::Indices::_0])>; + + using Comm = Dune::Amg::SequentialInformation; + using LinearOperator = Dune::MatrixAdapter; + using Smoother = Dune::SeqSSOR; + using AMGSolverForA = Dune::Amg::AMG; + + class SimpleOperator : public Dune::LinearOperator + { + // using LocalX = std::decay_t()[Dune::Indices::_0])>; + public: + + SimpleOperator(const M& m) : m_(m) {} + + /*! \brief apply operator to x: \f$ y = A(x) \f$ + The input vector is consistent and the output must also be + consistent on the interior+border partition. + */ + void apply (const U& x, U& y) const final + { + // convenience variables + const auto& deltaP = x; + auto& rhs = y; + + using namespace Dune::Indices; + auto& A = m_[_0][_0]; + auto& B = m_[_0][_1]; + auto& C = m_[_1][_0]; + auto& D = m_[_1][_1]; + + // Apply -(-D + C * diag(A)^-1 * B) * deltaP: + + // B * deltP + // rhsTmp has different size as rhs + auto rhsTmp = U(B.N()); + B.mv(deltaP, rhsTmp); + + // diag(A)^-1 * B * deltaP + for (std::size_t i = 0; i < A.N(); ++i) + { + auto tmp = A[i][i]; + tmp.invert(); + rhsTmp[i] *= tmp; + } + + // C * diag(A)^-1 * B *deltaP + C.mv(rhsTmp, rhs); + + // (-D + C * diag(A)^-1 * B) * deltaP + D.mmv(deltaP, rhs); + + rhs *= -1.0; + } + + //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$ + void applyscaleadd (typename U::field_type alpha, const U& x, U& y) const final + { + auto tmp = y; + apply(x,tmp); + tmp *= alpha; + y += tmp; + } + + //! Category of the linear operator (see SolverCategory::Category) + Dune::SolverCategory::Category category() const final + { + return Dune::SolverCategory::sequential; + }; + + const M& m_; + + }; + +public: + //! \brief The matrix type the preconditioner is for. + using matrix_type = M; + //! \brief The domain type of the preconditioner. + using domain_type = X; + //! \brief The range type of the preconditioner. + using range_type = Y; + //! \brief The field type of the preconditioner. + using field_type = typename X::field_type; + //! \brief Scalar type underlying the field_type. + using scalar_field_type = Dune::Simd::Scalar; + + /*! + * \brief Constructor + * + * \param mat The matrix to operate on. + * \param params Collection of paramters. + */ + SeqSimple(const std::shared_ptr>& mat, const Dune::ParameterTree& params) + : matrix_(mat->getmat()) + , numIterations_(params.get("iterations")) + , relaxationFactor_(params.get("relaxation")) + , verbosity_(params.get("verbosity")) + , paramGroup_(params.get("ParameterGroup")) + , useDirectVelocitySolverForA_(getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.DirectSolverForA", false)) + { + // AMG is needed for determination of omega + if (!useDirectVelocitySolverForA_) + initAMG_(params); + + if (useDirectVelocitySolverForA_) + initUMFPack_(); + } + + /*! + * \brief Prepare the preconditioner. + */ + virtual void pre(X& x, Y& b) {} + + /*! + * \brief Apply the preconditioner + * + * \param update The update to be computed. + * \param currentDefect The current defect. + */ + virtual void apply(X& update, const Y& currentDefect) + { + using namespace Dune::Indices; + + auto& A = matrix_[_0][_0]; + auto& B = matrix_[_0][_1]; + auto& C = matrix_[_1][_0]; + auto& D = matrix_[_1][_1]; + + const auto& f = currentDefect[_0]; + const auto& g = currentDefect[_1]; + auto& u = update[_0]; + auto& p = update[_1]; + + // incorporate Dirichlet cell values + // TODO: pass Dirichlet constraint handler from outside + for (std::size_t i = 0; i < D.N(); ++i) + { + const auto& block = D[i][i]; + for (auto rowIt = block.begin(); rowIt != block.end(); ++rowIt) + if (Dune::FloatCmp::eq(rowIt->one_norm(), 1.0)) + p[i][rowIt.index()] = g[i][rowIt.index()]; + } + + for (std::size_t k = 0; k < numIterations_; ++k) + { + + // the SIMPLE algorithm + // 1. Solve A*u + B*p = f for u + auto uRhs = f; + A.mmv(u, uRhs); + B.mmv(p, uRhs); + auto uNew = u; + applySolverForA_(uNew, uRhs); + + // 2. Solve -(D + C * diag(A)^-1 * B) * deltaP = C * uNew + D*p for deltaP + using SimpleSolver = Dune::RestartedGMResSolver; + // using SimpleSolverPreconditioner = Dune::SeqSSOR; + using SimpleSolverPreconditioner = Dune::Richardson; + + Dune::InverseOperatorResult result; + SimpleSolverPreconditioner precond; + auto linearOperator = SimpleOperator(matrix_); + static const int innersolverVerbosity = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverVerbosity", 0); + static const double innersolverResidualReduction = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverResidualReduction", 1e-12); + static const int innersolverMaxIter = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverMaxIter", 1000); + static const int innersolverRestart = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverRestart", 100); + SimpleSolver solver(linearOperator, precond, innersolverResidualReduction, innersolverRestart, innersolverMaxIter, innersolverVerbosity); + + auto pRhs = p; + pRhs = 0.0; + auto deltaP = pRhs; + C.mv(uNew, pRhs); + D.mmv(p, pRhs); + + solver.apply(deltaP, pRhs, result); + + auto deltaU = U(B.N()); + B.mv(deltaP, deltaU); + + for (std::size_t i = 0; i < A.N(); ++i) + { + auto tmp = A[i][i]; + tmp.invert(); + deltaU[i] *= -tmp; + } + + // 4. Update p + auto shift = deltaP; + shift *= relaxationFactor_; + p += shift; + + // 5. Update u + u += deltaU; + } + } + + /*! + * \brief Clean up. + */ + virtual void post(X& x) {} + + //! Category of the preconditioner (see SolverCategory::Category) + virtual Dune::SolverCategory::Category category() const + { + return Dune::SolverCategory::sequential; + } + +private: + + void initAMG_(const Dune::ParameterTree& params) + { + using namespace Dune::Indices; + auto linearOperator = std::make_shared(matrix_[_0][_0]); + amgSolverForA_ = std::make_unique(linearOperator, params); + } + + void initUMFPack_() + { +#if HAVE_UMFPACK + using namespace Dune::Indices; + umfPackSolverForA_ = std::make_unique>(matrix_[_0][_0]); +#else + DUNE_THROW(Dune::InvalidStateException, "UMFPack not available. Use LinearSolver.Preconditioner.DirectVelocitySolver = false."); +#endif + } + + template + void applySolverForA_(Sol& sol, Rhs& rhs) const + { + if (useDirectVelocitySolverForA_) + { +#if HAVE_UMFPACK + Dune::InverseOperatorResult res; + umfPackSolverForA_->apply(sol, rhs, res); +#endif + } + else + { + amgSolverForA_->pre(sol, rhs); + amgSolverForA_->apply(sol, rhs); + amgSolverForA_->post(sol); + } + } + + //! \brief The matrix we operate on. + const M& matrix_; + //! \brief The number of steps to do in apply + const std::size_t numIterations_; + //! \brief The relaxation factor to use + scalar_field_type relaxationFactor_; + //! \brief The verbosity level + const int verbosity_; + + std::unique_ptr amgSolverForA_; +#if HAVE_UMFPACK + std::unique_ptr> umfPackSolverForA_; +#endif + const std::string paramGroup_; + const bool useDirectVelocitySolverForA_; +}; + +DUMUX_REGISTER_PRECONDITIONER("simple", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator()); + + } // end namespace Dumux #endif // DUMUX_LINEAR_PRECONDITIONERS_HH diff --git a/test/freeflow/navierstokes/sincos/params.input b/test/freeflow/navierstokes/sincos/params.input index 0da3caf350..d688d8fbfb 100644 --- a/test/freeflow/navierstokes/sincos/params.input +++ b/test/freeflow/navierstokes/sincos/params.input @@ -38,10 +38,13 @@ UpwindWeight = 0.5 [LinearSolver] Type = restartedflexiblegmressolver -Verbosity = 1 +Verbosity = 2 GMResRestart = 50 [LinearSolver.Preconditioner] -Type = uzawa -Verbosity = 1 -Iterations = 5 +Type = simple +Verbosity = 2 +Iterations = 2 +Relaxation = 0.3 +SimpleSolverVerbosity = 1 +SimpleSolverResidualReduction = 1e-18 -- GitLab From ea412258f8b042fc69c4585dbad5d8e9f397758a Mon Sep 17 00:00:00 2001 From: Kilian Date: Thu, 16 Apr 2020 18:08:44 +0200 Subject: [PATCH 2/8] simple fixup --- dumux/linear/preconditioners.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index 5d22ccd7be..4a6a6271e4 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -400,8 +400,8 @@ class SeqSimple : public Dune::Preconditioner //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$ void applyscaleadd (typename U::field_type alpha, const U& x, U& y) const final { - auto tmp = y; - apply(x,tmp); + auto tmp = U(y.size()); + apply(x, tmp); tmp *= alpha; y += tmp; } -- GitLab From 40a3932348e289ffa9769ea63a1f6771641d3e88 Mon Sep 17 00:00:00 2001 From: Kilian Date: Thu, 16 Apr 2020 19:22:46 +0200 Subject: [PATCH 3/8] simple fixup --- dumux/linear/preconditioners.hh | 1 - 1 file changed, 1 deletion(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index 4a6a6271e4..b9afa8db37 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -491,7 +491,6 @@ public: // the SIMPLE algorithm // 1. Solve A*u + B*p = f for u auto uRhs = f; - A.mmv(u, uRhs); B.mmv(p, uRhs); auto uNew = u; applySolverForA_(uNew, uRhs); -- GitLab From cd7ef7af4da0b463614db800bdb80ee1661b873f Mon Sep 17 00:00:00 2001 From: Kilian Date: Sat, 18 Apr 2020 10:16:47 +0200 Subject: [PATCH 4/8] fixup simple --- dumux/linear/preconditioners.hh | 81 ++++++++++++------- .../freeflow/navierstokes/sincos/params.input | 9 ++- 2 files changed, 59 insertions(+), 31 deletions(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index b9afa8db37..2540a1367f 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -352,11 +352,14 @@ class SeqSimple : public Dune::Preconditioner class SimpleOperator : public Dune::LinearOperator { - // using LocalX = std::decay_t()[Dune::Indices::_0])>; public: SimpleOperator(const M& m) : m_(m) {} + template + void setSolver(const S& solver) + { solver_ = solver; } + /*! \brief apply operator to x: \f$ y = A(x) \f$ The input vector is consistent and the output must also be consistent on the interior+border partition. @@ -376,17 +379,12 @@ class SeqSimple : public Dune::Preconditioner // Apply -(-D + C * diag(A)^-1 * B) * deltaP: // B * deltP - // rhsTmp has different size as rhs + // rhsTmp has different size than rhs auto rhsTmp = U(B.N()); B.mv(deltaP, rhsTmp); // diag(A)^-1 * B * deltaP - for (std::size_t i = 0; i < A.N(); ++i) - { - auto tmp = A[i][i]; - tmp.invert(); - rhsTmp[i] *= tmp; - } + solver_(rhsTmp, rhsTmp); // C * diag(A)^-1 * B *deltaP C.mv(rhsTmp, rhs); @@ -412,8 +410,13 @@ class SeqSimple : public Dune::Preconditioner return Dune::SolverCategory::sequential; }; + private: + + //! The system matrix const M& m_; + //! Applies the effect of A^-1 or diag(A)^-1 (depending on what is chosen in setSolver) + std::function solver_; }; public: @@ -437,7 +440,8 @@ public: SeqSimple(const std::shared_ptr>& mat, const Dune::ParameterTree& params) : matrix_(mat->getmat()) , numIterations_(params.get("iterations")) - , relaxationFactor_(params.get("relaxation")) + , relaxationFactorP_(params.get("relaxation")) + , relaxationFactorU_(relaxationFactorP_) , verbosity_(params.get("verbosity")) , paramGroup_(params.get("ParameterGroup")) , useDirectVelocitySolverForA_(getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.DirectSolverForA", false)) @@ -448,6 +452,9 @@ public: if (useDirectVelocitySolverForA_) initUMFPack_(); + + relaxationFactorP_ = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.RelaxationP", relaxationFactorP_); + relaxationFactorU_ = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.RelaxationU", relaxationFactorU_); } /*! @@ -465,7 +472,6 @@ public: { using namespace Dune::Indices; - auto& A = matrix_[_0][_0]; auto& B = matrix_[_0][_1]; auto& C = matrix_[_1][_0]; auto& D = matrix_[_1][_1]; @@ -485,34 +491,39 @@ public: p[i][rowIt.index()] = g[i][rowIt.index()]; } + // the SIMPLE algorithm for (std::size_t k = 0; k < numIterations_; ++k) { - - // the SIMPLE algorithm - // 1. Solve A*u + B*p = f for u + // 1. Solve A*uNew + B*p = f for uNew auto uRhs = f; B.mmv(p, uRhs); auto uNew = u; applySolverForA_(uNew, uRhs); - // 2. Solve -(D + C * diag(A)^-1 * B) * deltaP = C * uNew + D*p for deltaP + // 2. Solve -(D + C * diag(A)^-1 * B) * deltaP = -g + C * uNew + D*p for deltaP using SimpleSolver = Dune::RestartedGMResSolver; - // using SimpleSolverPreconditioner = Dune::SeqSSOR; using SimpleSolverPreconditioner = Dune::Richardson; Dune::InverseOperatorResult result; SimpleSolverPreconditioner precond; + auto linearOperator = SimpleOperator(matrix_); + static const bool invertSchurComplementA = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.InvertSchurComplementA", false); + if (invertSchurComplementA) + linearOperator.setSolver([&](U& a, U& b) { applySolverForA_(a, b); }); + else + linearOperator.setSolver([&](U& a, U& b) { applyInverseOfDiagonalOfA_(a); }); + static const int innersolverVerbosity = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverVerbosity", 0); static const double innersolverResidualReduction = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverResidualReduction", 1e-12); static const int innersolverMaxIter = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverMaxIter", 1000); static const int innersolverRestart = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverRestart", 100); SimpleSolver solver(linearOperator, precond, innersolverResidualReduction, innersolverRestart, innersolverMaxIter, innersolverVerbosity); - auto pRhs = p; - pRhs = 0.0; - auto deltaP = pRhs; - C.mv(uNew, pRhs); + auto pRhs = g; + pRhs *= -1.0; + auto deltaP = p; + C.umv(uNew, pRhs); D.mmv(p, pRhs); solver.apply(deltaP, pRhs, result); @@ -520,19 +531,21 @@ public: auto deltaU = U(B.N()); B.mv(deltaP, deltaU); - for (std::size_t i = 0; i < A.N(); ++i) + if (invertSchurComplementA) + applySolverForA_(deltaU, deltaU); + else { - auto tmp = A[i][i]; - tmp.invert(); - deltaU[i] *= -tmp; + applyInverseOfDiagonalOfA_(deltaU); + deltaU *= -1.0; } // 4. Update p - auto shift = deltaP; - shift *= relaxationFactor_; - p += shift; + deltaP *= relaxationFactorP_; + p += deltaP; // 5. Update u + deltaU += uNew; + deltaU *= relaxationFactorU_; u += deltaU; } } @@ -567,6 +580,19 @@ private: #endif } + void applyInverseOfDiagonalOfA_(U& x) const + { + using namespace Dune::Indices; + const auto& A = matrix_[_0][_0]; + + for (std::size_t i = 0; i < A.N(); ++i) + { + auto tmp = A[i][i]; + tmp.invert(); + x[i] *= tmp; + } + } + template void applySolverForA_(Sol& sol, Rhs& rhs) const { @@ -590,7 +616,8 @@ private: //! \brief The number of steps to do in apply const std::size_t numIterations_; //! \brief The relaxation factor to use - scalar_field_type relaxationFactor_; + scalar_field_type relaxationFactorP_; + scalar_field_type relaxationFactorU_; //! \brief The verbosity level const int verbosity_; diff --git a/test/freeflow/navierstokes/sincos/params.input b/test/freeflow/navierstokes/sincos/params.input index d688d8fbfb..8bf8865877 100644 --- a/test/freeflow/navierstokes/sincos/params.input +++ b/test/freeflow/navierstokes/sincos/params.input @@ -44,7 +44,8 @@ GMResRestart = 50 [LinearSolver.Preconditioner] Type = simple Verbosity = 2 -Iterations = 2 -Relaxation = 0.3 -SimpleSolverVerbosity = 1 -SimpleSolverResidualReduction = 1e-18 +Iterations = 5 +Relaxation = 0.5 +SimpleSolverVerbosity = 0 +SimpleSolverResidualReduction = 1e-12 +InvertSchurComplementA = true -- GitLab From df7709fcc8ccbb5438385f3912e1e53cc0dff5c5 Mon Sep 17 00:00:00 2001 From: Kilian Date: Sat, 18 Apr 2020 15:52:00 +0200 Subject: [PATCH 5/8] fixup simple --- dumux/linear/preconditioners.hh | 1 - 1 file changed, 1 deletion(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index 2540a1367f..def26157f0 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -371,7 +371,6 @@ class SeqSimple : public Dune::Preconditioner auto& rhs = y; using namespace Dune::Indices; - auto& A = m_[_0][_0]; auto& B = m_[_0][_1]; auto& C = m_[_1][_0]; auto& D = m_[_1][_1]; -- GitLab From 8b17b9e789666d1870e3b21ccee5dce550b73005 Mon Sep 17 00:00:00 2001 From: Kilian Date: Thu, 23 Apr 2020 16:38:44 +0200 Subject: [PATCH 6/8] fixup SIMPLE --- dumux/linear/preconditioners.hh | 87 ++++++++++++++++++++------------- 1 file changed, 53 insertions(+), 34 deletions(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index def26157f0..5ca2a502ca 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -349,6 +349,7 @@ class SeqSimple : public Dune::Preconditioner using LinearOperator = Dune::MatrixAdapter; using Smoother = Dune::SeqSSOR; using AMGSolverForA = Dune::Amg::AMG; + using SchurComplementSolver = Dune::RestartedGMResSolver; class SimpleOperator : public Dune::LinearOperator { @@ -375,7 +376,7 @@ class SeqSimple : public Dune::Preconditioner auto& C = m_[_1][_0]; auto& D = m_[_1][_1]; - // Apply -(-D + C * diag(A)^-1 * B) * deltaP: + // Apply the Schur complement (-D + C * diag(A)^-1 * B) * deltaP: // B * deltP // rhsTmp has different size than rhs @@ -390,8 +391,6 @@ class SeqSimple : public Dune::Preconditioner // (-D + C * diag(A)^-1 * B) * deltaP D.mmv(deltaP, rhs); - - rhs *= -1.0; } //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$ @@ -444,6 +443,8 @@ public: , verbosity_(params.get("verbosity")) , paramGroup_(params.get("ParameterGroup")) , useDirectVelocitySolverForA_(getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.DirectSolverForA", false)) + , useDiagonalInSchurComplement_(getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.UseDiagonalInSchurComplement", false)) + , useDiagonalInUpdate_(getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.UseDiagonalInUpdate", false)) { // AMG is needed for determination of omega if (!useDirectVelocitySolverForA_) @@ -454,6 +455,8 @@ public: relaxationFactorP_ = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.RelaxationP", relaxationFactorP_); relaxationFactorU_ = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.RelaxationU", relaxationFactorU_); + + initSchurComplementSolver_(); } /*! @@ -471,6 +474,7 @@ public: { using namespace Dune::Indices; + auto& A = matrix_[_0][_0]; auto& B = matrix_[_0][_1]; auto& C = matrix_[_1][_0]; auto& D = matrix_[_1][_1]; @@ -480,8 +484,10 @@ public: auto& u = update[_0]; auto& p = update[_1]; + static const bool dirichletHack = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.DirichletHack", false); // incorporate Dirichlet cell values // TODO: pass Dirichlet constraint handler from outside + if (dirichletHack) for (std::size_t i = 0; i < D.N(); ++i) { const auto& block = D[i][i]; @@ -499,44 +505,26 @@ public: auto uNew = u; applySolverForA_(uNew, uRhs); - // 2. Solve -(D + C * diag(A)^-1 * B) * deltaP = -g + C * uNew + D*p for deltaP - using SimpleSolver = Dune::RestartedGMResSolver; - using SimpleSolverPreconditioner = Dune::Richardson; - - Dune::InverseOperatorResult result; - SimpleSolverPreconditioner precond; - - auto linearOperator = SimpleOperator(matrix_); - static const bool invertSchurComplementA = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.InvertSchurComplementA", false); - if (invertSchurComplementA) - linearOperator.setSolver([&](U& a, U& b) { applySolverForA_(a, b); }); - else - linearOperator.setSolver([&](U& a, U& b) { applyInverseOfDiagonalOfA_(a); }); - - static const int innersolverVerbosity = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverVerbosity", 0); - static const double innersolverResidualReduction = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverResidualReduction", 1e-12); - static const int innersolverMaxIter = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverMaxIter", 1000); - static const int innersolverRestart = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverRestart", 100); - SimpleSolver solver(linearOperator, precond, innersolverResidualReduction, innersolverRestart, innersolverMaxIter, innersolverVerbosity); - + // 2. Solve (D + C * diag(A)^-1 * B) * deltaP = -g + C * uNew + D*p for deltaP auto pRhs = g; pRhs *= -1.0; - auto deltaP = p; C.umv(uNew, pRhs); - D.mmv(p, pRhs); - - solver.apply(deltaP, pRhs, result); + D.umv(p, pRhs); + auto deltaP = p; + Dune::InverseOperatorResult result; + schurComplementSolver_->apply(deltaP, pRhs, result); + // 3. Calculate correction of u: + // deltaU = (-diag(A^-1) * B) * deltaP auto deltaU = U(B.N()); B.mv(deltaP, deltaU); - if (invertSchurComplementA) - applySolverForA_(deltaU, deltaU); - else - { + if (useDiagonalInUpdate_) applyInverseOfDiagonalOfA_(deltaU); - deltaU *= -1.0; - } + else + applySolverForA_(deltaU, deltaU); + + deltaU *= -1.0; // 4. Update p deltaP *= relaxationFactorP_; @@ -562,6 +550,13 @@ public: private: + // void apply + + void applyStandardSIMPLE_() + { + + } + void initAMG_(const Dune::ParameterTree& params) { using namespace Dune::Indices; @@ -579,6 +574,25 @@ private: #endif } + void initSchurComplementSolver_() + { + auto linearOperator = std::make_shared(matrix_); + if (useDiagonalInSchurComplement_) + linearOperator->setSolver([&](U& a, U& b) { applyInverseOfDiagonalOfA_(a); }); + else + linearOperator->setSolver([&](U& a, U& b) { applySolverForA_(a, b); }); + + Dune::ParameterTree tree; + tree["verbose"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverVerbosity", "0"); + tree["reduction"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverResidualReduction", "1e-12"); + tree["maxit"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverMaxIter", "1000"); + tree["restart"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverRestart", "100"); + + schurComplementPreconditioner_ = std::make_shared>(); + // schurComplementSolver_ = std::make_unique(linearOperator, schurComplementPreconditioner_, tree); + schurComplementSolver_ = std::make_unique(linearOperator, std::make_shared>(), tree); + } + void applyInverseOfDiagonalOfA_(U& x) const { using namespace Dune::Indices; @@ -599,7 +613,8 @@ private: { #if HAVE_UMFPACK Dune::InverseOperatorResult res; - umfPackSolverForA_->apply(sol, rhs, res); + auto rhsTmp = rhs; + umfPackSolverForA_->apply(sol, rhsTmp, res); #endif } else @@ -621,11 +636,15 @@ private: const int verbosity_; std::unique_ptr amgSolverForA_; + std::unique_ptr schurComplementSolver_; + std::shared_ptr> schurComplementPreconditioner_; #if HAVE_UMFPACK std::unique_ptr> umfPackSolverForA_; #endif const std::string paramGroup_; const bool useDirectVelocitySolverForA_; + const bool useDiagonalInSchurComplement_; + const bool useDiagonalInUpdate_; }; DUMUX_REGISTER_PRECONDITIONER("simple", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator()); -- GitLab From b40a499c03666689ddaec8ef48b852a9b600477e Mon Sep 17 00:00:00 2001 From: Kilian Date: Fri, 24 Apr 2020 09:33:01 +0200 Subject: [PATCH 7/8] [simple] Try Jacobi preconditioner on Schur complement --- dumux/linear/preconditioners.hh | 108 +++++++++++++++++++++++++++++--- 1 file changed, 100 insertions(+), 8 deletions(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index 5ca2a502ca..89d290c919 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -335,6 +335,90 @@ private: DUMUX_REGISTER_PRECONDITIONER("uzawa", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator()); +template +class SeqJacobiMatrixFree : public Dune::Preconditioner +{ +public: + //! \brief The domain type of the preconditioner. + using domain_type = X; + //! \brief The range type of the preconditioner. + using range_type = Y; + //! \brief The field type of the preconditioner. + using field_type = typename X::field_type; + //! \brief Scalar type underlying the field_type. + using scalar_field_type = Dune::Simd::Scalar; + + SeqJacobiMatrixFree(std::shared_ptr> op, const Dune::ParameterTree& params) + : op_(op) + {} + + /*! + * \brief Prepare the preconditioner. + */ + virtual void pre(X& x, Y& b) {} + + /*! + * \brief Apply the preconditioner + * + * \param update The update to be computed. + * \param currentDefect The current defect. + */ + virtual void apply(X& update, const Y& currentDefect) + { + + if (inverseDiagonalElements_.size() == 0) + init_(update, currentDefect); + + + for (std::size_t i = 0; i < update.size(); ++i) + { + auto tmp = currentDefect[i]; + tmp *= inverseDiagonalElements_[i]; + update[i] = tmp; + } + } + + /*! + * \brief Clean up. + */ + virtual void post(X& x) {} + + //! Category of the preconditioner (see SolverCategory::Category) + virtual Dune::SolverCategory::Category category() const + { + return Dune::SolverCategory::sequential; + } + + +private: + + void init_(X& update, const Y& currentDefect) + { + Dune::Timer timer; + auto tmp = X(update.size()); + inverseDiagonalElements_ = tmp; + auto res = tmp; + + for (std::size_t i = 0; i < inverseDiagonalElements_.size(); ++i) + { + tmp = 0; + res = 0; + tmp[i] = 1.0; + op_->apply(tmp, res); + inverseDiagonalElements_[i] = 1.0/res[i]; + } + + timer.stop(); + + std::cout << "initializing Jac precond took " << timer.elapsed() << " seconds" << std::endl; + } + + + std::shared_ptr> op_; + X inverseDiagonalElements_; +}; + + //! See https://www.cs.umd.edu/~elman/papers/tax.pdf template class SeqSimple : public Dune::Preconditioner @@ -582,15 +666,23 @@ private: else linearOperator->setSolver([&](U& a, U& b) { applySolverForA_(a, b); }); - Dune::ParameterTree tree; - tree["verbose"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverVerbosity", "0"); - tree["reduction"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverResidualReduction", "1e-12"); - tree["maxit"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverMaxIter", "1000"); - tree["restart"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverRestart", "100"); + Dune::ParameterTree treeSchurComplementSolver; + treeSchurComplementSolver["verbose"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverVerbosity", "0"); + treeSchurComplementSolver["reduction"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverResidualReduction", "1e-12"); + treeSchurComplementSolver["maxit"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverMaxIter", "1000"); + treeSchurComplementSolver["restart"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverRestart", "100"); + + static const bool usePreconditionerForSchurComplement = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.UsePreconditionerForSchurComplement", false); + if (usePreconditionerForSchurComplement) + { + Dune::ParameterTree treePreconditioner; + treePreconditioner["verbose"] = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverVerbosity", "0"); + schurComplementPreconditioner_ = std::make_shared>(linearOperator, treePreconditioner); + } + else + schurComplementPreconditioner_ = std::make_shared>(); - schurComplementPreconditioner_ = std::make_shared>(); - // schurComplementSolver_ = std::make_unique(linearOperator, schurComplementPreconditioner_, tree); - schurComplementSolver_ = std::make_unique(linearOperator, std::make_shared>(), tree); + schurComplementSolver_ = std::make_unique(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver); } void applyInverseOfDiagonalOfA_(U& x) const -- GitLab From ec6d0873039c78e7336b1ed19b141836d3578c22 Mon Sep 17 00:00:00 2001 From: Kilian Date: Tue, 28 Apr 2020 16:31:37 +0200 Subject: [PATCH 8/8] fixup simple --- dumux/linear/preconditioners.hh | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh index 89d290c919..bf51d155a2 100644 --- a/dumux/linear/preconditioners.hh +++ b/dumux/linear/preconditioners.hh @@ -433,7 +433,6 @@ class SeqSimple : public Dune::Preconditioner using LinearOperator = Dune::MatrixAdapter; using Smoother = Dune::SeqSSOR; using AMGSolverForA = Dune::Amg::AMG; - using SchurComplementSolver = Dune::RestartedGMResSolver; class SimpleOperator : public Dune::LinearOperator { @@ -482,8 +481,7 @@ class SeqSimple : public Dune::Preconditioner { auto tmp = U(y.size()); apply(x, tmp); - tmp *= alpha; - y += tmp; + y.axpy(alpha, tmp); } //! Category of the linear operator (see SolverCategory::Category) @@ -558,7 +556,6 @@ public: { using namespace Dune::Indices; - auto& A = matrix_[_0][_0]; auto& B = matrix_[_0][_1]; auto& C = matrix_[_1][_0]; auto& D = matrix_[_1][_1]; @@ -682,7 +679,16 @@ private: else schurComplementPreconditioner_ = std::make_shared>(); - schurComplementSolver_ = std::make_unique(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver); + static const auto schurComplementSolverType = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.SimpleSolverType", "restartedgmressolver"); + + if (schurComplementSolverType == "restartedgmressolver") + schurComplementSolver_ = std::make_unique>(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver); + else if (schurComplementSolverType == "bicgstabsolver") + schurComplementSolver_ = std::make_unique>(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver); + else if (schurComplementSolverType == "cgsolver") + schurComplementSolver_ = std::make_unique>(linearOperator, schurComplementPreconditioner_, treeSchurComplementSolver); + else + DUNE_THROW(Dune::InvalidStateException, schurComplementSolverType << "not supported. Use restartedgmressolver, bicgstabsolver or cgsolver"); } void applyInverseOfDiagonalOfA_(U& x) const @@ -728,7 +734,7 @@ private: const int verbosity_; std::unique_ptr amgSolverForA_; - std::unique_ptr schurComplementSolver_; + std::unique_ptr> schurComplementSolver_; std::shared_ptr> schurComplementPreconditioner_; #if HAVE_UMFPACK std::unique_ptr> umfPackSolverForA_; -- GitLab