diff --git a/dumux/linear/linearsolverproperties.hh b/dumux/linear/linearsolverproperties.hh index 32bdbbb1748364ca423829baaea3873c6b47ef51..66ef370c1e73a8497ade953baf655e49189b0bca 100644 --- a/dumux/linear/linearsolverproperties.hh +++ b/dumux/linear/linearsolverproperties.hh @@ -60,6 +60,9 @@ NEW_PROP_TAG(LinearSolverPreconditionerRelaxation); //! number of preconditioner iterations per solver iteration NEW_PROP_TAG(LinearSolverPreconditionerIterations); +//! Verbosity level of the preconditioner +NEW_PROP_TAG(LinearSolverPreconditionerVerbosity); + //! Block level depth for the preconditioner // Set this to more than one if the matrix to solve is nested multiple times // e.g. for Dune::MultiTypeBlockMatrix'es. @@ -85,6 +88,9 @@ SET_SCALAR_PROP(LinearSolverTypeTag, LinearSolverPreconditionerRelaxation, 1.0); //! set the preconditioner iterations to 1 by default SET_INT_PROP(LinearSolverTypeTag, LinearSolverPreconditionerIterations, 1); +//! set the preconditioner verbosity level to 0 by default +SET_INT_PROP(LinearSolverTypeTag, LinearSolverPreconditionerVerbosity, 0); + //! set the block level to 1, suitable for e.g. a simple Dune::BCRSMatrix. SET_INT_PROP(LinearSolverTypeTag, LinearSolverPreconditionerBlockLevel, 1); diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh index a729428ff8a47bad8bc9424b28b89c6881dddbf1..0362874ebfabc1d3c6c27652038334362d2ed393 100644 --- a/dumux/linear/seqsolverbackend.hh +++ b/dumux/linear/seqsolverbackend.hh @@ -32,6 +32,10 @@ #include #include +#include +#include +#include + namespace Dumux { @@ -644,6 +648,328 @@ class ILU0BiCGSTABBackend : public ILU0SolverBackend } }; + + +/*! + \brief Schur complement inverse operator. +*/ +template +class SchurComplement : public Dune::LinearOperator +{ +public: + // export types + typedef DType matrix_type; + typedef X domain_type; + typedef Y range_type; + typedef typename X::field_type field_type; + + //! constructor: just store a reference to a matrix + explicit SchurComplement (const AType& A, const BType& B, + const CType& C, const DType& D, + const std::shared_ptr& invA) + : A_(A), B_(B), C_(C), D_(D), invA_(invA) {} + + //! apply operator to x: \f$ y = A(x) \f$ + virtual void apply (const X& x, Y& y) const + { + // apply B, note x and aTmp1 have different size (Bx) + X aTmp1(A_.N()); + aTmp1 = 0.0; + B_.mv(x, aTmp1); + + // apply A^-1 (A^-1Bx) + auto aTmp2 = aTmp1; + + // if we use the exact Schur complement the eigenvalues are 1 and GMRes converges in 2 steps + // ... a single V-cycle of AMG instead + Dune::InverseOperatorResult result; + invA_->apply(aTmp2, aTmp1, result); + // invA_.pre(aTmp2, aTmp1); + // invA_.apply(aTmp2, aTmp1); + // invA_.post(aTmp2); + + // apply C (CA^-1Bx) + C_.mv(aTmp2, y); + // switch signs + y *= -1.0; + + // add Dx (Dx - CA^-1Bx) -> the full Schur complement + D_.umv(x, y); + } + + //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$ + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const + { + auto tmp = y; + this->apply(x, tmp); + y.axpy(alpha, tmp); + } + + //! Category of the preconditioner (see SolverCategory::Category) + virtual Dune::SolverCategory::Category category() const + { + return Dune::SolverCategory::sequential; + } + + private: + const AType& A_; + const BType& B_; + const CType& C_; + const DType& D_; + std::shared_ptr invA_; +}; + +/*! + \brief Schur complement inverse operator. +*/ +template +class SchurApproximate : public Dune::LinearOperator +{ +public: + // export types + typedef DType matrix_type; + typedef X domain_type; + typedef Y range_type; + typedef typename X::field_type field_type; + + //! constructor: just store a reference to a matrix + explicit SchurApproximate (const AType& A, const BType& B, + const CType& C, const DType& D) + : A_(A), B_(B), C_(C), D_(D) {} + + //! apply operator to x: \f$ y = A(x) \f$ + virtual void apply (const X& x, Y& y) const + { + auto viscosity = 1.0; + auto theta = 1.0; + auto rho = 1.0; + + // theta*rho*C*B*x + X tmp(A_.N()); + tmp = 0.0; + B_.mv(x, tmp); + C_.mv(tmp, y); + y *= theta*rho; + + // + D*x + D_.umv(x, y); + + // + I*mu*x + y.axpy(viscosity, x); + } + + //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$ + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const + { + auto tmp = y; + this->apply(x, tmp); + y.axpy(alpha, tmp); + } + + //! Category of the preconditioner (see SolverCategory::Category) + virtual Dune::SolverCategory::Category category() const + { + return Dune::SolverCategory::sequential; + } + + private: + const AType& A_; + const BType& B_; + const CType& C_; + const DType& D_; +}; + +/*! \brief The schur complement preconditioner + + \tparam M The matrix type to operate on + \tparam X Type of the update + \tparam Y Type of the defect + \tparam l The block level to invert. Default is 1 +*/ +template +class SchurComplementPreconditioner : public Dune::Preconditioner +{ +public: + //! \brief The matrix type the preconditioner is for. + typedef M matrix_type; + //! \brief The domain type of the preconditioner. + typedef X domain_type; + //! \brief The range type of the preconditioner. + typedef Y range_type; + //! \brief The field type of the preconditioner. + typedef typename X::field_type field_type; + + /*! \brief Constructor. + + Constructor gets all parameters to operate the prec. + \param A The matrix to operate on. + \note this is an upper triangular preconditioner of the form + | A B |^-1 + | 0 S | + */ + SchurComplementPreconditioner (const M& matrix, + const std::shared_ptr& aInvOp) + : M_(matrix), aInvOp_(aInvOp) + {} + + virtual void pre (X& x, Y& b) + { + DUNE_UNUSED_PARAMETER(x); + DUNE_UNUSED_PARAMETER(b); + } + + /*! + \brief Apply the preconditioner. + + \copydoc Preconditioner::apply(X&,const Y&) + */ + virtual void apply (X& v, const Y& b) + { + using namespace Dune::Indices; + using VVector = Dune::BlockVector; + using PVector = Dune::BlockVector; + + using BType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockCCToFace; + using AType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockFaceToFace; + using CType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockFaceToCC; + using DType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockCCToCC; + + // Currently Navier-Stokes system is assembled the other way around + // Consider flipping indices + const AType& A = M_[_1][_1]; + const BType& B = M_[_1][_0]; + const CType& C = M_[_0][_1]; + const DType& D = M_[_0][_0]; + + const VVector& f = b[_1]; + const PVector& g = b[_0]; + + // first solve pressure by applying the Schur inverse operator using a GMRes solver + static const int schurMaxIter = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, SchurIterations); + static const int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, PreconditionerVerbosity); + static const double schurResred = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, SchurResidualReduction); + static const double velResred = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, VelocityResidualReduction); + using Identity = Dune::Richardson; + Identity identity(1.0); + // SchurComplement schurOperator(A, B, C, D, aInvOp_); + SchurApproximate schurOperator(A, B, C, D); + Dune::RestartedGMResSolver invOpS(schurOperator, identity, schurResred*g.two_norm(), 100, schurMaxIter, verbosity); + // Dune::BiCGSTABSolver invOpS(schurOperator, identity, schurResred*g.two_norm(), schurMaxIter, verbosity); + + // apply S^-1 to solve pressure block + PVector gTmp(g); + Dune::InverseOperatorResult result; + invOpS.apply(v[_0], gTmp, result); + + // prepare rhs for application of A^-1 (f - Bp) + VVector vrhs(f); + B.mmv(v[_0], vrhs); + + // then solve for velocity block v using the action of A^-1 + // approximate the inverse operator of A or use direct solver + // here e.g. GMRes with AMG preconditioner + Dune::InverseOperatorResult resultTmp; + aInvOp_->apply(v[_1], vrhs, velResred*vrhs.two_norm(), resultTmp); + } + + virtual void post (X& x) + { DUNE_UNUSED_PARAMETER(x); } + + //! Category of the preconditioner (see SolverCategory::Category) + virtual Dune::SolverCategory::Category category() const + { return Dune::SolverCategory::sequential; } + + private: + //! \brief The matrix we operate on. + const M& M_; + std::shared_ptr aInvOp_; + }; + +template +class SchurComplementSolver +{ + using Problem = typename GET_PROP_TYPE(TypeTag, Problem); + +public: + SchurComplementSolver() = default; + SchurComplementSolver(const Problem& problem) {} + + // expects a system as a multi-type matrix + // | A B | + // | C D | + + // Solve saddle-point problem using a Schur complement based preconditioner + template + bool solve(const Matrix& A, Vector& x, const Vector& b) + { + int verbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, Verbosity); + const int maxIter = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, MaxIterations); + const double residReduction = GET_PARAM_FROM_GROUP(TypeTag, double, LinearSolver, ResidualReduction); + const int restartGMRes = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, GMResRestart); + + Vector bTmp(b); + + using MatrixAdapter = Dune::MatrixAdapter; + MatrixAdapter operatorA(A); + + // LU decompose vel sub matrix and hand it to the preconditioner + using namespace Dune::Indices; + using VelMatrixType = typename GET_PROP(TypeTag, JacobianMatrix)::MatrixBlockFaceToFace; + using VVector = Dune::BlockVector; + // using InnerSolver = Dune::UMFPack; + // auto invOpVel = std::make_shared(A[_1][_1], false); + // using InnerPreconditioner = Dune::SeqJac; + // auto innerPrec = std::make_shared(A[_1][_1], 1, 1.0); + + // or use an AMG approximation instead instead + using Comm = Dune::Amg::SequentialInformation; + using LinearOperator = Dune::MatrixAdapter; + using ScalarProduct = Dune::SeqScalarProduct; + using Smoother = Dune::SeqSSOR; + auto comm = std::make_shared(); + auto fop = std::make_shared(A[_1][_1]); + auto sp = std::make_shared(); + + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + using Criterion = Dune::Amg::CoarsenCriterion >; + Dune::Amg::Parameters params(15, 2000, 1.2, 1.6, Dune::Amg::atOnceAccu); + params.setDefaultValuesIsotropic(GET_PROP_TYPE(TypeTag, GridView)::Traits::Grid::dimension); + params.setDebugLevel(verbosity); + params.setGamma(1); // one V-cycle + Criterion criterion(params); + SmootherArgs smootherArgs; + smootherArgs.iterations = 1; + smootherArgs.relaxationFactor = 1.0; + + using InnerPreconditioner = Dune::Amg::AMG; + auto innerPrec = std::make_shared(*fop, criterion, smootherArgs, *comm); + static const int precVerbosity = GET_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, PreconditionerVerbosity); + static const int precMaxIter = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, LinearSolver, VelocityIterations); + using InnerSolver = Dune::RestartedGMResSolver; + auto invOpVel = std::make_shared(*fop, *sp, *innerPrec, 1.0, 100, precMaxIter, precVerbosity); + + using Preconditioner = SchurComplementPreconditioner; + Preconditioner schurPrecond(A, invOpVel); + + // using Solver = Dune::MINRESSolver; + using Solver = Dune::RestartedfGMResSolver; + // using Solver = Dune::BiCGSTABSolver; + Solver solver(operatorA, schurPrecond, residReduction, restartGMRes, maxIter, verbosity); + // Solver solver(operatorA, schurPrecond, residReduction, maxIter, verbosity); + solver.apply(x, bTmp, result_); + + return result_.converged; + } + + const Dune::InverseOperatorResult& result() const + { + return result_; + } + +private: + Dune::InverseOperatorResult result_; +}; + /*! * \ingroup Linear * \brief Sequential ILU(0)-preconditioned CG solver. diff --git a/dumux/porousmediumflow/sequential/pressureproperties.hh b/dumux/porousmediumflow/sequential/pressureproperties.hh index 668ed17b82bded003dcf8e9bea27ce2eb908a163..d48db6c3ebc9c04a4e3c88b5aad9639ed0f381e3 100644 --- a/dumux/porousmediumflow/sequential/pressureproperties.hh +++ b/dumux/porousmediumflow/sequential/pressureproperties.hh @@ -70,6 +70,9 @@ NEW_PROP_TAG(Velocity); namespace Dumux { +template +class ILU0BiCGSTABBackend; + namespace Properties { //! Faces are only regarded from one side and not from both cells diff --git a/test/freeflow/staggered/doneatestproblem.hh b/test/freeflow/staggered/doneatestproblem.hh index ba22ea2b507d174444be45c8317a629f537f1e55..43386a1c693b151e131abbe9503df79e25144663 100644 --- a/test/freeflow/staggered/doneatestproblem.hh +++ b/test/freeflow/staggered/doneatestproblem.hh @@ -64,6 +64,8 @@ SET_TYPE_PROP(DoneaTestProblem, Grid, Dune::YaspGrid<2>); // Set the problem property SET_TYPE_PROP(DoneaTestProblem, Problem, Dumux::DoneaTestProblem ); +SET_TYPE_PROP(DoneaTestProblem, LinearSolver, SchurComplementSolver ); + SET_BOOL_PROP(DoneaTestProblem, EnableFVGridGeometryCache, true); SET_BOOL_PROP(DoneaTestProblem, EnableGlobalFluxVariablesCache, true); diff --git a/test/freeflow/staggered/test_donea.input b/test/freeflow/staggered/test_donea.input index 59941d981064392712132f82fdbf37cf8aee7950..a4d64837cd91f7db837ba6bffb7e4f049c9ac6a0 100644 --- a/test/freeflow/staggered/test_donea.input +++ b/test/freeflow/staggered/test_donea.input @@ -1,6 +1,6 @@ [TimeManager] DtInitial = 1 # [s] -TEnd = 2 # [s] +TEnd = 1 # [s] [Grid] UpperRight = 1 1 @@ -15,3 +15,14 @@ PrintL2Error = false [ Newton ] MaxSteps = 10 MaxRelativeShift = 1e-5 + +[LinearSolver] +GMResRestart = 100 +ResidualReduction = 1e-7 +Verbosity = 0 +MaxIterations = 2000 +PreconditionerVerbosity = 0 +SchurIterations = 1 +VelocityIterations = 3 +SchurResidualReduction = 1e-10 +VelocityResidualReduction = 1e-10