My Project
Loading...
Searching...
No Matches
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
28namespace 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
48template <class X, class Y>
49class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
50{
51public:
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
58template <class WellModel, class X, class Y>
60{
61public:
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
106private:
107 const WellModel& wellMod_;
108};
109
121template<class M, class X, class Y, bool overlapping >
122class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
123{
124public:
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
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
192protected:
193 const matrix_type& A_ ;
194 const Opm::LinearOperatorExtra<X, Y>& wellOper_;
195 std::shared_ptr< communication_type > comm_;
196};
197
198
207template<class M, class X, class Y, bool overlapping >
208class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
209{
210public:
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
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
282protected:
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
Definition AquiferInterface.hpp:35
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