22#ifndef OPM_WELLOPERATORS_HEADER_INCLUDED
23#define OPM_WELLOPERATORS_HEADER_INCLUDED
25#include <dune/common/parallel/communication.hh>
26#include <dune/istl/operators.hh>
27#include <dune/istl/bcrsmatrix.hh>
29#include <opm/common/TimingMacros.hpp>
31#include <opm/simulators/linalg/matrixblock.hh>
32#include <dune/common/shared_ptr.hh>
33#include <dune/istl/paamg/smoother.hh>
55template <
class X,
class Y>
59 using field_type =
typename X::field_type;
60 using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
61 virtual void addWellPressureEquations(PressureMatrix& jacobian,
64 virtual void addWellPressureEquationsStruct(PressureMatrix& jacobian)
const = 0;
65 virtual int getNumberOfExtraEquations()
const = 0;
68template <
class WellModel,
class X,
class Y>
73 using field_type =
typename Base::field_type;
74 using PressureMatrix =
typename Base::PressureMatrix;
84 void apply(
const X& x, Y&
y)
const override
87 for (
const auto& well : this->wellMod_) {
88 this->applySingleWell(x,
y, well, well->cells());
96 if (this->wellMod_.empty()) {
100 if (scaleAddRes_.size() !=
y.size()) {
101 scaleAddRes_.resize(
y.size());
106 apply(x, scaleAddRes_);
108 y.axpy(alpha, scaleAddRes_);
116 Dune::SolverCategory::Category
category()
const override
118 return Dune::SolverCategory::sequential;
121 void addWellPressureEquations(PressureMatrix& jacobian,
129 void addWellPressureEquationsStruct(PressureMatrix& jacobian)
const override
132 wellMod_.addWellPressureEquationsStruct(jacobian);
135 int getNumberOfExtraEquations()
const override
137 return wellMod_.numLocalWellsEnd();
141 const WellModel& wellMod_;
143 template<
class WellType,
class ArrayType>
144 void applySingleWell(
const X& x, Y&
y,
145 const WellType& well,
149 x_local_.resize(cells.size());
150 Ax_local_.resize(cells.size());
152 for (
size_t i = 0; i < cells.size(); ++i) {
153 x_local_[i] = x[cells[i]];
154 Ax_local_[i] =
y[cells[i]];
157 well->apply(x_local_, Ax_local_);
159 for (
size_t i = 0; i < cells.size(); ++i) {
161 y[cells[i]] = Ax_local_[i];
169 mutable X x_local_{};
170 mutable Y Ax_local_{};
171 mutable Y scaleAddRes_{};
174template <
class WellModel,
class X,
class Y>
180 using field_type =
typename WBase::field_type;
181 using PressureMatrix =
typename WBase::PressureMatrix;
183 void setDomainIndex(
int index) { domainIndex_ = index; }
185 void apply(
const X& x, Y&
y)
const override
188 std::size_t well_index = 0;
189 for (
const auto& well : this->wellMod_) {
190 if (this->wellMod_.well_domain().at(well->name()) == domainIndex_) {
191 this->applySingleWell(x,
y, well,
192 this->wellMod_.well_local_cells()[well_index]);
198 void addWellPressureEquations(PressureMatrix& jacobian,
203 this->wellMod_.addWellPressureEquationsDomain(jacobian,
210 int domainIndex_ = -1;
223template<
class M,
class X,
class Y>
227 using matrix_type = M;
228 using domain_type = X;
229 using range_type = Y;
230 using field_type =
typename X::field_type;
231 using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
233 Dune::SolverCategory::Category category()
const override
235 return Dune::SolverCategory::sequential;
244 void apply(
const X& x, Y&
y )
const override
250 wellOper_.apply(x,
y);
254 void applyscaleadd (field_type alpha,
const X& x, Y&
y)
const override
257 A_.usmv(alpha, x,
y);
260 wellOper_.applyscaleadd(alpha, x,
y);
263 const matrix_type& getmat()
const override {
return A_; }
265 void addWellPressureEquations(PressureMatrix& jacobian,
273 void addWellPressureEquationsStruct(PressureMatrix& jacobian)
const
276 wellOper_.addWellPressureEquationsStruct(jacobian);
279 int getNumberOfExtraEquations()
const
281 return wellOper_.getNumberOfExtraEquations();
285 const matrix_type& A_ ;
297template<
class M,
class X,
class Y,
bool overlapping >
301 using matrix_type = M;
302 using domain_type = X;
303 using range_type = Y;
304 using field_type =
typename X::field_type;
305 using PressureMatrix = Dune::BCRSMatrix<MatrixBlock<field_type, 1, 1>>;
307 using communication_type = Dune::OwnerOverlapCopyCommunication<int,int>;
309 using communication_type = Dune::Communication<int>;
312 Dune::SolverCategory::Category category()
const override
315 Dune::SolverCategory::overlapping : Dune::SolverCategory::sequential;
325 void apply(
const X& x, Y&
y)
const override
328 for (
auto row = A_.begin();
row.index() < interiorSize_; ++
row)
331 auto endc = (*row).end();
333 (*col).umv(x[
col.index()],
y[
row.index()]);
337 wellOper_.apply(x,
y);
343 void applyscaleadd (field_type alpha,
const X& x, Y&
y)
const override
346 for (
auto row = A_.begin();
row.index() < interiorSize_; ++
row)
348 auto endc = (*row).end();
350 (*col).usmv(alpha, x[
col.index()],
y[
row.index()]);
353 wellOper_.applyscaleadd(alpha, x,
y);
358 const matrix_type& getmat()
const override {
return A_; }
360 void addWellPressureEquations(PressureMatrix& jacobian,
368 void addWellPressureEquationsStruct(PressureMatrix& jacobian)
const
371 wellOper_.addWellPressureEquationsStruct(jacobian);
374 int getNumberOfExtraEquations()
const
376 return wellOper_.getNumberOfExtraEquations();
380 void ghostLastProject(Y&
y)
const
382 std::size_t end =
y.size();
383 for (std::size_t i = interiorSize_; i < end; ++i)
387 const matrix_type& A_ ;
389 std::size_t interiorSize_;
400template<
class M,
class X,
class Y,
class C>
404 typedef M matrix_type;
405 typedef X domain_type;
406 typedef Y range_type;
407 typedef typename X::field_type field_type;
410 typedef C communication_type;
412 Dune::SolverCategory::Category category()
const override
414 return Dune::SolverCategory::overlapping;
419 const communication_type& comm)
422 interiorSize_ = setInteriorSize(comm_);
426 const communication_type& comm)
427 : A_(
A ), comm_(comm)
429 interiorSize_ = setInteriorSize(comm_);
432 virtual void apply(
const X& x, Y&
y )
const override
434 for (
auto row = A_->begin();
row.index() < interiorSize_; ++
row)
437 auto endc = (*row).end();
439 (*col).umv(x[
col.index()],
y[
row.index()]);
442 ghostLastProject(
y );
446 virtual void applyscaleadd (field_type alpha,
const X& x, Y&
y)
const override
448 for (
auto row = A_->begin();
row.index() < interiorSize_; ++
row)
450 auto endc = (*row).end();
452 (*col).usmv(alpha, x[
col.index()],
y[
row.index()]);
455 ghostLastProject(
y );
458 virtual const matrix_type& getmat()
const override {
return *A_; }
460 size_t getInteriorSize()
const {
return interiorSize_;}
463 void ghostLastProject(Y&
y)
const
465 size_t end =
y.size();
466 for (
size_t i = interiorSize_; i < end; ++i)
470 size_t setInteriorSize(
const communication_type& comm)
const
478 if (
idx->local().attribute()==1) {
480 auto loc =
idx->local().local();
489 const std::shared_ptr<const matrix_type> A_ ;
490 const communication_type& comm_;
491 size_t interiorSize_;
499template<
class M,
class X,
class Y,
class C>
500class ConstructionTraits<
Opm::GhostLastMatrixAdapter<M,X,Y,C> >
503 typedef ParallelOperatorArgs<M,C> Arguments;
505 static inline std::shared_ptr<Opm::GhostLastMatrixAdapter<M,X,Y,C>> construct(
const Arguments& args)
507 return std::make_shared<Opm::GhostLastMatrixAdapter<M,X,Y,C>>
508 (args.matrix_, args.comm_);
Definition WellOperators.hpp:176
Dune linear operator that assumes ghost rows are ordered after interior rows.
Definition WellOperators.hpp:402
GhostLastMatrixAdapter(const M &A, const communication_type &comm)
constructor: just store a reference to a matrix
Definition WellOperators.hpp:418
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:93
Dune::SolverCategory::Category category() const override
Category for operator.
Definition WellOperators.hpp:116
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:299
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:319
Adapter to combine a matrix and another linear operator into a combined linear operator.
Definition WellOperators.hpp:225
WellModelMatrixAdapter(const M &A, const LinearOperatorExtra< X, Y > &wellOper)
constructor: just store a reference to a matrix
Definition WellOperators.hpp:239
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242