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/common/parallel/communication.hh>
26#include <dune/istl/operators.hh>
27#include <dune/istl/bcrsmatrix.hh>
28
29#include <opm/common/TimingMacros.hpp>
30
31#include <opm/simulators/linalg/matrixblock.hh>
32#include <dune/common/shared_ptr.hh>
33#include <dune/istl/paamg/smoother.hh>
34
35#include <cstddef>
36
37namespace Opm {
38
39//=====================================================================
40// Implementation for ISTL-matrix based operators
41// Note: the classes WellModelMatrixAdapter and
42// WellModelGhostLastMatrixAdapter were moved from ISTLSolver.hpp
43// and subsequently modified.
44//=====================================================================
45
55template <class X, class Y>
56class LinearOperatorExtra : public Dune::LinearOperator<X, Y>
57{
58public:
59 using field_type = typename X::field_type;
60 using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
61 virtual void addWellPressureEquations(PressureMatrix& jacobian,
62 const X& weights,
63 const bool use_well_weights) const = 0;
64 virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian) const = 0;
65 virtual int getNumberOfExtraEquations() const = 0;
66};
67
68template <class WellModel, class X, class Y>
70{
71public:
73 using field_type = typename Base::field_type;
74 using PressureMatrix = typename Base::PressureMatrix;
75 explicit WellModelAsLinearOperator(const WellModel& wm)
76 : wellMod_(wm)
77 {
78 }
79
84 void apply(const X& x, Y& y) const override
85 {
86 OPM_TIMEBLOCK(apply);
87 wellMod_.apply(x, y);
88 }
89
91 void applyscaleadd(field_type alpha, const X& x, Y& y) const override
92 {
93 OPM_TIMEBLOCK(applyscaleadd);
94 wellMod_.applyScaleAdd(alpha, x, y);
95 }
96
102 Dune::SolverCategory::Category category() const override
103 {
104 return Dune::SolverCategory::sequential;
105 }
106
107 void addWellPressureEquations(PressureMatrix& jacobian,
108 const X& weights,
109 const bool use_well_weights) const override
110 {
111 OPM_TIMEBLOCK(addWellPressureEquations);
112 wellMod_.addWellPressureEquations(jacobian, weights, use_well_weights);
113 }
114
115 void addWellPressureEquationsStruct(PressureMatrix& jacobian) const override
116 {
117 OPM_TIMEBLOCK(addWellPressureEquationsStruct);
118 wellMod_.addWellPressureEquationsStruct(jacobian);
119 }
120
121 int getNumberOfExtraEquations() const override
122 {
123 return wellMod_.numLocalWellsEnd();
124 }
125
126protected:
127 const WellModel& wellMod_;
128};
129
130
131template <class WellModel, class X, class Y>
133{
134public:
136 using WBase::WBase; // inherit all constructors from the base class
137 using field_type = typename WBase::field_type;
138 using PressureMatrix = typename WBase::PressureMatrix;
139
140 void setDomainIndex(int index) { domainIndex_ = index; }
141
142 void apply(const X& x, Y& y) const override
143 {
144 OPM_TIMEBLOCK(apply);
145 this->wellMod_.applyDomain(x, y, domainIndex_);
146 }
147
148 void applyscaleadd(field_type alpha, const X& x, Y& y) const override
149 {
150 OPM_TIMEBLOCK(applyscaleadd);
151 this->wellMod_.applyScaleAddDomain(alpha, x, y, domainIndex_);
152 }
153
154 void addWellPressureEquations(PressureMatrix& jacobian,
155 const X& weights,
156 const bool use_well_weights) const override
157 {
158 OPM_TIMEBLOCK(addWellPressureEquations);
159 this->wellMod_.addWellPressureEquationsDomain(jacobian, weights, use_well_weights, domainIndex_);
160 }
161
162private:
163 int domainIndex_ = -1;
164};
165
177template<class M, class X, class Y, bool overlapping >
178class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
179{
180public:
181 using matrix_type = M;
182 using domain_type = X;
183 using range_type = Y;
184 using field_type = typename X::field_type;
185 using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
186#if HAVE_MPI
187 using communication_type = Dune::OwnerOverlapCopyCommunication<int,int>;
188#else
189 using communication_type = Dune::Communication<int>;
190#endif
191
192 Dune::SolverCategory::Category category() const override
193 {
194 return overlapping ?
195 Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
196 }
197
200 const LinearOperatorExtra<X, Y>& wellOper,
201 const std::shared_ptr<communication_type>& comm = {})
202 : A_( A ), wellOper_( wellOper ), comm_(comm)
203 {}
204
205 void apply( const X& x, Y& y ) const override
206 {
207 OPM_TIMEBLOCK(apply);
208 A_.mv(x, y);
209
210 // add well model modification to y
211 wellOper_.apply(x, y);
212
213 #if HAVE_MPI
214 if (comm_) {
215 comm_->project(y);
216 }
217 #endif
218 }
219
220 // y += \alpha * A * x
221 void applyscaleadd (field_type alpha, const X& x, Y& y) const override
222 {
223 OPM_TIMEBLOCK(applyscaleadd);
224 A_.usmv(alpha, x, y);
225
226 // add scaled well model modification to y
227 wellOper_.applyscaleadd(alpha, x, y);
228
229 #if HAVE_MPI
230 if (comm_) {
231 comm_->project( y );
232 }
233 #endif
234 }
235
236 const matrix_type& getmat() const override { return A_; }
237
238 void addWellPressureEquations(PressureMatrix& jacobian,
239 const X& weights,
240 const bool use_well_weights) const
241 {
242 OPM_TIMEBLOCK(addWellPressureEquations);
243 wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
244 }
245
246 void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
247 {
248 OPM_TIMEBLOCK(addWellPressureEquations);
249 wellOper_.addWellPressureEquationsStruct(jacobian);
250 }
251
252 int getNumberOfExtraEquations() const
253 {
254 return wellOper_.getNumberOfExtraEquations();
255 }
256
257protected:
258 const matrix_type& A_ ;
259 const LinearOperatorExtra<X, Y>& wellOper_;
260 std::shared_ptr<communication_type> comm_;
261};
262
271template<class M, class X, class Y, bool overlapping >
272class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
273{
274public:
275 using matrix_type = M;
276 using domain_type = X;
277 using range_type = Y;
278 using field_type = typename X::field_type;
279 using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
280#if HAVE_MPI
281 using communication_type = Dune::OwnerOverlapCopyCommunication<int,int>;
282#else
283 using communication_type = Dune::Communication<int>;
284#endif
285
286 Dune::SolverCategory::Category category() const override
287 {
288 return overlapping ?
289 Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
290 }
291
294 const LinearOperatorExtra<X, Y>& wellOper,
295 const std::size_t interiorSize )
296 : A_( A ), wellOper_( wellOper ), interiorSize_(interiorSize)
297 {}
298
299 void apply(const X& x, Y& y) const override
300 {
301 OPM_TIMEBLOCK(apply);
302 for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
303 {
304 y[row.index()]=0;
305 auto endc = (*row).end();
306 for (auto col = (*row).begin(); col != endc; ++col)
307 (*col).umv(x[col.index()], y[row.index()]);
308 }
309
310 // add well model modification to y
311 wellOper_.apply(x, y);
312
313 ghostLastProject(y);
314 }
315
316 // y += \alpha * A * x
317 void applyscaleadd (field_type alpha, const X& x, Y& y) const override
318 {
319 OPM_TIMEBLOCK(applyscaleadd);
320 for (auto row = A_.begin(); row.index() < interiorSize_; ++row)
321 {
322 auto endc = (*row).end();
323 for (auto col = (*row).begin(); col != endc; ++col)
324 (*col).usmv(alpha, x[col.index()], y[row.index()]);
325 }
326 // add scaled well model modification to y
327 wellOper_.applyscaleadd(alpha, x, y);
328
329 ghostLastProject(y);
330 }
331
332 const matrix_type& getmat() const override { return A_; }
333
334 void addWellPressureEquations(PressureMatrix& jacobian,
335 const X& weights,
336 const bool use_well_weights) const
337 {
338 OPM_TIMEBLOCK(addWellPressureEquations);
339 wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights);
340 }
341
342 void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
343 {
344 OPM_TIMEBLOCK(addWellPressureEquationsStruct);
345 wellOper_.addWellPressureEquationsStruct(jacobian);
346 }
347
348 int getNumberOfExtraEquations() const
349 {
350 return wellOper_.getNumberOfExtraEquations();
351 }
352
353protected:
354 void ghostLastProject(Y& y) const
355 {
356 std::size_t end = y.size();
357 for (std::size_t i = interiorSize_; i < end; ++i)
358 y[i] = 0;
359 }
360
361 const matrix_type& A_ ;
362 const LinearOperatorExtra<X, Y>& wellOper_;
363 std::size_t interiorSize_;
364};
365
374template<class M, class X, class Y, class C>
375class GhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
376{
377public:
378 typedef M matrix_type;
379 typedef X domain_type;
380 typedef Y range_type;
381 typedef typename X::field_type field_type;
382
383
384 typedef C communication_type;
385
386 Dune::SolverCategory::Category category() const override
387 {
388 return Dune::SolverCategory::overlapping;
389 }
390
393 const communication_type& comm)
394 : A_( Dune::stackobject_to_shared_ptr(A) ), comm_(comm)
395 {
396 interiorSize_ = setInteriorSize(comm_);
397 }
398
399 GhostLastMatrixAdapter (const std::shared_ptr<M> A,
400 const communication_type& comm)
401 : A_( A ), comm_(comm)
402 {
403 interiorSize_ = setInteriorSize(comm_);
404 }
405
406 virtual void apply( const X& x, Y& y ) const override
407 {
408 for (auto row = A_->begin(); row.index() < interiorSize_; ++row)
409 {
410 y[row.index()]=0;
411 auto endc = (*row).end();
412 for (auto col = (*row).begin(); col != endc; ++col)
413 (*col).umv(x[col.index()], y[row.index()]);
414 }
415
416 ghostLastProject( y );
417 }
418
419 // y += \alpha * A * x
420 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override
421 {
422 for (auto row = A_->begin(); row.index() < interiorSize_; ++row)
423 {
424 auto endc = (*row).end();
425 for (auto col = (*row).begin(); col != endc; ++col)
426 (*col).usmv(alpha, x[col.index()], y[row.index()]);
427 }
428
429 ghostLastProject( y );
430 }
431
432 virtual const matrix_type& getmat() const override { return *A_; }
433
434 size_t getInteriorSize() const { return interiorSize_;}
435
436private:
437 void ghostLastProject(Y& y) const
438 {
439 size_t end = y.size();
440 for (size_t i = interiorSize_; i < end; ++i)
441 y[i] = 0; //project to interiorsize, i.e. ignore ghost
442 }
443
444 size_t setInteriorSize(const communication_type& comm) const
445 {
446 auto indexSet = comm.indexSet();
447
448 size_t is = 0;
449 // Loop over index set
450 for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) {
451 //Only take "owner" indices
452 if (idx->local().attribute()==1) {
453 //get local index
454 auto loc = idx->local().local();
455 // if loc is higher than "old interior size", update it
456 if (loc > is) {
457 is = loc;
458 }
459 }
460 }
461 return is + 1; //size is plus 1 since we start at 0
462 }
463 const std::shared_ptr<const matrix_type> A_ ;
464 const communication_type& comm_;
465 size_t interiorSize_;
466};
467
468} // namespace Opm
469
470namespace Dune {
471namespace Amg {
472
473template<class M, class X, class Y, class C>
474class ConstructionTraits<Opm::GhostLastMatrixAdapter<M,X,Y,C> >
475{
476public:
477 typedef ParallelOperatorArgs<M,C> Arguments;
478
479 static inline std::shared_ptr<Opm::GhostLastMatrixAdapter<M,X,Y,C>> construct(const Arguments& args)
480 {
481 return std::make_shared<Opm::GhostLastMatrixAdapter<M,X,Y,C>>
482 (args.matrix_, args.comm_);
483 }
484};
485
486} // end namespace Amg
487} // end namespace Dune
488
489#endif // OPM_WELLOPERATORS_HEADER_INCLUDED
490
Definition WellOperators.hpp:133
Dune linear operator that assumes ghost rows are ordered after interior rows.
Definition WellOperators.hpp:376
GhostLastMatrixAdapter(const M &A, const communication_type &comm)
constructor: just store a reference to a matrix
Definition WellOperators.hpp:392
Linear operator wrapper for well model.
Definition WellOperators.hpp:57
Definition WellOperators.hpp:70
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:84
void applyscaleadd(field_type alpha, const X &x, Y &y) const override
apply operator to x, scale and add:
Definition WellOperators.hpp:91
Dune::SolverCategory::Category category() const override
Category for operator.
Definition WellOperators.hpp:102
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:273
WellModelGhostLastMatrixAdapter(const M &A, const LinearOperatorExtra< X, Y > &wellOper, const std::size_t interiorSize)
constructor: just store a reference to a matrix
Definition WellOperators.hpp:293
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:179
WellModelMatrixAdapter(const M &A, const LinearOperatorExtra< X, Y > &wellOper, const std::shared_ptr< communication_type > &comm={})
constructor: just store a reference to a matrix
Definition WellOperators.hpp:199
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37