From 22c59060ea04df72e4a7240493cb014cd0c95fd7 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Wed, 22 May 2019 17:41:11 +0200 Subject: [PATCH] [feature][solver] This adds a FemSolverBackend and the AMGXSolverBackend and various smaller fixes to accommodate this. --- ewoms-prereqs.cmake | 10 ++ ewoms/common/basicproperties.hh | 4 + ewoms/disc/common/fvbasediscretization.hh | 135 ++++++++++++++++++---- ewoms/disc/common/fvbaselinearizer.hh | 12 +- ewoms/disc/common/fvbaseproperties.hh | 15 ++- ewoms/disc/ecfv/ecfvdiscretization.hh | 18 +++ ewoms/disc/vcfv/vcfvdiscretization.hh | 18 +++ ewoms/linear/overlappingbcrsmatrix.hh | 36 +++++- ewoms/linear/parallelbasebackend.hh | 2 +- ewoms/nonlinear/newtonmethod.hh | 3 + tests/co2injection_flash_ni_ecfv.cc | 3 + tests/co2injection_flash_ni_vcfv.cc | 3 + tests/co2injection_flash_vcfv.cc | 3 + 13 files changed, 222 insertions(+), 40 deletions(-) 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/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..4f53584b15 100644 --- a/tests/co2injection_flash_ni_ecfv.cc +++ b/tests/co2injection_flash_ni_ecfv.cc @@ -32,6 +32,9 @@ #include #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..320fec9269 100644 --- a/tests/co2injection_flash_ni_vcfv.cc +++ b/tests/co2injection_flash_ni_vcfv.cc @@ -32,6 +32,9 @@ #include #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..68b5aed8c1 100644 --- a/tests/co2injection_flash_vcfv.cc +++ b/tests/co2injection_flash_vcfv.cc @@ -35,6 +35,9 @@ #include #include #include + +#include + #include "problems/co2injectionflash.hh" #include "problems/co2injectionproblem.hh"