diff --git a/ewoms/common/basicproperties.hh b/ewoms/common/basicproperties.hh index 2195f13635..680fc91ab1 100644 --- a/ewoms/common/basicproperties.hh +++ b/ewoms/common/basicproperties.hh @@ -78,10 +78,6 @@ 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 9f8221d3e3..4b8dbaa77a 100644 --- a/ewoms/disc/common/fvbasediscretization.hh +++ b/ewoms/disc/common/fvbasediscretization.hh @@ -50,7 +50,6 @@ #include #include #include -#include #include #include #include @@ -76,13 +75,7 @@ #include #include #include - -#if HAVE_PETSC -#include #endif -#include -#include -#endif // endif HAVE_DUNE_FEM #include #include @@ -136,78 +129,6 @@ 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 -#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: - typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; - enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; - typedef Ewoms::MatrixBlock Block; - -public: - typedef typename Ewoms::Linear::IstlSparseMatrixAdapter type; -}; -#endif - //! The maximum allowed number of timestep divisions for the //! Newton solver SET_INT_PROP(FvBaseDiscretization, MaxTimeStepDivisions, 10); @@ -408,15 +329,13 @@ 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 DiscreteFunctionSpace& space) - : blockVector_(space.size()) + BlockVectorWrapper(const std::string& name OPM_UNUSED, const size_t size) + : blockVector_(size) {} SolutionVector& blockVector() @@ -426,12 +345,14 @@ class FvBaseDiscretization }; #if HAVE_DUNE_FEM + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; + // discrete function storing solution data - typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction 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; @@ -440,6 +361,7 @@ 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 @@ -459,14 +381,14 @@ public: , elementMapper_(gridView_) , vertexMapper_(gridView_) #endif + , newtonMethod_(simulator) + , localLinearizer_(ThreadManager::maxThreads()) + , linearizer_(new Linearizer()) #if HAVE_DUNE_FEM - , discreteFunctionSpace_( simulator.vanguard().gridPart() ) + , space_( simulator.vanguard().gridPart() ) #else - , discreteFunctionSpace_( asImp_().numGridDof() ) + , space_( 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)) @@ -491,7 +413,7 @@ public: size_t numDof = asImp_().numGridDof(); for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) { - solution_[timeIdx].reset(new DiscreteFunction("solution", discreteFunctionSpace_)); + solution_[timeIdx].reset(new DiscreteFunction("solution", space_)); if (storeIntensiveQuantities()) { intensiveQuantityCache_[timeIdx].resize(numDof); @@ -1156,16 +1078,6 @@ 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 @@ -1583,7 +1495,7 @@ public: void resetLinearizer () { delete linearizer_; - linearizer_ = new Linearizer(); + linearizer_ = new Linearizer; linearizer_->init(simulator_); } @@ -1817,11 +1729,6 @@ public: solution(timeIdx).resize(numDof); auxMod->applyInitial(); - -#if DUNE_VERSION_NEWER( DUNE_FEM, 2, 7 ) - discreteFunctionSpace_.extendSize( asImp_().numAuxiliaryDof() ); -#endif - } /*! @@ -1888,11 +1795,6 @@ public: const Ewoms::Timer& updateTimer() const { return updateTimer_; } -#if HAVE_DUNE_FEM - const DiscreteFunctionSpace& space() const { return discreteFunctionSpace_; } -#endif - - protected: void resizeAndResetIntensiveQuantitiesCache_() { @@ -1963,18 +1865,9 @@ 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_; @@ -1993,6 +1886,15 @@ 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 6a055e6842..dfe0ed4b80 100644 --- a/ewoms/disc/common/fvbaselinearizer.hh +++ b/ewoms/disc/common/fvbaselinearizer.hh @@ -54,7 +54,6 @@ namespace Ewoms { template class EcfvDiscretization; - /*! * \ingroup FiniteVolumeDiscretizations * diff --git a/ewoms/disc/common/fvbaseproperties.hh b/ewoms/disc/common/fvbaseproperties.hh index d372b7483b..e1bf3a6570 100644 --- a/ewoms/disc/common/fvbaseproperties.hh +++ b/ewoms/disc/common/fvbaseproperties.hh @@ -37,9 +37,6 @@ #include #include #include -#include -#include -#include BEGIN_PROPERTIES @@ -57,19 +54,10 @@ 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); @@ -92,6 +80,9 @@ 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 b2a789c769..9ae5805e5b 100644 --- a/ewoms/disc/ecfv/ecfvdiscretization.hh +++ b/ewoms/disc/ecfv/ecfvdiscretization.hh @@ -91,24 +91,6 @@ 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 042c194ce6..47d25fb3a2 100644 --- a/ewoms/disc/vcfv/vcfvdiscretization.hh +++ b/ewoms/disc/vcfv/vcfvdiscretization.hh @@ -101,24 +101,6 @@ 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 index 40af005968..58acca222d 100644 --- a/ewoms/linear/amgxsolverbackend.hh +++ b/ewoms/linear/amgxsolverbackend.hh @@ -31,7 +31,7 @@ #if USE_AMGX_SOLVERS -#if ! HAVE_PETSC || ! HAVE_AMGXSOLVER +#if !HAVE_PETSC || !HAVE_AMGXSOLVER #error "PETSc and AmgXSolver is needed for the AMGX solver backend" #endif @@ -130,28 +130,7 @@ 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 + * \brief Provides a linear solver backend that utilizes the AMG-X package. */ template class AmgXSolverBackend @@ -166,37 +145,36 @@ protected: typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace; - typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; // discrete function to wrap what is used as Vector in eWoms - typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpace > + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction VectorWrapperDiscreteFunction; - typedef Dune::Fem::PetscDiscreteFunction< DiscreteFunctionSpace > + typedef Dune::Fem::PetscDiscreteFunction PetscDiscreteFunctionType; enum { dimWorld = GridView::dimensionworld }; public: - AmgXSolverBackend(const Simulator& simulator) + AmgXSolverBackend(Simulator& simulator) : simulator_(simulator) + , space_(simulator.vanguard().gridPart()) , amgxSolver_() - , rhs_( nullptr ) - , iterations_( 0 ) + , rhs_(nullptr) + , iterations_(0) { std::string paramFileName = EWOMS_GET_PARAM(TypeTag, std::string, FemSolverParameterFileName); - if( paramFileName != "" ) - { - Dune::Fem::Parameter::append( paramFileName ); - } + if (paramFileName != "") + Dune::Fem::Parameter::append(paramFileName); } ~AmgXSolverBackend() - { cleanup_(); - if( A_ ) - { - ::Dune::Petsc::MatDestroy( A_.operator ->() ); - A_.reset(); - } + { + cleanup_(); + if (A_) { + ::Dune::Petsc::MatDestroy(A_.operator ->()); + A_.reset(); + } } /*! @@ -246,23 +224,22 @@ public: //Scalar linearSolverTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverTolerance); //Scalar linearSolverAbsTolerance = EWOMS_GET_PARAM(TypeTag, Scalar, LinearSolverAbsTolerance); - if( ! amgxSolver_ ) - { + 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_.reset(new AmgXSolver()); amgxSolver_->initialize(PETSC_COMM_WORLD, mode, solverconfig); - if(! A_) { - A_.reset( new Mat() ); + 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 ->() ) ); + ::Dune::Petsc::ErrorCheck(::MatConvert(op.petscMatrix(), MATAIJ, MAT_INITIAL_MATRIX, A_.operator ->())); // set up the matrix used by AmgX - amgxSolver_->setA( *A_ ); + amgxSolver_->setA(*A_); } // store pointer to right hand side @@ -277,35 +254,31 @@ public: bool solve(Vector& x) { // wrap x into discrete function X (no copy) - VectorWrapperDiscreteFunction X( "FSB::x", space(), x ); - VectorWrapperDiscreteFunction B( "FSB::rhs", space(), *rhs_ ); + 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()) ); - } + if (!petscRhs_) + petscRhs_.reset(new PetscDiscreteFunctionType("AMGX::rhs", space_)); + + if (!petscX_) + petscX_.reset(new PetscDiscreteFunctionType("AMGX::X", space_)); - petscRhs_->assign( B ); - petscX_->assign( X ); + 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 ); + amgxSolver_->solve(p, rhs); } - catch (...) - { + catch (...) { OPM_THROW(Opm::NumericalIssue, "AmgXSolver: no convergence of linear solver!"); } // copy result to ewoms solution - X.assign( *petscX_ ); + X.assign(*petscX_); amgxSolver_->getIters(iterations_); @@ -327,14 +300,9 @@ protected: const Implementation& asImp_() const { return *static_cast(this); } - const DiscreteFunctionSpace& space() const { - return simulator_.model().space(); - } - void cleanup_() { - if( amgxSolver_ ) - { + if (amgxSolver_) { amgxSolver_->finalize(); amgxSolver_.reset(); } @@ -347,12 +315,13 @@ protected: const Simulator& simulator_; - std::unique_ptr< Mat > A_; + DiscreteFunctionSpace space_; + std::unique_ptr A_; - std::unique_ptr< PetscDiscreteFunctionType > petscRhs_; - std::unique_ptr< PetscDiscreteFunctionType > petscX_; + std::unique_ptr petscRhs_; + std::unique_ptr petscX_; - std::unique_ptr< AmgXSolver > amgxSolver_; + std::unique_ptr amgxSolver_; Vector* rhs_; int iterations_; diff --git a/ewoms/linear/femsolverbackend.hh b/ewoms/linear/femsolverbackend.hh index 460fa47bf3..3d54edbf6f 100644 --- a/ewoms/linear/femsolverbackend.hh +++ b/ewoms/linear/femsolverbackend.hh @@ -29,30 +29,24 @@ #include -#if USE_DUNE_FEM_SOLVERS - #define DISABLE_AMG_DIRECTSOLVER 1 #include -#include +//#include #include #if HAVE_VIENNACL #include #endif - #include #include #include +#include +#include #include - #include - -#include -#include - #include #include #include @@ -72,6 +66,8 @@ SET_TYPE_PROP(FemSolverBackend, LinearSolverBackend, Ewoms::Linear::FemSolverBackend); +NEW_PROP_TAG(DiscreteFunction); + //NEW_PROP_TAG(LinearSolverTolerance); NEW_PROP_TAG(LinearSolverMaxIterations); NEW_PROP_TAG(LinearSolverVerbosity); @@ -109,6 +105,64 @@ SET_STRING_PROP(FemSolverBackend, FemSolverParameterFileName, ""); //! make the linear solver shut up by default //SET_SCALAR_PROP(FemSolverBackend, LinearSolverTolerance, 0.01); +SET_PROP(FemSolverBackend, 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; +}; + +SET_PROP(FemSolverBackend, 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 LinearOperator; +#elif USE_DUNE_FEM_VIENNACL_SOLVERS +#warning "Using Dune-Fem ViennaCL solvers" + typedef Dune::Fem::SparseRowLinearOperator LinearOperator; +#else +#warning "Using Dune-Fem ISTL solvers" + typedef Dune::Fem::ISTLLinearOperator LinearOperator; +#endif + + struct FemMatrixBackend : public LinearOperator + { + typedef LinearOperator ParentType; + typedef typename LinearOperator::MatrixType Matrix; + typedef typename ParentType::MatrixBlockType MatrixBlock; + template + FemMatrixBackend(Simulator& simulator) + : LinearOperator("eWoms::Jacobian", space_, space_) + , space_(simulator.vanguard().gridPart()) + {} + + void commit() + { this->flushAssembly(); } + + template + 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); } + + DiscreteFunctionSpace space_; + }; + +public: + typedef FemMatrixBackend type; +}; + + END_PROPERTIES namespace Ewoms { @@ -116,28 +170,7 @@ 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 + * \brief Uses the dune-fem infrastructure to solve the linear system of equations. */ template class FemSolverBackend @@ -147,74 +180,74 @@ protected: 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, SparseMatrixAdapter) LinearOperator; + typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter; 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; + typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction) DiscreteFunction; // discrete function to wrap what is used as Vector in eWoms - typedef Dune::Fem::ISTLBlockVectorDiscreteFunction< DiscreteFunctionSpace > + typedef Dune::Fem::ISTLBlockVectorDiscreteFunction VectorWrapperDiscreteFunction; template struct SolverSelector { - typedef Dune::Fem::KrylovInverseOperator< VectorWrapperDiscreteFunction > type; + typedef Dune::Fem::KrylovInverseOperator type; }; #if HAVE_PETSC template - struct SolverSelector< d, Dune::Fem::PetscLinearOperator< VectorWrapperDiscreteFunction, VectorWrapperDiscreteFunction > > + struct SolverSelector> { - typedef Dune::Fem::PetscInverseOperator< VectorWrapperDiscreteFunction, LinearOperator > type; + typedef Dune::Fem::PetscInverseOperator type; }; #endif #if HAVE_VIENNACL template - struct SolverSelector< d, Dune::Fem::SparseRowLinearOperator< VectorWrapperDiscreteFunction, VectorWrapperDiscreteFunction > > + struct SolverSelector > { - typedef Dune::Fem::ViennaCLInverseOperator< VectorWrapperDiscreteFunction > type; + typedef Dune::Fem::ViennaCLInverseOperator type; }; #endif template - struct SolverSelector< d, Dune::Fem::ISTLLinearOperator< VectorWrapperDiscreteFunction, VectorWrapperDiscreteFunction > > + struct SolverSelector > { - typedef Dune::Fem::ISTLBICGSTABOp< VectorWrapperDiscreteFunction, LinearOperator > type; + typedef Dune::Fem::ISTLBICGSTABOp type; }; // select solver type depending on linear operator type - typedef typename SolverSelector<0, typename LinearOperator::ParentType > :: type InverseLinearOperator; + typedef typename LinearOperator::ParentType Bla; + typedef typename SolverSelector<0, Bla>::type InverseLinearOperator; enum { dimWorld = GridView::dimensionworld }; public: - FemSolverBackend(const Simulator& simulator) + FemSolverBackend(Simulator& simulator) : simulator_(simulator) , invOp_() - , rhs_( nullptr ) + , rhs_(nullptr) + , space_(simulator.vanguard().gridPart()) { std::string paramFileName = EWOMS_GET_PARAM(TypeTag, std::string, FemSolverParameterFileName); - if( paramFileName != "" ) - { - Dune::Fem::Parameter::append( paramFileName ); - } - else - { + 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" ); + 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" ); + 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" ); + 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"); } } @@ -262,13 +295,11 @@ public: Scalar linearSolverAbsTolerance = this->simulator_.model().newtonMethod().tolerance() / 100000.0; // reset linear solver - if( ! invOp_ ) - { - invOp_.reset( new InverseLinearOperator( linearSolverTolerance, linearSolverAbsTolerance ) ); - } + if (!invOp_) + invOp_.reset(new InverseLinearOperator(linearSolverTolerance, linearSolverAbsTolerance)); - setMatrix( op ); - setResidual( b ); + setMatrix(op); + setResidual(b); // not needed asImp_().rescale_(); @@ -280,11 +311,7 @@ public: * 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 ; - } + { rhs_ = &b ; } /*! * \brief Retrieve the synchronized internal residual vector. @@ -293,8 +320,7 @@ public: */ void getResidual(Vector& b) const { - assert( rhs_ ); - // copy residual + assert(rhs_); b = *rhs_; } @@ -306,7 +332,7 @@ public: */ void setMatrix(const SparseMatrixAdapter& op) { - invOp_.bind( op ); + invOp_->bind(op); } @@ -318,12 +344,12 @@ public: 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_ ); + 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 ); + (*invOp_)(B, X); // return the result of the solver return invOp_->iterations() < 0 ? false : true; @@ -332,8 +358,9 @@ public: /*! * \brief Return number of iterations used during last solve. */ - size_t iterations () const { - assert( invOp_); + size_t iterations () const + { + assert(invOp_); return invOp_->iterations(); } @@ -344,32 +371,8 @@ protected: 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_() { @@ -378,13 +381,10 @@ protected: } const Simulator& simulator_; - - std::unique_ptr< InverseLinearOperator > invOp_; - + std::unique_ptr invOp_; const Vector* rhs_; + DiscreteFunctionSpace space_; }; }} // namespace Linear, Ewoms -#endif // HAVE_DUNE_FEM - #endif diff --git a/ewoms/linear/overlappingbcrsmatrix.hh b/ewoms/linear/overlappingbcrsmatrix.hh index 0658378256..9adb168476 100644 --- a/ewoms/linear/overlappingbcrsmatrix.hh +++ b/ewoms/linear/overlappingbcrsmatrix.hh @@ -72,15 +72,8 @@ public: : ParentType(other) {} - template - OverlappingBCRSMatrix(const JacobianMatrix& jacobian, - const BorderList& borderList, - const BlackList& blackList, - unsigned overlapSize) - : OverlappingBCRSMatrix( jacobian.matrix(), borderList, blackList, overlapSize ) - {} - - OverlappingBCRSMatrix(const ParentType& nativeMatrix, + template + OverlappingBCRSMatrix(const NativeBCRSMatrix& nativeMatrix, const BorderList& borderList, const BlackList& blackList, unsigned overlapSize) @@ -158,20 +151,8 @@ public: * * The non-master entries are copied from the master */ - 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) + template + void assignCopy(const NativeBCRSMatrix& nativeMatrix) { // copy the native entries assignFromNative(nativeMatrix); @@ -236,13 +217,8 @@ public: "row"); } - template - void assignFromNative(const JacobianMatrix& jacobian) - { - assignFromNative( jacobian.matrix() ); - } - - void assignFromNative(const ParentType& nativeMatrix) + template + void assignFromNative(const NativeBCRSMatrix& nativeMatrix) { // first, set everything to 0, BCRSMatrix::operator=(0.0); diff --git a/ewoms/linear/parallelbasebackend.hh b/ewoms/linear/parallelbasebackend.hh index dbc2b20a1c..b12c653cec 100644 --- a/ewoms/linear/parallelbasebackend.hh +++ b/ewoms/linear/parallelbasebackend.hh @@ -27,6 +27,7 @@ #ifndef EWOMS_PARALLEL_BASE_BACKEND_HH #define EWOMS_PARALLEL_BASE_BACKEND_HH +#include #include #include #include @@ -114,6 +115,20 @@ NEW_PROP_TAG(PreconditionerOrder); //! The relaxation factor of the preconditioner NEW_PROP_TAG(PreconditionerRelaxation); + +//! Set the type of a global jacobian matrix for linear solvers that are based on +//! dune-istl. +SET_PROP(ParallelBaseLinearSolver, SparseMatrixAdapter) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + typedef Ewoms::MatrixBlock Block; + +public: + typedef typename Ewoms::Linear::IstlSparseMatrixAdapter type; +}; + END_PROPERTIES namespace Ewoms { @@ -453,7 +468,7 @@ SET_PROP(ParallelBaseLinearSolver, OverlappingMatrix) private: static constexpr int numEq = GET_PROP_VALUE(TypeTag, NumEq); typedef typename GET_PROP_TYPE(TypeTag, LinearSolverScalar) LinearSolverScalar; - typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) :: 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 a4fc83149b..c5da6638eb 100644 --- a/ewoms/nonlinear/newtonmethod.hh +++ b/ewoms/nonlinear/newtonmethod.hh @@ -373,6 +373,8 @@ public: linearizer.finalize(); linearizeTimer_.stop(); + // finalize the linearization process if not done already + linearizer.finalize(); solveTimer_.start(); auto& residual = linearizer.residual(); const auto& jacobian = linearizer.jacobian(); diff --git a/tests/lens_immiscible_ecfv_ad.hh b/tests/lens_immiscible_ecfv_ad.hh index 276a77e52c..bad1a7e0d2 100644 --- a/tests/lens_immiscible_ecfv_ad.hh +++ b/tests/lens_immiscible_ecfv_ad.hh @@ -43,8 +43,9 @@ SET_TAG_PROP(LensProblemEcfvAd, SpatialDiscretizationSplice, EcfvDiscretization) // use automatic differentiation for this simulator SET_TAG_PROP(LensProblemEcfvAd, LocalLinearizerSplice, AutoDiffLocalLinearizer); -// this problem works fine if the linear solver uses single precision scalars -SET_TYPE_PROP(LensProblemEcfvAd, LinearSolverScalar, float); +#warning WIP: only for testing! +// use the dune-fem linear solver backend +SET_TAG_PROP(LensProblemEcfvAd, LinearSolverSplice, FemSolverBackend); END_PROPERTIES diff --git a/tests/problems/lensproblem.hh b/tests/problems/lensproblem.hh index 766f2744f6..b573e93d73 100644 --- a/tests/problems/lensproblem.hh +++ b/tests/problems/lensproblem.hh @@ -28,6 +28,7 @@ #ifndef EWOMS_LENS_PROBLEM_HH #define EWOMS_LENS_PROBLEM_HH +#include #include #include #include