From c2c79f09032388100ea9c063f6c090b3bbacadf7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 24 Jun 2020 15:55:43 +0200 Subject: [PATCH] Move well operators to separate file. Also introduce new class WellModelAsLinearOperator making a well model into an actual Dune::LinearOperator, this prevents the TypeTag dependent type from leaking into the type of the WellModelMatrixAdapter instantiation. As a side benefit, the adapter classes can now adapt (i.e. combine with a matrix operator) any linear operator. --- CMakeLists_files.cmake | 1 + opm/simulators/linalg/ISTLSolverEbos.hpp | 178 +--------------- opm/simulators/linalg/WellOperators.hpp | 249 +++++++++++++++++++++++ 3 files changed, 258 insertions(+), 170 deletions(-) create mode 100644 opm/simulators/linalg/WellOperators.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index e963fe54f..b579cc7cc 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -172,6 +172,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/linalg/PressureTransferPolicy.hpp opm/simulators/linalg/PreconditionerFactory.hpp opm/simulators/linalg/PreconditionerWithUpdate.hpp + opm/simulators/linalg/WellOperators.hpp opm/simulators/linalg/WriteSystemMatrixHelper.hpp opm/simulators/linalg/findOverlapRowsAndColumns.hpp opm/simulators/linalg/getQuasiImpesWeights.hpp diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 6815b298e..75d9ea855 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -21,6 +21,7 @@ #ifndef OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED #define OPM_ISTLSOLVER_EBOS_HEADER_INCLUDED +#include #include #include #include @@ -88,165 +89,6 @@ DenseMatrix transposeDenseMatrix(const DenseMatrix& M) return tmp; } -//===================================================================== -// Implementation for ISTL-matrix based operator -//===================================================================== - - -/*! - \brief Adapter to turn a matrix into a linear operator. - - Adapts a matrix to the assembled linear operator interface - */ -template -class WellModelMatrixAdapter : public Dune::AssembledLinearOperator -{ - typedef Dune::AssembledLinearOperator BaseType; - -public: - typedef M matrix_type; - typedef X domain_type; - typedef Y range_type; - typedef typename X::field_type field_type; - -#if HAVE_MPI - typedef Dune::OwnerOverlapCopyCommunication communication_type; -#else - typedef Dune::CollectiveCommunication< int > communication_type; -#endif - - Dune::SolverCategory::Category category() const override - { - return overlapping ? - Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential; - } - - //! constructor: just store a reference to a matrix - WellModelMatrixAdapter (const M& A, - const WellModel& wellMod, - const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >()) - : A_( A ), wellMod_( wellMod ), comm_(comm) - {} - - - virtual void apply( const X& x, Y& y ) const override - { - A_.mv( x, y ); - - // add well model modification to y - wellMod_.apply(x, y ); - -#if HAVE_MPI - if( comm_ ) - comm_->project( y ); -#endif - } - - // y += \alpha * A * x - virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override - { - A_.usmv(alpha,x,y); - - // add scaled well model modification to y - wellMod_.applyScaleAdd( alpha, x, y ); - -#if HAVE_MPI - if( comm_ ) - comm_->project( y ); -#endif - } - - virtual const matrix_type& getmat() const override { return A_; } - -protected: - const matrix_type& A_ ; - const WellModel& wellMod_; - std::shared_ptr< communication_type > comm_; -}; - - -/*! - \brief Adapter to turn a matrix into a linear operator. - Adapts a matrix to the assembled linear operator interface. - We assume parallel ordering, where ghost rows are located after interior rows - */ -template -class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator -{ - typedef Dune::AssembledLinearOperator BaseType; - -public: - typedef M matrix_type; - typedef X domain_type; - typedef Y range_type; - typedef typename X::field_type field_type; - -#if HAVE_MPI - typedef Dune::OwnerOverlapCopyCommunication communication_type; -#else - typedef Dune::CollectiveCommunication< int > communication_type; -#endif - - - Dune::SolverCategory::Category category() const override - { - return overlapping ? - Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential; - } - - //! constructor: just store a reference to a matrix - WellModelGhostLastMatrixAdapter (const M& A, - const WellModel& wellMod, - const size_t interiorSize ) - : A_( A ), wellMod_( wellMod ), interiorSize_(interiorSize) - {} - - virtual void apply( const X& x, Y& y ) const override - { - for (auto row = A_.begin(); row.index() < interiorSize_; ++row) - { - y[row.index()]=0; - auto endc = (*row).end(); - for (auto col = (*row).begin(); col != endc; ++col) - (*col).umv(x[col.index()], y[row.index()]); - } - - // add well model modification to y - wellMod_.apply(x, y ); - - ghostLastProject( y ); - } - - // y += \alpha * A * x - virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override - { - for (auto row = A_.begin(); row.index() < interiorSize_; ++row) - { - auto endc = (*row).end(); - for (auto col = (*row).begin(); col != endc; ++col) - (*col).usmv(alpha, x[col.index()], y[row.index()]); - } - // add scaled well model modification to y - wellMod_.applyScaleAdd( alpha, x, y ); - - ghostLastProject( y ); - } - - virtual const matrix_type& getmat() const override { return A_; } - -protected: - void ghostLastProject(Y& y) const - { - size_t end = y.size(); - for (size_t i = interiorSize_; i < end; ++i) - y[i] = 0; - } - - const matrix_type& A_ ; - const WellModel& wellMod_; - size_t interiorSize_; -}; - /// This class solves the fully implicit black-oil system by /// solving the reduced system (after eliminating well variables) /// as a block-structured matrix (one block for all cell variables) for a fixed @@ -503,31 +345,27 @@ protected: } const WellModel& wellModel = simulator_.problem().wellModel(); + const WellModelAsLinearOperator wellOp(wellModel); if( isParallel() ) { if ( ownersFirst_ && !parameters_.linear_solver_use_amg_ && !useFlexible_) { - typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator; + typedef WellModelGhostLastMatrixAdapter< Matrix, Vector, Vector, true > Operator; assert(matrix_); - Operator opA(getMatrix(), wellModel, interiorCellNum_); - + Operator opA(getMatrix(), wellOp, interiorCellNum_); solve( opA, x, *rhs_, *comm_ ); } else { - - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, true > Operator; + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, true > Operator; assert (noGhostMat_); - Operator opA(getMatrix(), wellModel, - comm_ ); - + Operator opA(getMatrix(), wellOp, comm_ ); solve( opA, x, *rhs_, *comm_ ); - } } else { - typedef WellModelMatrixAdapter< Matrix, Vector, Vector, WellModel, false > Operator; - Operator opA(getMatrix(), wellModel); + typedef WellModelMatrixAdapter< Matrix, Vector, Vector, false > Operator; + Operator opA(getMatrix(), wellOp); solve( opA, x, *rhs_ ); } diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp new file mode 100644 index 000000000..f161f0dc2 --- /dev/null +++ b/opm/simulators/linalg/WellOperators.hpp @@ -0,0 +1,249 @@ +/* + Copyright 2016 IRIS AS + Copyright 2019, 2020 Equinor ASA + Copyright 2020 SINTEF + + 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 3 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 . +*/ + +#ifndef OPM_WELLOPERATORS_HEADER_INCLUDED +#define OPM_WELLOPERATORS_HEADER_INCLUDED + +#include + + +namespace Opm +{ + + +//===================================================================== +// Implementation for ISTL-matrix based operators +// Note: the classes WellModelMatrixAdapter and +// WellModelGhostLastMatrixAdapter were moved from ISTLSolverEbos.hpp +// and subsequently modified. +//===================================================================== + +/// Linear operator wrapper for well model. +/// +/// This class is intended to hide the actual type of the well model +/// (which depends on a TypeTag) by making a simple linear operator +/// wrapper. That way the WellModelMatrixAdapter template does not need +/// the concrete WellModel type, and we can avoid instantiating +/// WellModelMatrixAdapter separately for each TypeTag, as it will only +/// depend on the matrix and vector types involved, which typically are +/// just one for each block size with block sizes 1-4. +template +class WellModelAsLinearOperator : public Dune::LinearOperator +{ +public: + using Base = Dune::LinearOperator; + using field_type = typename Base::field_type; + + explicit WellModelAsLinearOperator(const WellModel& wm) + : wellMod_(wm) + { + } + /*! \brief apply operator to x: \f$ y = A(x) \f$ + The input vector is consistent and the output must also be + consistent on the interior+border partition. + */ + void apply (const X& x, Y& y) const override + { + wellMod_.apply(x, y ); + } + + //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$ + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override + { + wellMod_.applyScaleAdd( alpha, x, y ); + } + + /// Category for operator. + /// This is somewhat tricky, I consider this operator sequential + /// since (unlike WellModelMatrixAdapter) it does not do any + /// projections etc. but only forwards the calls to the sequential + /// well model. + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::sequential; + } +private: + const WellModel& wellMod_; +}; + +/*! + \brief Adapter to combine a matrix and another linear operator into + a combined linear operator. + + Adapts a matrix A plus another linear operator W (typically from + wells) to the assembled linear operator interface by returning S + from getmat() and making apply() and applyscaleadd() apply both A + and W to the input vector. In addition this is a parallel-aware + adapter, that does not require the W operator to be parallel, but + makes it into one by making the proper projections. + */ +template +class WellModelMatrixAdapter : public Dune::AssembledLinearOperator +{ +public: + typedef M matrix_type; + typedef X domain_type; + typedef Y range_type; + typedef typename X::field_type field_type; + +#if HAVE_MPI + typedef Dune::OwnerOverlapCopyCommunication communication_type; +#else + typedef Dune::CollectiveCommunication< int > communication_type; +#endif + + Dune::SolverCategory::Category category() const override + { + return overlapping ? + Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential; + } + + //! constructor: just store a reference to a matrix + WellModelMatrixAdapter (const M& A, + const Dune::LinearOperator& wellOper, + const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >()) + : A_( A ), wellOper_( wellOper ), comm_(comm) + {} + + + virtual void apply( const X& x, Y& y ) const override + { + A_.mv( x, y ); + + // add well model modification to y + wellOper_.apply(x, y ); + +#if HAVE_MPI + if( comm_ ) + comm_->project( y ); +#endif + } + + // y += \alpha * A * x + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override + { + A_.usmv(alpha,x,y); + + // add scaled well model modification to y + wellOper_.applyscaleadd( alpha, x, y ); + +#if HAVE_MPI + if( comm_ ) + comm_->project( y ); +#endif + } + + virtual const matrix_type& getmat() const override { return A_; } + +protected: + const matrix_type& A_ ; + const Dune::LinearOperator& wellOper_; + std::shared_ptr< communication_type > comm_; +}; + + +/*! + \brief Adapter to combine a matrix and another linear operator into + a combined linear operator. + + This is similar to WellModelMatrixAdapter, with the difference that + here we assume a parallel ordering of rows, where ghost rows are + located after interior rows. + */ +template +class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator +{ +public: + typedef M matrix_type; + typedef X domain_type; + typedef Y range_type; + typedef typename X::field_type field_type; + +#if HAVE_MPI + typedef Dune::OwnerOverlapCopyCommunication communication_type; +#else + typedef Dune::CollectiveCommunication< int > communication_type; +#endif + + + Dune::SolverCategory::Category category() const override + { + return overlapping ? + Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential; + } + + //! constructor: just store a reference to a matrix + WellModelGhostLastMatrixAdapter (const M& A, + const Dune::LinearOperator& wellOper, + const size_t interiorSize ) + : A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize) + {} + + virtual void apply( const X& x, Y& y ) const override + { + for (auto row = A_.begin(); row.index() < interiorSize_; ++row) + { + y[row.index()]=0; + auto endc = (*row).end(); + for (auto col = (*row).begin(); col != endc; ++col) + (*col).umv(x[col.index()], y[row.index()]); + } + + // add well model modification to y + wellOper_.apply(x, y ); + + ghostLastProject( y ); + } + + // y += \alpha * A * x + virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override + { + for (auto row = A_.begin(); row.index() < interiorSize_; ++row) + { + auto endc = (*row).end(); + for (auto col = (*row).begin(); col != endc; ++col) + (*col).usmv(alpha, x[col.index()], y[row.index()]); + } + // add scaled well model modification to y + wellOper_.applyscaleadd( alpha, x, y ); + + ghostLastProject( y ); + } + + virtual const matrix_type& getmat() const override { return A_; } + +protected: + void ghostLastProject(Y& y) const + { + size_t end = y.size(); + for (size_t i = interiorSize_; i < end; ++i) + y[i] = 0; + } + + const matrix_type& A_ ; + const Dune::LinearOperator& wellOper_; + size_t interiorSize_; +}; + +} // namespace Opm + +#endif // OPM_WELLOPERATORS_HEADER_INCLUDED +