diff --git a/doc/handbook/dumux-handbook.bib b/doc/handbook/dumux-handbook.bib index 88860ce45fe6aca8bff3da2bc3c0d3e7b22764e9..9b5588c9ac46388b98565dc4233d2b0c85e0192e 100644 --- a/doc/handbook/dumux-handbook.bib +++ b/doc/handbook/dumux-handbook.bib @@ -1882,3 +1882,28 @@ year={1999} pages={37}, year={1964} } + +@article{benzi2005, + doi = {10.1017/s0962492904000212}, + url = {https://doi.org/10.1017/s0962492904000212}, + year = {2005}, + publisher = {Cambridge University Press ({CUP})}, + volume = {14}, + pages = {1--137}, + author = {Michele Benzi and Gene H. Golub and J\"{o}rg Liesen}, + title = {Numerical solution of saddle point problems}, + journal = {Acta Numerica} +} + +@article{ho2017, + doi = {10.1137/16m1076770}, + url = {https://doi.org/10.1137/16m1076770}, + year = {2017}, + publisher = {Society for Industrial {\&} Applied Mathematics ({SIAM})}, + volume = {39}, + number = {5}, + pages = {S461--S476}, + author = {Nguyenho Ho and Sarah D. Olson and Homer F. Walker}, + title = {Accelerating the Uzawa Algorithm}, + journal = {{SIAM} Journal on Scientific Computing} +} diff --git a/dumux/linear/istlsolverfactorybackend.hh b/dumux/linear/istlsolverfactorybackend.hh index 776ba71309b41d0197704719ef06a8a7afbf79ef..2de5db96946ab333b7c6c9b383a4c1d452251bdb 100644 --- a/dumux/linear/istlsolverfactorybackend.hh +++ b/dumux/linear/istlsolverfactorybackend.hh @@ -37,16 +37,63 @@ // preconditioners #include #include +#include "preconditioners.hh" // solvers #include #include +#include #include #include namespace Dumux { +// initSolverFactoriesForMultiTypeBlockMatrix differs in different compilation units, +// so we have it in an anonymous namespace +namespace { + +/*! + * \brief Initializes the direct solvers, preconditioners and iterative solvers in + * the factories with the corresponding Matrix and Vector types. + * \note We currently consider only direct solvers and preconditioners provided by + * Dumux which, unlike the ones implemented in Dune, also support MultiTypeBlockMatrices. + * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices. + * \tparam LinearOperator the linear operator type + */ +template +int initSolverFactoriesForMultiTypeBlockMatrix() +{ + using M = typename LinearOperator::matrix_type; + using X = typename LinearOperator::range_type; + using Y = typename LinearOperator::domain_type; + using TL = Dune::TypeList; + auto& dsfac = Dune::DirectSolverFactory::instance(); + Dune::addRegistryToFactory(dsfac, Dumux::MultiTypeBlockMatrixDirectSolverTag{}); + auto& pfac = Dune::PreconditionerFactory::instance(); + Dune::addRegistryToFactory(pfac, Dumux::MultiTypeBlockMatrixPreconditionerTag{}); + using TLS = Dune::TypeList; + auto& isfac = Dune::IterativeSolverFactory::instance(); + return Dune::addRegistryToFactory(isfac, Dune::IterativeSolverTag{}); +} +} // end namespace + +/*! + * \brief Initialize the solver factories for regular matrices or MultiTypeBlockMatrices + * \tparam Matrix the matrix + * \tparam LinearOperator the linear operator + * + * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices. + */ +template +void initSolverFactories() +{ + if constexpr (isMultiTypeBlockMatrix::value) + initSolverFactoriesForMultiTypeBlockMatrix(); + else + Dune::initSolverFactories(); +} + /*! * \ingroup Linear * \brief A linear solver using the dune-istl solver factory @@ -148,7 +195,10 @@ private: template void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b) { - if constexpr (LinearSolverTraits::canCommunicate) + // 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::value) { if (isParallel_) { @@ -181,7 +231,7 @@ private: if (firstCall_) { - Dune::initSolverFactories(); + initSolverFactories(); parallelHelper_->initGhostsAndOwners(); } @@ -207,7 +257,7 @@ private: auto linearOperator = std::make_shared(A); if (firstCall_) - Dune::initSolverFactories(); + initSolverFactories(); // construct solver auto solver = getSolverFromFactory_(linearOperator); diff --git a/dumux/linear/istlsolverregistry.hh b/dumux/linear/istlsolverregistry.hh new file mode 100644 index 0000000000000000000000000000000000000000..6caa4c4b2588b3b87b25a526e2360df8bc6c26b0 --- /dev/null +++ b/dumux/linear/istlsolverregistry.hh @@ -0,0 +1,56 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Linear + * \brief The specialized Dumux macro and tag for the ISTL registry to choose the + * solver and preconditioner at runtime + */ +#ifndef DUMUX_LINEAR_ISTL_SOLVER_REGISTRY_HH +#define DUMUX_LINEAR_ISTL_SOLVER_REGISTRY_HH + +#include +#if DUNE_VERSION_NEWER_REV(DUNE_ISTL,2,7,1) + +#include +#include + +/*! + * \brief Register a Dumux preconditioner + * + * Use this macro in namespace Dumux. + * Example: + * DUMUX_REGISTER_PRECONDITIONER("mypreconditioner", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator()); + * Expicitly specifying the namespaces is required. + * Set parameter Preconditioner.Type to "mypreconditioner" to use it through the factory. + */ +#define DUMUX_REGISTER_PRECONDITIONER(name, tag, ...) \ +} namespace Dune { \ +DUNE_REGISTRY_PUT(tag, name, __VA_ARGS__); \ +} namespace Dumux { \ + +namespace Dumux { +namespace { +struct MultiTypeBlockMatrixPreconditionerTag {}; +struct MultiTypeBlockMatrixDirectSolverTag {}; +} // end namespace +} // end namespace Dumux + +#endif // DUNE version check +#endif // DUMUX_LINEAR_ISTL_SOLVER_REGISTRY_HH diff --git a/dumux/linear/linearsolvertraits.hh b/dumux/linear/linearsolvertraits.hh index 154071021c4551d148cc1c349c73ee30a5ee5f8d..b61ba21e14b26e44f58906c62bac35afaf6afd00 100644 --- a/dumux/linear/linearsolvertraits.hh +++ b/dumux/linear/linearsolvertraits.hh @@ -177,6 +177,11 @@ template struct LinearSolverTraitsImpl : public LinearSolverTraitsImpl {}; +//! staggered: use overlapping model TODO provide staggered-specific traits, combining overlapping/non-overlapping +template +struct LinearSolverTraitsImpl +: public LinearSolverTraitsImpl {}; + } // end namespace Dumux #endif // DUMUX_LINEAR_SOLVER_TRAITS_HH diff --git a/dumux/linear/preconditioners.hh b/dumux/linear/preconditioners.hh new file mode 100644 index 0000000000000000000000000000000000000000..16df85cfbadae070aa70ab252116b34516e052b5 --- /dev/null +++ b/dumux/linear/preconditioners.hh @@ -0,0 +1,340 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup Linear + * \brief Dumux preconditioners for iterative solvers + */ +#ifndef DUMUX_LINEAR_PRECONDITIONERS_HH +#define DUMUX_LINEAR_PRECONDITIONERS_HH + +#include +#include +#include +#include +#include + +#if HAVE_UMFPACK +#include +#endif + +#include +#include + +#if DUNE_VERSION_NEWER_REV(DUNE_ISTL,2,7,1) +#include +#endif + +namespace Dumux { + +#if DUNE_VERSION_GTE(DUNE_ISTL,2,7) +/*! + * \ingroup Linear + * \brief A preconditioner based on the Uzawa algorithm for saddle-point problems of the form + * \f$ + \begin{pmatrix} + A & B \\ + C & D + \end{pmatrix} + + \begin{pmatrix} + u\\ + p + \end{pmatrix} + + = + + \begin{pmatrix} + f\\ + g + \end{pmatrix} + * \f$ + * + * This preconditioner is especially suited for solving the incompressible (Navier-)Stokes equations. + * Here, \f$D = 0\f$ and \f$B = C^T\f$ if \f$\rho = 1\f$. + * We do not expect good convergence if energy or mass transport is considered. + * + * See: Benzi, M., Golub, G. H., & Liesen, J. (2005). Numerical solution of saddle point problems. Acta numerica, 14, 1-137 \cite benzi2005 and
+ * Ho, N., Olson, S. D., & Walker, H. F. (2017). Accelerating the Uzawa algorithm. SIAM Journal on Scientific Computing, 39(5), S461-S476 \cite ho2017 + * + * \tparam M Type of the matrix. + * \tparam X Type of the update. + * \tparam Y Type of the defect. + * \tparam l Preconditioner block level (for compatibility reasons, unused). + */ +template +class SeqUzawa : public Dune::Preconditioner +{ + static_assert(Dumux::isMultiTypeBlockMatrix::value && M::M() == 2 && M::N() == 2, "SeqUzawa expects a 2x2 MultiTypeBlockMatrix."); + static_assert(l== 1, "SeqUzawa 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; + using UMFPackSolverForA = Dune::UMFPack; + +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. + */ + SeqUzawa(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)) + { + const bool determineRelaxationFactor = getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.DetermineRelaxationFactor", true); + + // AMG is needed for determination of omega + if (determineRelaxationFactor || !useDirectVelocitySolverForA_) + initAMG_(params); + + if (useDirectVelocitySolverForA_) + initUMFPack_(); + + if (determineRelaxationFactor) + relaxationFactor_ = estimateOmega_(); + } + + /*! + * \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()]; + } + + // the actual Uzawa iteration + for (std::size_t k = 0; k < numIterations_; ++k) + { + // u_k+1 = u_k + Q_A^−1*(f − (A*u_k + B*p_k)), + auto uRhs = f; + A.mmv(u, uRhs); + B.mmv(p, uRhs); + auto uIncrement = u; + applySolverForA_(uIncrement, uRhs); + u += uIncrement; + + // p_k+1 = p_k + omega*(g - C*u_k+1 - D*p_k) + auto pIncrement = g; + C.mmv(u, pIncrement); + D.mmv(p, pIncrement); + pIncrement *= relaxationFactor_; + p += pIncrement; + + if (verbosity_ > 1) + { + std::cout << "Uzawa iteration " << k + << ", residual: " << uRhs.two_norm() + pIncrement.two_norm()/relaxationFactor_ << std::endl; + } + } + } + + /*! + * \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 + } + + /*! + * \brief Estimate the relaxation factor omega + * + * The optimal relaxation factor is omega = 2/(lambdaMin + lambdaMax), where lambdaMin and lambdaMax are the smallest and largest + * eigenvalues of the Schur complement -C*Ainv*B (assuming D = 0). + * lambdaMax can be easily determined using the power iteration algorithm (https://en.wikipedia.org/wiki/Power_iteration) and lambdaMin + * could be estimated in a similar manner. We do not consider lambdaMin because for certain cases, e.g., when C contains some rows of zeroes only, + * this estimate will fail. + * + * Instead we assume that lambdaMin is sufficiently close to lambdaMax such that omega = 1/lambdaMax. + * This seems to work rather well for various applications. + * We will underestimate omega by a factor of 2 in the worst case (i.e, lambdaMin = 0). + * + * When facing convergence issues, you may set LinearSolver.Preconditioner.Verbosity = 1 to see the estimate of lambdaMax. + * In a new simulation run, you can then set LinearSolver.Preconditioner.DetermineRelaxationFactor = false and set some other value + * for LinearSolver.Preconditioner.Relaxation based on the estimate of lambdaMax. + * + * See: Benzi, M., Golub, G. H., & Liesen, J. (2005). Numerical solution of saddle point problems. Acta numerica, 14, 1-137. + */ + scalar_field_type estimateOmega_() const + { + using namespace Dune::Indices; + auto& A = matrix_[_0][_0]; + auto& B = matrix_[_0][_1]; + auto& C = matrix_[_1][_0]; + + U x(A.M()); + x = 1.0; + + scalar_field_type omega = 0.0; + scalar_field_type lambdaMax = 0.0; + + static const auto iterations = Dumux::getParamFromGroup(paramGroup_, "LinearSolver.Preconditioner.PowerLawIterations", 5); + + // apply power iteration x_k+1 = M*x_k/|M*x_k| for the matrix M = -C*Ainv*B + for (std::size_t i = 0; i < iterations; ++i) + { + // bx = B*x + U bx(x.size()); + B.mv(x, bx); + + // ainvbx = Ainv*(B*x) + auto ainvbx = x; + applySolverForA_(ainvbx, bx); + + // v = M*x = -C*(Ainv*B*x) + U v(x.size()); + C.mv(ainvbx, v); + v *= -1.0; + + // eigenvalue lambdaMax = xt*M*x/(xt*x) = xt*v/(xt*x); + lambdaMax = x.dot(v)/(x.dot(x)); + + // relaxation factor omega = 1/lambda; + omega = 1.0/lambdaMax; + + // new iterate x = M*x/|M*x| = v/|v| + x = v; + x /= v.two_norm(); + } + + if (verbosity_ > 0) + { + std::cout << "\n*** Uzawa Preconditioner ***" << std::endl; + std::cout << "Estimating relaxation factor based on Schur complement" << std::endl; + std::cout << "Largest estimated eigenvalue lambdaMax = " << lambdaMax << std::endl; + std::cout << "Relaxation factor omega = 1/lambdaMax = " << omega << std::endl; + } + + return omega; + } + + 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_; +}; +#endif // DUNE_VERSION_GTE(DUNE_ISTL,2,7) + +#if DUNE_VERSION_NEWER_REV(DUNE_ISTL,2,7,1) +DUMUX_REGISTER_PRECONDITIONER("uzawa", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator()); +#endif + +} // end namespace Dumux + +#endif // DUMUX_LINEAR_PRECONDITIONERS_HH diff --git a/dumux/linear/seqsolverbackend.hh b/dumux/linear/seqsolverbackend.hh index d06a0a737f23b913f3e70590acd16bb71ca4aaf0..00db276e88d6979ddfdf50a1c7f02936e2926a3f 100644 --- a/dumux/linear/seqsolverbackend.hh +++ b/dumux/linear/seqsolverbackend.hh @@ -32,6 +32,8 @@ #include #include #include +#include +#include #include #include @@ -40,6 +42,8 @@ #include #include #include +#include +#include namespace Dumux { @@ -147,6 +151,29 @@ public: return result.converged; } + +#if DUNE_VERSION_GTE(DUNE_ISTL,2,7) + // solve with generic parameter tree + template + static bool solveWithParamTree(const Matrix& A, Vector& x, const Vector& b, + const Dune::ParameterTree& params) + { + + // make a linear operator from a matrix + using MatrixAdapter = Dune::MatrixAdapter; + const auto linearOperator = std::make_shared(A); + + auto precond = std::make_shared(linearOperator, params.sub("preconditioner")); + Solver solver(linearOperator, precond, params); + + Vector bTmp(b); + + Dune::InverseOperatorResult result; + solver.apply(x, bTmp, result); + + return result.converged; + } +#endif }; /*! @@ -876,6 +903,33 @@ private: */ // \{ +#if DUNE_VERSION_GTE(DUNE_ISTL,2,7) +/*! + * \ingroup Linear + * \brief A Uzawa preconditioned BiCGSTAB solver for saddle-point problems + */ +template +class UzawaBiCGSTABBackend : public LinearSolver +{ +public: + using LinearSolver::LinearSolver; + + template + bool solve(const Matrix& A, Vector& x, const Vector& b) + { + using Preconditioner = SeqUzawa; + using Solver = Dune::BiCGSTABSolver; + static const auto solverParams = LinearSolverParameters::createParameterTree(this->paramGroup()); + return IterativePreconditionedSolverImpl::template solveWithParamTree(A, x, b, solverParams); + } + + std::string name() const + { + return "Uzawa preconditioned BiCGSTAB solver"; + } +}; +#endif + /*! * \ingroup Linear * \brief A simple ilu0 block diagonal preconditioner diff --git a/test/freeflow/navierstokes/kovasznay/main.cc b/test/freeflow/navierstokes/kovasznay/main.cc index c21d9a0936fbfe11e162a5ee14395a4ca4a09df6..77175fa6c4b7b32e89d4b540194cc2e12beeb034 100644 --- a/test/freeflow/navierstokes/kovasznay/main.cc +++ b/test/freeflow/navierstokes/kovasznay/main.cc @@ -32,6 +32,7 @@ #include #include #include +#include #include #include @@ -39,11 +40,15 @@ #include #include #include -#include +#include #include #include #include +#if DUNE_VERSION_GTE(DUNE_ISTL,2,7) +#include +#endif + #include "problem.hh" int main(int argc, char** argv) try @@ -108,7 +113,11 @@ int main(int argc, char** argv) try auto assembler = std::make_shared(problem, gridGeometry, gridVariables); // the linear solver +#if DUNE_VERSION_GTE(DUNE_ISTL,2,7) + using LinearSolver = Dumux::UzawaBiCGSTABBackend>; +#else using LinearSolver = Dumux::UMFPackBackend; +#endif auto linearSolver = std::make_shared(); // the non-linear solver diff --git a/test/freeflow/navierstokes/sincos/CMakeLists.txt b/test/freeflow/navierstokes/sincos/CMakeLists.txt index b11ea227015f096e71078e60e353dfd0e9282422..2be99bb669d38aa9bb91975c06a3ace4fcc1a14f 100644 --- a/test/freeflow/navierstokes/sincos/CMakeLists.txt +++ b/test/freeflow/navierstokes/sincos/CMakeLists.txt @@ -6,6 +6,7 @@ dumux_add_test(NAME test_ff_navierstokes_sincos LABELS freeflow TIMEOUT 1000 CMAKE_GUARD HAVE_UMFPACK + COMPILE_DEFINITIONS LINEARSOLVER=UMFPackBackend COMMAND ./convergencetest.py CMD_ARGS test_ff_navierstokes_sincos params.input -Grid.UpperRight "6.28 6.28" @@ -18,6 +19,7 @@ dumux_add_test(NAME test_ff_navierstokes_sincos_instationary TARGET test_ff_navierstokes_sincos LABELS freeflow CMAKE_GUARD HAVE_UMFPACK + COMPILE_DEFINITIONS LINEARSOLVER=UMFPackBackend COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py CMD_ARGS --script fuzzy --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_instationary-reference.vtu @@ -28,3 +30,20 @@ dumux_add_test(NAME test_ff_navierstokes_sincos_instationary -Problem.Name test_ff_navierstokes_sincos_instationary -Problem.IsStationary false -Component.LiquidKinematicViscosity 0.1") + +dumux_add_test(NAME test_ff_navierstokes_sincos_uzawapreconditioner_factory + SOURCES main.cc + LABELS freeflow + TIMEOUT 5000 + CMAKE_GUARD "( HAVE_UMFPACK AND DUNE_ISTL_VERSION GREATER_EQUAL 2.7 )" + COMPILE_DEFINITIONS LINEARSOLVER=IstlSolverFactoryBackend> + COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py + CMD_ARGS --script fuzzy + --files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_sincos_instationary-reference.vtu + ${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_uzawapreconditioner-00017.vtu + --command "${CMAKE_CURRENT_BINARY_DIR}/test_ff_navierstokes_sincos_uzawapreconditioner_factory params.input + -Grid.UpperRight '1 1' + -Grid.Cells '50 50' + -Problem.Name test_ff_navierstokes_sincos_uzawapreconditioner + -Problem.IsStationary false + -Component.LiquidKinematicViscosity 0.1") diff --git a/test/freeflow/navierstokes/sincos/main.cc b/test/freeflow/navierstokes/sincos/main.cc index 507ff765ca274320d2dcdab3e33761ca222aa6aa..12f0bc9facbf1da6b6931c34b5be5ae7f1202acf 100644 --- a/test/freeflow/navierstokes/sincos/main.cc +++ b/test/freeflow/navierstokes/sincos/main.cc @@ -41,11 +41,16 @@ #include #include #include -#include +#include #include #include #include +#include +#if DUNE_VERSION_NEWER_REV(DUNE_ISTL,2,7,1) +#include +#endif + #include "problem.hh" /*! @@ -252,7 +257,7 @@ int main(int argc, char** argv) try : std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld); // the linear solver - using LinearSolver = Dumux::UMFPackBackend; + using LinearSolver = LINEARSOLVER; auto linearSolver = std::make_shared(); // the non-linear solver diff --git a/test/freeflow/navierstokes/sincos/params.input b/test/freeflow/navierstokes/sincos/params.input index 80f5d2becadc5c3abb0bff3c8bf8cf310fb9afd0..0da3caf350b99c444f236021777946ddeb951f23 100644 --- a/test/freeflow/navierstokes/sincos/params.input +++ b/test/freeflow/navierstokes/sincos/params.input @@ -4,7 +4,8 @@ TEnd = 1.0 # [s] [Grid] LowerLeft = 0 0 -UpperRight = 1 1 +UpperRight = 6.28 6.28 +Cells = 20 20 [Problem] Name = test_ff_sincos @@ -34,3 +35,13 @@ AddProcessRank = false [Flux] UpwindWeight = 0.5 + +[LinearSolver] +Type = restartedflexiblegmressolver +Verbosity = 1 +GMResRestart = 50 + +[LinearSolver.Preconditioner] +Type = uzawa +Verbosity = 1 +Iterations = 5