diff --git a/ewoms-prereqs.cmake b/ewoms-prereqs.cmake index 6443f99984..21dfe66952 100644 --- a/ewoms-prereqs.cmake +++ b/ewoms-prereqs.cmake @@ -13,8 +13,12 @@ set (ewoms_CONFIG_VAR HAVE_DUNE_ISTL HAVE_DUNE_ALUGRID HAVE_DUNE_FEM + HAVE_PETSC + HAVE_AMGXSOLVER HAVE_ECL_INPUT HAVE_ECL_OUTPUT + HAVE_VIENNACL + HAVE_OPENCL DUNE_AVOID_CAPABILITIES_IS_PARALLEL_DEPRECATION_WARNING ) @@ -36,6 +40,12 @@ set (ewoms_DEPS "dune-alugrid" "dune-fem" "opm-grid" + # PETSc numerical backend + "PETSc" + # AMGX wrapper using PETSc + "AmgXSolver" + # ViennaCL + "ViennaCL" # valgrind client requests "Valgrind" # quadruple precision floating point calculations diff --git a/ewoms/common/basicproperties.hh b/ewoms/common/basicproperties.hh index 680fc91ab1..2195f13635 100644 --- a/ewoms/common/basicproperties.hh +++ b/ewoms/common/basicproperties.hh @@ -78,6 +78,10 @@ NEW_PROP_TAG(GridView); #if HAVE_DUNE_FEM NEW_PROP_TAG(GridPart); +//! The class describing the discrete function space when dune-fem is used, otherwise it points to the stencil class +NEW_PROP_TAG(DiscreteFunctionSpace); +NEW_PROP_TAG(DiscreteFunction); +NEW_PROP_TAG(LinearOperator); #endif //! Property which tells the Vanguard how often the grid should be refined diff --git a/ewoms/disc/common/fvbasediscretization.hh b/ewoms/disc/common/fvbasediscretization.hh index 0d72f0241e..9f8221d3e3 100644 --- a/ewoms/disc/common/fvbasediscretization.hh +++ b/ewoms/disc/common/fvbasediscretization.hh @@ -76,7 +76,13 @@ #include #include #include + +#if HAVE_PETSC +#include #endif +#include +#include +#endif // endif HAVE_DUNE_FEM #include #include @@ -130,7 +136,66 @@ SET_TYPE_PROP(FvBaseDiscretization, DiscExtensiveQuantities, Ewoms::FvBaseExtens //! Calculates the gradient of any quantity given the index of a flux approximation point SET_TYPE_PROP(FvBaseDiscretization, GradientCalculator, Ewoms::FvBaseGradientCalculator); -//! Set the type of a global jacobian matrix from the solution types +//! Set the type of a global Jacobian matrix from the solution types +#if HAVE_DUNE_FEM +SET_PROP(FvBaseDiscretization, DiscreteFunction) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; +public: + // discrete function storing solution data + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction type; +}; +#endif + +#if USE_DUNE_FEM_SOLVERS +SET_PROP(FvBaseDiscretization, SparseMatrixAdapter) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + // discrete function storing solution data + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction DiscreteFunction; + +#if USE_DUNE_FEM_PETSC_SOLVERS +#warning "Using Dune-Fem PETSc solvers" + typedef Dune::Fem::PetscLinearOperator< DiscreteFunction, DiscreteFunction > LinearOperator; +#elif USE_DUNE_FEM_VIENNACL_SOLVERS +#warning "Using Dune-Fem ViennaCL solvers" + typedef Dune::Fem::SparseRowLinearOperator < DiscreteFunction, DiscreteFunction > LinearOperator; +#else +#warning "Using Dune-Fem ISTL solvers" + typedef Dune::Fem::ISTLLinearOperator < DiscreteFunction, DiscreteFunction > LinearOperator; +#endif + + struct FemMatrixBackend : public LinearOperator + { + typedef LinearOperator ParentType; + typedef typename LinearOperator :: MatrixType Matrix; + typedef typename ParentType :: MatrixBlockType MatrixBlock; + template + FemMatrixBackend( const Simulator& simulator ) + : LinearOperator("eWoms::Jacobian", simulator.model().space(), simulator.model().space() ) + {} + + void commit() + { + this->flushAssembly(); + } + + template< class LocalBlock > + void addToBlock ( const size_t row, const size_t col, const LocalBlock& block ) + { + this->addBlock( row, col, block ); + } + + void clearRow( const size_t row, const Scalar diag = 1.0 ) { this->unitRow( row ); } + }; +public: + typedef FemMatrixBackend type; +}; +#else SET_PROP(FvBaseDiscretization, SparseMatrixAdapter) { private: @@ -141,6 +206,7 @@ private: public: typedef typename Ewoms::Linear::IstlSparseMatrixAdapter type; }; +#endif //! The maximum allowed number of timestep divisions for the //! Newton solver @@ -342,13 +408,15 @@ class FvBaseDiscretization typedef typename LocalResidual::LocalEvalBlockVector LocalEvalBlockVector; + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + class BlockVectorWrapper { protected: SolutionVector blockVector_; public: - BlockVectorWrapper(const std::string& name OPM_UNUSED, const size_t size) - : blockVector_(size) + BlockVectorWrapper(const std::string& name OPM_UNUSED, const DiscreteFunctionSpace& space) + : blockVector_(space.size()) {} SolutionVector& blockVector() @@ -358,14 +426,12 @@ class FvBaseDiscretization }; #if HAVE_DUNE_FEM - typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; - // discrete function storing solution data - typedef Dune::Fem::ISTLBlockVectorDiscreteFunction DiscreteFunction; + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; // problem restriction and prolongation operator for adaptation - typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; - typedef typename Problem :: RestrictProlongOperator ProblemRestrictProlongOperator; + typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem; + typedef typename Problem :: RestrictProlongOperator ProblemRestrictProlongOperator; // discrete function restriction and prolongation operator for adaptation typedef Dune::Fem::RestrictProlongDefault< DiscreteFunction > DiscreteFunctionRestrictProlong; @@ -374,7 +440,6 @@ class FvBaseDiscretization typedef Dune::Fem::AdaptationManager AdaptationManager; #else typedef BlockVectorWrapper DiscreteFunction; - typedef size_t DiscreteFunctionSpace; #endif // copying a discretization object is not a good idea @@ -394,14 +459,14 @@ public: , elementMapper_(gridView_) , vertexMapper_(gridView_) #endif - , newtonMethod_(simulator) - , localLinearizer_(ThreadManager::maxThreads()) - , linearizer_(new Linearizer()) #if HAVE_DUNE_FEM - , space_( simulator.vanguard().gridPart() ) + , discreteFunctionSpace_( simulator.vanguard().gridPart() ) #else - , space_( asImp_().numGridDof() ) + , discreteFunctionSpace_( asImp_().numGridDof() ) #endif + , newtonMethod_(simulator) + , localLinearizer_(ThreadManager::maxThreads()) + , linearizer_(new Linearizer( )) , enableGridAdaptation_( EWOMS_GET_PARAM(TypeTag, bool, EnableGridAdaptation) ) , enableIntensiveQuantityCache_(EWOMS_GET_PARAM(TypeTag, bool, EnableIntensiveQuantityCache)) , enableStorageCache_(EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache)) @@ -426,7 +491,7 @@ public: size_t numDof = asImp_().numGridDof(); for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) { - solution_[timeIdx].reset(new DiscreteFunction("solution", space_)); + solution_[timeIdx].reset(new DiscreteFunction("solution", discreteFunctionSpace_)); if (storeIntensiveQuantities()) { intensiveQuantityCache_[timeIdx].resize(numDof); @@ -1091,6 +1156,16 @@ public: SolutionVector& solution(unsigned timeIdx) { return solution_[timeIdx]->blockVector(); } + template + void communicate( BVector& x ) const + { +#if HAVE_DUNE_FEM + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction DF; + DF tmpX("temp-x", discreteFunctionSpace_, x); + tmpX.communicate(); +#endif + } + protected: /*! * \copydoc solution(int) const @@ -1508,7 +1583,7 @@ public: void resetLinearizer () { delete linearizer_; - linearizer_ = new Linearizer; + linearizer_ = new Linearizer(); linearizer_->init(simulator_); } @@ -1742,6 +1817,11 @@ public: solution(timeIdx).resize(numDof); auxMod->applyInitial(); + +#if DUNE_VERSION_NEWER( DUNE_FEM, 2, 7 ) + discreteFunctionSpace_.extendSize( asImp_().numAuxiliaryDof() ); +#endif + } /*! @@ -1808,6 +1888,11 @@ public: const Ewoms::Timer& updateTimer() const { return updateTimer_; } +#if HAVE_DUNE_FEM + const DiscreteFunctionSpace& space() const { return discreteFunctionSpace_; } +#endif + + protected: void resizeAndResetIntensiveQuantitiesCache_() { @@ -1878,9 +1963,18 @@ protected: ElementMapper elementMapper_; VertexMapper vertexMapper_; + DiscreteFunctionSpace discreteFunctionSpace_; + // a vector with all auxiliary equations to be considered std::vector*> auxEqModules_; + mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_; + +#if HAVE_DUNE_FEM + std::unique_ptr< RestrictProlong > restrictProlong_; + std::unique_ptr< AdaptationManager> adaptationManager_; +#endif + NewtonMethod newtonMethod_; Ewoms::Timer prePostProcessTimer_; @@ -1899,15 +1993,6 @@ protected: mutable IntensiveQuantitiesVector intensiveQuantityCache_[historySize]; mutable std::vector intensiveQuantityCacheUpToDate_[historySize]; - DiscreteFunctionSpace space_; - mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_; - -#if HAVE_DUNE_FEM - std::unique_ptr restrictProlong_; - std::unique_ptr adaptationManager_; -#endif - - std::list*> outputModules_; Scalar gridTotalVolume_; diff --git a/ewoms/disc/common/fvbaselinearizer.hh b/ewoms/disc/common/fvbaselinearizer.hh index 0bfe1eacf1..6a055e6842 100644 --- a/ewoms/disc/common/fvbaselinearizer.hh +++ b/ewoms/disc/common/fvbaselinearizer.hh @@ -54,6 +54,7 @@ namespace Ewoms { template class EcfvDiscretization; + /*! * \ingroup FiniteVolumeDiscretizations * @@ -95,8 +96,6 @@ class FvBaseLinearizer typedef GlobalEqVector Vector; - typedef typename SparseMatrixAdapter::IstlMatrix IstlMatrix; - enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; enum { historySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) }; @@ -229,9 +228,6 @@ public: */ void linearizeAuxiliaryEquations() { - // flush possible local caches into matrix structure - jacobian_->commit(); - auto& model = model_(); const auto& comm = simulator_().gridView().comm(); for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) { @@ -262,6 +258,9 @@ public: if (!succeeded) throw Opm::NumericalIssue("linearization of an auxilary equation failed"); } + + // flush possible local caches into matrix structure + jacobian_->commit(); } /*! @@ -499,6 +498,9 @@ private: std::rethrow_exception(exceptionPtr); } + // flush possible local caches into matrix structure + jacobian_->commit(); + applyConstraintsToLinearization_(); } diff --git a/ewoms/disc/common/fvbaseproperties.hh b/ewoms/disc/common/fvbaseproperties.hh index e1bf3a6570..d372b7483b 100644 --- a/ewoms/disc/common/fvbaseproperties.hh +++ b/ewoms/disc/common/fvbaseproperties.hh @@ -37,6 +37,9 @@ #include #include #include +#include +#include +#include BEGIN_PROPERTIES @@ -54,10 +57,19 @@ NEW_PROP_TAG(ParallelBiCGStabLinearSolver); NEW_PROP_TAG(LocalLinearizerSplice); NEW_PROP_TAG(FiniteDifferenceLocalLinearizer); +NEW_PROP_TAG(DiscreteFunctionSpace); + SET_SPLICES(FvBaseDiscretization, LinearSolverSplice, LocalLinearizerSplice); //! use a parallel BiCGStab linear solver by default +#if USE_AMGX_SOLVERS +SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, AmgXSolverBackend); +#elif USE_DUNE_FEM_SOLVERS +SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, FemSolverBackend); +#else SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, ParallelBiCGStabLinearSolver); +//SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, ParallelIstlLinearSolver); +#endif //! by default, use finite differences to linearize the system of PDEs SET_TAG_PROP(FvBaseDiscretization, LocalLinearizerSplice, FiniteDifferenceLocalLinearizer); @@ -80,9 +92,6 @@ NEW_PROP_TAG(GridView); //! The class describing the stencil of the spatial discretization NEW_PROP_TAG(Stencil); -//! The class describing the discrete function space when dune-fem is used, otherwise it points to the stencil class -NEW_PROP_TAG(DiscreteFunctionSpace); - //! The type of the problem NEW_PROP_TAG(Problem); //! The type of the base class for all problems which use this model diff --git a/ewoms/disc/ecfv/ecfvdiscretization.hh b/ewoms/disc/ecfv/ecfvdiscretization.hh index 9ae5805e5b..b2a789c769 100644 --- a/ewoms/disc/ecfv/ecfvdiscretization.hh +++ b/ewoms/disc/ecfv/ecfvdiscretization.hh @@ -91,6 +91,24 @@ private: public: typedef Dune::Fem::FiniteVolumeSpace< FunctionSpace, GridPart, 0 > type; }; +#else +SET_PROP(EcfvDiscretization, DiscreteFunctionSpace) +{ +private: + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + struct DiscreteFunctionSpace + { + static const int dimRange = numEq ; + size_t numGridDofs_; + size_t extension_; + DiscreteFunctionSpace( const size_t numGridDofs ) + : numGridDofs_( numGridDofs ), extension_(0) {} + void extendSize( const size_t extension ) { extension_ = extension; } + size_t size() const { return dimRange * (numGridDofs_ + extension_); } + }; +public: + typedef DiscreteFunctionSpace type; +}; #endif //! Set the border list creator for to the one of an element based diff --git a/ewoms/disc/vcfv/vcfvdiscretization.hh b/ewoms/disc/vcfv/vcfvdiscretization.hh index 47d25fb3a2..042c194ce6 100644 --- a/ewoms/disc/vcfv/vcfvdiscretization.hh +++ b/ewoms/disc/vcfv/vcfvdiscretization.hh @@ -101,6 +101,24 @@ public: // Lagrange discrete function space with unknowns at the cell vertices typedef Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpace, GridPart, 1 > type; }; +#else +SET_PROP(VcfvDiscretization, DiscreteFunctionSpace) +{ +private: + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + struct DiscreteFunctionSpace + { + static const int dimRange = numEq ; + size_t numGridDofs_; + size_t extension_; + DiscreteFunctionSpace( const size_t numGridDofs ) + : numGridDofs_( numGridDofs ), extension_(0) {} + void extendSize( const size_t extension ) { extension_ = extension; } + size_t size() const { return dimRange * (numGridDofs_ + extension_); } + }; +public: + typedef DiscreteFunctionSpace type; +}; #endif //! Set the border list creator for vertices diff --git a/ewoms/linear/amgxsolverbackend.hh b/ewoms/linear/amgxsolverbackend.hh new file mode 100644 index 0000000000..40af005968 --- /dev/null +++ b/ewoms/linear/amgxsolverbackend.hh @@ -0,0 +1,364 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM 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 2 of the License, or + (at your option) any later version. + + OPM 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 OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Ewoms::Linear::AmgXSolverBackend + */ +#ifndef EWOMS_AMGX_SOLVER_BACKEND_HH +#define EWOMS_AMGX_SOLVER_BACKEND_HH + +#include + +#if USE_AMGX_SOLVERS + +#if ! HAVE_PETSC || ! HAVE_AMGXSOLVER +#error "PETSc and AmgXSolver is needed for the AMGX solver backend" +#endif + +#define DISABLE_AMG_DIRECTSOLVER 1 +#include +#include +#include +#include + +#include + +#include + + +#include +#include + +#include + +#include +#include +#include + +namespace Ewoms { +namespace Linear { +template +class AmgXSolverBackend; +}} // namespace Linear, Ewoms + + +BEGIN_PROPERTIES + +NEW_TYPE_TAG(AmgXSolverBackend); + +SET_TYPE_PROP(AmgXSolverBackend, + LinearSolverBackend, + Ewoms::Linear::AmgXSolverBackend); + +NEW_PROP_TAG(LinearSolverTolerance); +NEW_PROP_TAG(LinearSolverAbsTolerance); +NEW_PROP_TAG(LinearSolverMaxIterations); +NEW_PROP_TAG(LinearSolverVerbosity); +NEW_PROP_TAG(LinearSolverMaxError); +NEW_PROP_TAG(LinearSolverOverlapSize); +//! The order of the sequential preconditioner +NEW_PROP_TAG(PreconditionerOrder); + +//! The relaxation factor of the preconditioner +NEW_PROP_TAG(PreconditionerRelaxation); + +//! Solver mode for AMGX solver +NEW_PROP_TAG(AmgxSolverMode); + +//! Filename for AMGX solver configuration +NEW_PROP_TAG(AmgxSolverConfigFileName); + +//! Filename for DUNE-FEM solver parameters +NEW_PROP_TAG(FemSolverParameterFileName); + +//! make the linear solver shut up by default +SET_INT_PROP(AmgXSolverBackend, LinearSolverVerbosity, 0); + +//! set the default number of maximum iterations for the linear solver +SET_INT_PROP(AmgXSolverBackend, LinearSolverMaxIterations, 1000); + +SET_SCALAR_PROP(AmgXSolverBackend, LinearSolverMaxError, 1e7); + +//! set the default overlap size to 2 +SET_INT_PROP(AmgXSolverBackend, LinearSolverOverlapSize, 2); + +//! set the preconditioner order to 0 by default +SET_INT_PROP(AmgXSolverBackend, PreconditionerOrder, 0); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_SCALAR_PROP(AmgXSolverBackend, PreconditionerRelaxation, 1.0); + +//! make the linear solver shut up by default +SET_SCALAR_PROP(AmgXSolverBackend, LinearSolverTolerance, 0.01); + +//! make the linear solver shut up by default +SET_SCALAR_PROP(AmgXSolverBackend, LinearSolverAbsTolerance, 0.01); + +//! set default solver mode for AMGX solver configuratrion +SET_STRING_PROP(AmgXSolverBackend, AmgxSolverMode, "dDDI"); + +//! set default filename for AMGX solver configuratrion +SET_STRING_PROP(AmgXSolverBackend, AmgxSolverConfigFileName, "amgxconfig.json"); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_STRING_PROP(AmgXSolverBackend, FemSolverParameterFileName, ""); + +END_PROPERTIES + +namespace Ewoms { +namespace Linear { +/*! + * \ingroup Linear + * + * \brief Provides the common code which is required by most linear solvers. + * + * This class provides access to all preconditioners offered by dune-istl using the + * PreconditionerWrapper property: + * \code + * SET_TYPE_PROP(YourTypeTag, PreconditionerWrapper, + * Ewoms::Linear::PreconditionerWrapper$PRECONDITIONER); + * \endcode + * + * Where the choices possible for '\c $PRECONDITIONER' are: + * - \c Jacobi: A Jacobi preconditioner + * - \c GaussSeidel: A Gauss-Seidel preconditioner + * - \c SSOR: A symmetric successive overrelaxation (SSOR) preconditioner + * - \c SOR: A successive overrelaxation (SOR) preconditioner + * - \c ILUn: An ILU(n) preconditioner + * - \c ILU0: An ILU(0) preconditioner. The results of this + * preconditioner are the same as setting the + * PreconditionerOrder property to 0 and using the ILU(n) + * preconditioner. The reason for the existence of ILU0 is + * that it is computationally cheaper because it does not + * need to consider things which are only required for + * higher orders + */ +template +class AmgXSolverBackend +{ +protected: + typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) LinearOperator; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; + + // discrete function to wrap what is used as Vector in eWoms + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpace > + VectorWrapperDiscreteFunction; + typedef Dune::Fem::PetscDiscreteFunction< DiscreteFunctionSpace > + PetscDiscreteFunctionType; + + enum { dimWorld = GridView::dimensionworld }; + +public: + AmgXSolverBackend(const Simulator& simulator) + : simulator_(simulator) + , amgxSolver_() + , rhs_( nullptr ) + , iterations_( 0 ) + { + std::string paramFileName = EWOMS_GET_PARAM(TypeTag, std::string, FemSolverParameterFileName); + if( paramFileName != "" ) + { + Dune::Fem::Parameter::append( paramFileName ); + } + } + + ~AmgXSolverBackend() + { cleanup_(); + if( A_ ) + { + ::Dune::Petsc::MatDestroy( A_.operator ->() ); + A_.reset(); + } + } + + /*! + * \brief Register all run-time parameters for the linear solver. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverTolerance, + "The maximum allowed error between of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverAbsTolerance, + "The maximum allowed error between of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIterations, + "The maximum number of iterations of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity, + "The verbosity level of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, LinearSolverOverlapSize, + "The size of the algebraic overlap for the linear solver"); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverMaxError, + "The maximum residual error which the linear solver tolerates" + " without giving up"); + + EWOMS_REGISTER_PARAM(TypeTag, int, PreconditionerOrder, + "The order of the preconditioner"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, PreconditionerRelaxation, + "The relaxation factor of the preconditioner"); + + EWOMS_REGISTER_PARAM(TypeTag, std::string, AmgxSolverMode, + "The name of the solver mode for AMGX"); + EWOMS_REGISTER_PARAM(TypeTag, std::string, AmgxSolverConfigFileName, + "The name of the file which contains the AMGX solver configuration"); + + EWOMS_REGISTER_PARAM(TypeTag, std::string, FemSolverParameterFileName, + "The name of the file which contains the parameters for the DUNE-FEM solvers"); + + } + + /*! + * \brief Causes the solve() method to discared the structure of the linear system of + * equations the next time it is called. + */ + void eraseMatrix() + { cleanup_(); } + + void prepare(const LinearOperator& op, Vector& b) + { + //Scalar linearSolverTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance); + //Scalar linearSolverAbsTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverAbsTolerance); + + if( ! amgxSolver_ ) + { + // reset linear solver + std::string mode = EWOMS_GET_PARAM(TypeTag, std::string, AmgxSolverMode); + std::string solverconfig = EWOMS_GET_PARAM(TypeTag, std::string, AmgxSolverConfigFileName); + + amgxSolver_.reset( new AmgXSolver() ); + amgxSolver_->initialize(PETSC_COMM_WORLD, mode, solverconfig); + + if(! A_) { + A_.reset( new Mat() ); + } + + // convert MATBAIJ to MATAIJ which is needed by AmgXSolver + ::Dune::Petsc::ErrorCheck( ::MatConvert( op.petscMatrix(), MATAIJ, MAT_INITIAL_MATRIX, A_.operator ->() ) ); + // set up the matrix used by AmgX + amgxSolver_->setA( *A_ ); + } + + // store pointer to right hand side + rhs_ = &b; + } + + /*! + * \brief Actually solve the linear system of equations. + * + * \return true if the residual reduction could be achieved, else false. + */ + bool solve(Vector& x) + { + // wrap x into discrete function X (no copy) + VectorWrapperDiscreteFunction X( "FSB::x", space(), x ); + VectorWrapperDiscreteFunction B( "FSB::rhs", space(), *rhs_ ); + + if( ! petscRhs_ ) + { + petscRhs_.reset( new PetscDiscreteFunctionType( "AMGX::rhs", space() ) ); + } + if( ! petscX_ ) + { + petscX_.reset( new PetscDiscreteFunctionType("AMGX::X", space()) ); + } + + petscRhs_->assign( B ); + petscX_->assign( X ); + + // solve with right hand side rhs and store in x + Vec& p = *(petscX_->petscVec()); + Vec& rhs = *(petscRhs_->petscVec()); + + try { + amgxSolver_->solve( p, rhs ); + } + catch (...) + { + OPM_THROW(Opm::NumericalIssue, "AmgXSolver: no convergence of linear solver!"); + } + + // copy result to ewoms solution + X.assign( *petscX_ ); + + amgxSolver_->getIters(iterations_); + + // return the result of the solver + return true; + } + + /*! + * \brief Return number of iterations used during last solve. + */ + size_t iterations () const { + return iterations_; + } + +protected: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + + const DiscreteFunctionSpace& space() const { + return simulator_.model().space(); + } + + void cleanup_() + { + if( amgxSolver_ ) + { + amgxSolver_->finalize(); + amgxSolver_.reset(); + } + + rhs_ = nullptr; + + petscRhs_.reset(); + petscX_.reset(); + } + + const Simulator& simulator_; + + std::unique_ptr< Mat > A_; + + std::unique_ptr< PetscDiscreteFunctionType > petscRhs_; + std::unique_ptr< PetscDiscreteFunctionType > petscX_; + + std::unique_ptr< AmgXSolver > amgxSolver_; + + Vector* rhs_; + int iterations_; +}; +}} // namespace Linear, Ewoms + +#endif // HAVE_DUNE_FEM + +#endif diff --git a/ewoms/linear/femsolverbackend.hh b/ewoms/linear/femsolverbackend.hh new file mode 100644 index 0000000000..460fa47bf3 --- /dev/null +++ b/ewoms/linear/femsolverbackend.hh @@ -0,0 +1,390 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM 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 2 of the License, or + (at your option) any later version. + + OPM 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 OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Ewoms::Linear::FemSolverBackend + */ +#ifndef EWOMS_FEM_SOLVER_BACKEND_HH +#define EWOMS_FEM_SOLVER_BACKEND_HH + +#include + +#if USE_DUNE_FEM_SOLVERS + +#define DISABLE_AMG_DIRECTSOLVER 1 +#include +#include +#include + +#if HAVE_VIENNACL +#include +#endif + + +#include +#include +#include + +#include + +#include + + +#include +#include + +#include +#include +#include + +namespace Ewoms { +namespace Linear { +template +class FemSolverBackend; +}} // namespace Linear, Ewoms + + +BEGIN_PROPERTIES + +NEW_TYPE_TAG(FemSolverBackend); + +SET_TYPE_PROP(FemSolverBackend, + LinearSolverBackend, + Ewoms::Linear::FemSolverBackend); + +//NEW_PROP_TAG(LinearSolverTolerance); +NEW_PROP_TAG(LinearSolverMaxIterations); +NEW_PROP_TAG(LinearSolverVerbosity); +NEW_PROP_TAG(LinearSolverMaxError); +NEW_PROP_TAG(LinearSolverOverlapSize); +//! The order of the sequential preconditioner +NEW_PROP_TAG(PreconditionerOrder); + +//! The relaxation factor of the preconditioner +NEW_PROP_TAG(PreconditionerRelaxation); + +//! Filename for DUNE-FEM solver parameters +NEW_PROP_TAG(FemSolverParameterFileName); + +//! make the linear solver shut up by default +SET_INT_PROP(FemSolverBackend, LinearSolverVerbosity, 0); + +//! set the default number of maximum iterations for the linear solver +SET_INT_PROP(FemSolverBackend, LinearSolverMaxIterations, 1000); + +SET_SCALAR_PROP(FemSolverBackend, LinearSolverMaxError, 1e7); + +//! set the default overlap size to 2 +SET_INT_PROP(FemSolverBackend, LinearSolverOverlapSize, 2); + +//! set the preconditioner order to 0 by default +SET_INT_PROP(FemSolverBackend, PreconditionerOrder, 0); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_SCALAR_PROP(FemSolverBackend, PreconditionerRelaxation, 1.0); + +//! set the preconditioner relaxation parameter to 1.0 by default +SET_STRING_PROP(FemSolverBackend, FemSolverParameterFileName, ""); + +//! make the linear solver shut up by default +//SET_SCALAR_PROP(FemSolverBackend, LinearSolverTolerance, 0.01); + +END_PROPERTIES + +namespace Ewoms { +namespace Linear { +/*! + * \ingroup Linear + * + * \brief Provides the common code which is required by most linear solvers. + * + * This class provides access to all preconditioners offered by dune-istl using the + * PreconditionerWrapper property: + * \code + * SET_TYPE_PROP(YourTypeTag, PreconditionerWrapper, + * Ewoms::Linear::PreconditionerWrapper$PRECONDITIONER); + * \endcode + * + * Where the choices possible for '\c $PRECONDITIONER' are: + * - \c Jacobi: A Jacobi preconditioner + * - \c GaussSeidel: A Gauss-Seidel preconditioner + * - \c SSOR: A symmetric successive overrelaxation (SSOR) preconditioner + * - \c SOR: A successive overrelaxation (SOR) preconditioner + * - \c ILUn: An ILU(n) preconditioner + * - \c ILU0: An ILU(0) preconditioner. The results of this + * preconditioner are the same as setting the + * PreconditionerOrder property to 0 and using the ILU(n) + * preconditioner. The reason for the existence of ILU0 is + * that it is computationally cheaper because it does not + * need to consider things which are only required for + * higher orders + */ +template +class FemSolverBackend +{ +protected: + typedef typename GET_PROP_TYPE(TypeTag, LinearSolverBackend) Implementation; + + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) LinearOperator; + typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) Vector; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; + + // discrete function to wrap what is used as Vector in eWoms + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpace > + VectorWrapperDiscreteFunction; + + template + struct SolverSelector + { + typedef Dune::Fem::KrylovInverseOperator< VectorWrapperDiscreteFunction > type; + }; + +#if HAVE_PETSC + template + struct SolverSelector< d, Dune::Fem::PetscLinearOperator< VectorWrapperDiscreteFunction, VectorWrapperDiscreteFunction > > + { + typedef Dune::Fem::PetscInverseOperator< VectorWrapperDiscreteFunction, LinearOperator > type; + }; +#endif + +#if HAVE_VIENNACL + template + struct SolverSelector< d, Dune::Fem::SparseRowLinearOperator< VectorWrapperDiscreteFunction, VectorWrapperDiscreteFunction > > + { + typedef Dune::Fem::ViennaCLInverseOperator< VectorWrapperDiscreteFunction > type; + }; +#endif + + template + struct SolverSelector< d, Dune::Fem::ISTLLinearOperator< VectorWrapperDiscreteFunction, VectorWrapperDiscreteFunction > > + { + typedef Dune::Fem::ISTLBICGSTABOp< VectorWrapperDiscreteFunction, LinearOperator > type; + }; + + // select solver type depending on linear operator type + typedef typename SolverSelector<0, typename LinearOperator::ParentType > :: type InverseLinearOperator; + + enum { dimWorld = GridView::dimensionworld }; + +public: + FemSolverBackend(const Simulator& simulator) + : simulator_(simulator) + , invOp_() + , rhs_( nullptr ) + { + std::string paramFileName = EWOMS_GET_PARAM(TypeTag, std::string, FemSolverParameterFileName); + if( paramFileName != "" ) + { + Dune::Fem::Parameter::append( paramFileName ); + } + else + { + // default parameters + Dune::Fem::Parameter::append("fem.solver.errormeasure", "residualreduction" ); + Dune::Fem::Parameter::append("fem.solver.verbose", "false" ); + + // Krylov solver + Dune::Fem::Parameter::append("fem.solver.method", "bicgstab" ); + + // Preconditioner + Dune::Fem::Parameter::append("fem.solver.preconditioning.method", "ilu" ); + Dune::Fem::Parameter::append("fem.solver.preconditioning.level", "0" ); + Dune::Fem::Parameter::append("fem.solver.preconditioning.relaxation", "0.9" ); + } + } + + ~FemSolverBackend() + { cleanup_(); } + + /*! + * \brief Register all run-time parameters for the linear solver. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverTolerance, + "The maximum allowed error between of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverMaxIterations, + "The maximum number of iterations of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, int, LinearSolverVerbosity, + "The verbosity level of the linear solver"); + EWOMS_REGISTER_PARAM(TypeTag, unsigned, LinearSolverOverlapSize, + "The size of the algebraic overlap for the linear solver"); + + EWOMS_REGISTER_PARAM(TypeTag, Scalar, LinearSolverMaxError, + "The maximum residual error which the linear solver tolerates" + " without giving up"); + + EWOMS_REGISTER_PARAM(TypeTag, int, PreconditionerOrder, + "The order of the preconditioner"); + EWOMS_REGISTER_PARAM(TypeTag, Scalar, PreconditionerRelaxation, + "The relaxation factor of the preconditioner"); + + EWOMS_REGISTER_PARAM(TypeTag, std::string, FemSolverParameterFileName, + "The name of the file which contains the parameters for the DUNE-FEM solvers"); + + } + + /*! + * \brief Causes the solve() method to discared the structure of the linear system of + * equations the next time it is called. + */ + void eraseMatrix() + { cleanup_(); } + + void prepare(const LinearOperator& op, Vector& b) + { + Scalar linearSolverTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance); + Scalar linearSolverAbsTolerance = this->simulator_.model().newtonMethod().tolerance() / 100000.0; + + // reset linear solver + if( ! invOp_ ) + { + invOp_.reset( new InverseLinearOperator( linearSolverTolerance, linearSolverAbsTolerance ) ); + } + + setMatrix( op ); + setResidual( b ); + + // not needed + asImp_().rescale_(); + } + + /*! + * \brief Assign values to the internal data structure for the residual vector. + * + * This method also cares about synchronizing that vector with the peer processes. + */ + void setResidual(const Vector& b) + { + // copy the interior values of the non-overlapping residual vector to the + // overlapping one + rhs_ = &b ; + } + + /*! + * \brief Retrieve the synchronized internal residual vector. + * + * This only deals with entries which are local to the current process. + */ + void getResidual(Vector& b) const + { + assert( rhs_ ); + // copy residual + b = *rhs_; + } + + /*! + * \brief Sets the values of the residual's Jacobian matrix. + * + * This method also synchronizes the data structure across the processes which are + * involved in the simulation run. + */ + void setMatrix(const SparseMatrixAdapter& op) + { + invOp_.bind( op ); + } + + + /*! + * \brief Actually solve the linear system of equations. + * + * \return true if the residual reduction could be achieved, else false. + */ + bool solve(Vector& x) + { + // wrap x into discrete function X (no copy) + VectorWrapperDiscreteFunction X( "FSB::x", space(), x ); + assert( rhs_ ); + VectorWrapperDiscreteFunction B( "FSB::rhs", space(), *rhs_ ); + + // solve with right hand side rhs and store in x + (*invOp_)( B, X ); + + // return the result of the solver + return invOp_->iterations() < 0 ? false : true; + } + + /*! + * \brief Return number of iterations used during last solve. + */ + size_t iterations () const { + assert( invOp_); + return invOp_->iterations(); + } + +protected: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + + const DiscreteFunctionSpace& space() const { + return simulator_.model().space(); + } + + void rescale_() + { + /* + const auto& overlap = overlappingMatrix_->overlap(); + for (unsigned domesticRowIdx = 0; domesticRowIdx < overlap.numLocal(); ++domesticRowIdx) { + Index nativeRowIdx = overlap.domesticToNative(static_cast(domesticRowIdx)); + auto& row = (*overlappingMatrix_)[domesticRowIdx]; + + auto colIt = row.begin(); + const auto& colEndIt = row.end(); + for (; colIt != colEndIt; ++ colIt) { + auto& entry = *colIt; + for (unsigned i = 0; i < entry.rows; ++i) + entry[i] *= simulator_.model().eqWeight(nativeRowIdx, i); + } + + auto& rhsEntry = (*overlappingb_)[domesticRowIdx]; + for (unsigned i = 0; i < rhsEntry.size(); ++i) + rhsEntry[i] *= simulator_.model().eqWeight(nativeRowIdx, i); + } + */ + } + + void cleanup_() + { + invOp_.reset(); + rhs_ = nullptr; + } + + const Simulator& simulator_; + + std::unique_ptr< InverseLinearOperator > invOp_; + + const Vector* rhs_; +}; +}} // namespace Linear, Ewoms + +#endif // HAVE_DUNE_FEM + +#endif diff --git a/ewoms/linear/overlappingbcrsmatrix.hh b/ewoms/linear/overlappingbcrsmatrix.hh index 9adb168476..0658378256 100644 --- a/ewoms/linear/overlappingbcrsmatrix.hh +++ b/ewoms/linear/overlappingbcrsmatrix.hh @@ -72,8 +72,15 @@ public: : ParentType(other) {} - template - OverlappingBCRSMatrix(const NativeBCRSMatrix& nativeMatrix, + template + OverlappingBCRSMatrix(const JacobianMatrix& jacobian, + const BorderList& borderList, + const BlackList& blackList, + unsigned overlapSize) + : OverlappingBCRSMatrix( jacobian.matrix(), borderList, blackList, overlapSize ) + {} + + OverlappingBCRSMatrix(const ParentType& nativeMatrix, const BorderList& borderList, const BlackList& blackList, unsigned overlapSize) @@ -151,8 +158,20 @@ public: * * The non-master entries are copied from the master */ - template - void assignCopy(const NativeBCRSMatrix& nativeMatrix) + template + void assignCopy(const JacobianMatrix& jacobian) + { + assignCopy( jacobian.matrix() ); + } + + + /*! + * \brief Assign and syncronize the overlapping matrix from a + * non-overlapping one. + * + * The non-master entries are copied from the master + */ + void assignCopy(const ParentType& nativeMatrix) { // copy the native entries assignFromNative(nativeMatrix); @@ -217,8 +236,13 @@ public: "row"); } - template - void assignFromNative(const NativeBCRSMatrix& nativeMatrix) + template + void assignFromNative(const JacobianMatrix& jacobian) + { + assignFromNative( jacobian.matrix() ); + } + + void assignFromNative(const ParentType& nativeMatrix) { // first, set everything to 0, BCRSMatrix::operator=(0.0); diff --git a/ewoms/linear/parallelbasebackend.hh b/ewoms/linear/parallelbasebackend.hh index 4f92900263..dbc2b20a1c 100644 --- a/ewoms/linear/parallelbasebackend.hh +++ b/ewoms/linear/parallelbasebackend.hh @@ -453,7 +453,7 @@ SET_PROP(ParallelBaseLinearSolver, OverlappingMatrix) private: static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq); typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar; - typedef Ewoms::MatrixBlock MatrixBlock; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) :: MatrixBlock MatrixBlock; typedef Dune::BCRSMatrix NonOverlappingMatrix; public: diff --git a/ewoms/nonlinear/newtonmethod.hh b/ewoms/nonlinear/newtonmethod.hh index f86dcfc309..15bd95a71d 100644 --- a/ewoms/nonlinear/newtonmethod.hh +++ b/ewoms/nonlinear/newtonmethod.hh @@ -371,9 +371,12 @@ public: asImp_().linearizeAuxiliaryEquations_(); linearizeTimer_.stop(); + // finalize the linearization process if not done already + linearizer.finalize(); solveTimer_.start(); auto& residual = linearizer.residual(); const auto& jacobian = linearizer.jacobian(); + linearSolver_.prepare(jacobian, residual); linearSolver_.setResidual(residual); linearSolver_.getResidual(residual); diff --git a/tests/co2injection_flash_ni_ecfv.cc b/tests/co2injection_flash_ni_ecfv.cc index 1072a9f272..bb1d343ffc 100644 --- a/tests/co2injection_flash_ni_ecfv.cc +++ b/tests/co2injection_flash_ni_ecfv.cc @@ -32,6 +32,7 @@ #include #include #include + #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh" diff --git a/tests/co2injection_flash_ni_vcfv.cc b/tests/co2injection_flash_ni_vcfv.cc index 8a890e2b00..2f43adcd1b 100644 --- a/tests/co2injection_flash_ni_vcfv.cc +++ b/tests/co2injection_flash_ni_vcfv.cc @@ -32,6 +32,7 @@ #include #include #include + #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh" diff --git a/tests/co2injection_flash_vcfv.cc b/tests/co2injection_flash_vcfv.cc index b7ac5842d0..4ca01b138f 100644 --- a/tests/co2injection_flash_vcfv.cc +++ b/tests/co2injection_flash_vcfv.cc @@ -35,6 +35,7 @@ #include #include #include + #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh"