My Project
WellOperators.hpp
1 /*
2  Copyright 2016 IRIS AS
3  Copyright 2019, 2020 Equinor ASA
4  Copyright 2020 SINTEF
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #ifndef OPM_WELLOPERATORS_HEADER_INCLUDED
23 #define OPM_WELLOPERATORS_HEADER_INCLUDED
24 
25 #include <dune/istl/operators.hh>
26 #include <opm/simulators/linalg/matrixblock.hh>
27 
28 namespace Opm
29 {
30 
31 
32 //=====================================================================
33 // Implementation for ISTL-matrix based operators
34 // Note: the classes WellModelMatrixAdapter and
35 // WellModelGhostLastMatrixAdapter were moved from ISTLSolverEbos.hpp
36 // and subsequently modified.
37 //=====================================================================
38 
48 template <class X, class Y>
49 class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
50 {
51 public:
52  using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
53  virtual void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const = 0;
54  virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian) const = 0;
55  virtual int getNumberOfExtraEquations() const = 0;
56 };
57 
58 template <class WellModel, class X, class Y>
60 {
61 public:
63  using field_type = typename Base::field_type;
64  using PressureMatrix = typename Base::PressureMatrix;
65  explicit WellModelAsLinearOperator(const WellModel& wm)
66  : wellMod_(wm)
67  {
68  }
73  void apply(const X& x, Y& y) const override
74  {
75  wellMod_.apply(x, y);
76  }
77 
79  virtual void applyscaleadd(field_type alpha, const X& x, Y& y) const override
80  {
81  wellMod_.applyScaleAdd(alpha, x, y);
82  }
83 
89  Dune::SolverCategory::Category category() const override
90  {
91  return Dune::SolverCategory::sequential;
92  }
93  void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const override
94  {
95  wellMod_.addWellPressureEquations(jacobian, weights, use_well_weights);
96  }
97  void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
98  {
99  wellMod_.addWellPressureEquationsStruct(jacobian);
100  }
101  int getNumberOfExtraEquations() const override
102  {
103  return wellMod_.numLocalWellsEnd();
104  }
105 
106 private:
107  const WellModel& wellMod_;
108 };
109 
121 template<class M, class X, class Y, bool overlapping >
122 class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
123 {
124 public:
125  typedef M matrix_type;
126  typedef X domain_type;
127  typedef Y range_type;
128  typedef typename X::field_type field_type;
129  using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
130 #if HAVE_MPI
131  typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
132 #else
133  typedef Dune::CollectiveCommunication< int > communication_type;
134 #endif
135 
136  Dune::SolverCategory::Category category() const override
137  {
138  return overlapping ?
139  Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
140  }
141 
144  const Opm::LinearOperatorExtra<X, Y>& wellOper,
145  const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >())
146  : A_( A ), wellOper_( wellOper ), comm_(comm)
147  {}
148 
149 
150  virtual void apply( const X& x, Y& y ) const override
151  {
152  A_.mv( x, y );
153 
154  // add well model modification to y
155  wellOper_.apply(x, y );
156 
157 #if HAVE_MPI
158  if( comm_ )
159  comm_->project( y );
160 #endif
161  }
162 
163  // y += \alpha * A * x
164  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
165  {
166  A_.usmv(alpha,x,y);
167 
168  // add scaled well model modification to y
169  wellOper_.applyscaleadd( alpha, x, y );
170 
171 #if HAVE_MPI
172  if( comm_ )
173  comm_->project( y );
174 #endif
175  }
176 
177  virtual const matrix_type& getmat() const override { return A_; }
178 
179  void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
180  {
181  wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
182  }
183  void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
184  {
185  wellOper_.addWellPressureEquationsStruct(jacobian);
186  }
187  int getNumberOfExtraEquations() const
188  {
189  return wellOper_.getNumberOfExtraEquations();
190  }
191 
192 protected:
193  const matrix_type& A_ ;
194  const Opm::LinearOperatorExtra<X, Y>& wellOper_;
195  std::shared_ptr< communication_type > comm_;
196 };
197 
198 
207 template<class M, class X, class Y, bool overlapping >
208 class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
209 {
210 public:
211  typedef M matrix_type;
212  typedef X domain_type;
213  typedef Y range_type;
214  typedef typename X::field_type field_type;
215  using PressureMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<double, 1, 1>>;
216 #if HAVE_MPI
217  typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
218 #else
219  typedef Dune::CollectiveCommunication< int > communication_type;
220 #endif
221 
222 
223  Dune::SolverCategory::Category category() const override
224  {
225  return overlapping ?
226  Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
227  }
228 
231  const Opm::LinearOperatorExtra<X, Y>& wellOper,
232  const size_t interiorSize )
233  : A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize)
234  {}
235 
236  virtual void apply( const X& x, Y& y ) const override
237  {
238  for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
239  {
240  y[row.index()]=0;
241  auto endc = (*row).end();
242  for (auto col = (*row).begin(); col != endc; ++col)
243  (*col).umv(x[col.index()], y[row.index()]);
244  }
245 
246  // add well model modification to y
247  wellOper_.apply(x, y );
248 
249  ghostLastProject( y );
250  }
251 
252  // y += \alpha * A * x
253  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
254  {
255  for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
256  {
257  auto endc = (*row).end();
258  for (auto col = (*row).begin(); col != endc; ++col)
259  (*col).usmv(alpha, x[col.index()], y[row.index()]);
260  }
261  // add scaled well model modification to y
262  wellOper_.applyscaleadd( alpha, x, y );
263 
264  ghostLastProject( y );
265  }
266 
267  virtual const matrix_type& getmat() const override { return A_; }
268 
269  void addWellPressureEquations(PressureMatrix& jacobian, const X& weights,const bool use_well_weights) const
270  {
271  wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
272  }
273  void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
274  {
275  wellOper_.addWellPressureEquationsStruct(jacobian);
276  }
277  int getNumberOfExtraEquations() const
278  {
279  return wellOper_.getNumberOfExtraEquations();
280  }
281 
282 protected:
283  void ghostLastProject(Y& y) const
284  {
285  size_t end = y.size();
286  for (size_t i = interiorSize_; i < end; ++i)
287  y[i] = 0;
288  }
289 
290  const matrix_type& A_ ;
291  const Opm::LinearOperatorExtra< X, Y>& wellOper_;
292  size_t interiorSize_;
293 };
294 
295 } // namespace Opm
296 
297 #endif // OPM_WELLOPERATORS_HEADER_INCLUDED
298 
Linear operator wrapper for well model.
Definition: WellOperators.hpp:50
Definition: WellOperators.hpp:60
void apply(const X &x, Y &y) const override
apply operator to x: The input vector is consistent and the output must also be consistent on the in...
Definition: WellOperators.hpp:73
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const override
apply operator to x, scale and add:
Definition: WellOperators.hpp:79
Dune::SolverCategory::Category category() const override
Category for operator.
Definition: WellOperators.hpp:89
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:209
WellModelGhostLastMatrixAdapter(const M &A, const Opm::LinearOperatorExtra< X, Y > &wellOper, const size_t interiorSize)
constructor: just store a reference to a matrix
Definition: WellOperators.hpp:230
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:123
WellModelMatrixAdapter(const M &A, const Opm::LinearOperatorExtra< X, Y > &wellOper, const std::shared_ptr< communication_type > &comm=std::shared_ptr< communication_type >())
constructor: just store a reference to a matrix
Definition: WellOperators.hpp:143
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27