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 
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 WellModel, class X, class Y>
49 class WellModelAsLinearOperator : public Dune::LinearOperator<X, Y>
50 {
51 public:
52  using Base = Dune::LinearOperator<X, Y>;
53  using field_type = typename Base::field_type;
54 
55  explicit WellModelAsLinearOperator(const WellModel& wm)
56  : wellMod_(wm)
57  {
58  }
63  void apply (const X& x, Y& y) const override
64  {
65  wellMod_.apply(x, y );
66  }
67 
69  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
70  {
71  wellMod_.applyScaleAdd( alpha, x, y );
72  }
73 
79  Dune::SolverCategory::Category category() const override
80  {
81  return Dune::SolverCategory::sequential;
82  }
83 private:
84  const WellModel& wellMod_;
85 };
86 
98 template<class M, class X, class Y, bool overlapping >
99 class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
100 {
101 public:
102  typedef M matrix_type;
103  typedef X domain_type;
104  typedef Y range_type;
105  typedef typename X::field_type field_type;
106 
107 #if HAVE_MPI
108  typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
109 #else
110  typedef Dune::CollectiveCommunication< int > communication_type;
111 #endif
112 
113  Dune::SolverCategory::Category category() const override
114  {
115  return overlapping ?
116  Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
117  }
118 
121  const Dune::LinearOperator<X, Y>& wellOper,
122  const std::shared_ptr< communication_type >& comm = std::shared_ptr< communication_type >())
123  : A_( A ), wellOper_( wellOper ), comm_(comm)
124  {}
125 
126 
127  virtual void apply( const X& x, Y& y ) const override
128  {
129  A_.mv( x, y );
130 
131  // add well model modification to y
132  wellOper_.apply(x, y );
133 
134 #if HAVE_MPI
135  if( comm_ )
136  comm_->project( y );
137 #endif
138  }
139 
140  // y += \alpha * A * x
141  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
142  {
143  A_.usmv(alpha,x,y);
144 
145  // add scaled well model modification to y
146  wellOper_.applyscaleadd( alpha, x, y );
147 
148 #if HAVE_MPI
149  if( comm_ )
150  comm_->project( y );
151 #endif
152  }
153 
154  virtual const matrix_type& getmat() const override { return A_; }
155 
156 protected:
157  const matrix_type& A_ ;
158  const Dune::LinearOperator<X, Y>& wellOper_;
159  std::shared_ptr< communication_type > comm_;
160 };
161 
162 
171 template<class M, class X, class Y, bool overlapping >
172 class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
173 {
174 public:
175  typedef M matrix_type;
176  typedef X domain_type;
177  typedef Y range_type;
178  typedef typename X::field_type field_type;
179 
180 #if HAVE_MPI
181  typedef Dune::OwnerOverlapCopyCommunication<int,int> communication_type;
182 #else
183  typedef Dune::CollectiveCommunication< int > communication_type;
184 #endif
185 
186 
187  Dune::SolverCategory::Category category() const override
188  {
189  return overlapping ?
190  Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
191  }
192 
195  const Dune::LinearOperator<X, Y>& wellOper,
196  const size_t interiorSize )
197  : A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize)
198  {}
199 
200  virtual void apply( const X& x, Y& y ) const override
201  {
202  for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
203  {
204  y[row.index()]=0;
205  auto endc = (*row).end();
206  for (auto col = (*row).begin(); col != endc; ++col)
207  (*col).umv(x[col.index()], y[row.index()]);
208  }
209 
210  // add well model modification to y
211  wellOper_.apply(x, y );
212 
213  ghostLastProject( y );
214  }
215 
216  // y += \alpha * A * x
217  virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
218  {
219  for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
220  {
221  auto endc = (*row).end();
222  for (auto col = (*row).begin(); col != endc; ++col)
223  (*col).usmv(alpha, x[col.index()], y[row.index()]);
224  }
225  // add scaled well model modification to y
226  wellOper_.applyscaleadd( alpha, x, y );
227 
228  ghostLastProject( y );
229  }
230 
231  virtual const matrix_type& getmat() const override { return A_; }
232 
233 protected:
234  void ghostLastProject(Y& y) const
235  {
236  size_t end = y.size();
237  for (size_t i = interiorSize_; i < end; ++i)
238  y[i] = 0;
239  }
240 
241  const matrix_type& A_ ;
242  const Dune::LinearOperator<X, Y>& wellOper_;
243  size_t interiorSize_;
244 };
245 
246 } // namespace Opm
247 
248 #endif // OPM_WELLOPERATORS_HEADER_INCLUDED
249 
Linear operator wrapper for well model.
Definition: WellOperators.hpp:50
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:63
virtual void applyscaleadd(field_type alpha, const X &x, Y &y) const override
apply operator to x, scale and add:
Definition: WellOperators.hpp:69
Dune::SolverCategory::Category category() const override
Category for operator.
Definition: WellOperators.hpp:79
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:173
WellModelGhostLastMatrixAdapter(const M &A, const Dune::LinearOperator< X, Y > &wellOper, const size_t interiorSize)
constructor: just store a reference to a matrix
Definition: WellOperators.hpp:194
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition: WellOperators.hpp:100
WellModelMatrixAdapter(const M &A, const Dune::LinearOperator< 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:120
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27