21 #ifndef OPM_AUTODIFFHELPERS_HEADER_INCLUDED
22 #define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
26 #include <opm/core/grid.h>
27 #include <opm/common/ErrorMacros.hpp>
28 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
29 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
42 typedef Eigen::SparseMatrix<double>
M;
46 typedef Eigen::Array<int, Eigen::Dynamic, 1>
IFaces;
64 typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
TwoColInt;
72 HelperOps(
const Grid& grid, Opm::EclipseStateConstPtr eclState = EclipseStateConstPtr () )
74 using namespace AutoDiffGrid;
75 const int nc = numCells(grid);
76 const int nf = numFaces(grid);
81 int num_internal = internal_faces.size();
83 int num_connections = num_internal;
87 std::shared_ptr<const NNC> nnc = eclState ? eclState->getNNC()
88 : std::shared_ptr<const Opm::NNC>();
89 const bool has_nnc = nnc && nnc->hasNNC();
91 numNNC = nnc->numNNC();
92 num_connections += numNNC;
94 nbi.resize(num_internal, 2);
97 size_t cartesianSize = eclState->getEclipseGrid()->getCartesianSize();
98 std::vector<int> global2localIdx(cartesianSize,0);
99 for (
int i = 0; i< nc; ++i) {
100 global2localIdx[ globalCell( grid )[i] ] = i;
102 const std::vector<size_t>& NNC1 = nnc->nnc1();
103 const std::vector<size_t>& NNC2 = nnc->nnc2();
104 nnc_cells.resize(numNNC,2);
105 for (
int i = 0; i < numNNC; ++i) {
106 nnc_cells(i,0) = global2localIdx[NNC1[i]];
107 nnc_cells(i,1) = global2localIdx[NNC2[i]];
110 nnc_trans = Eigen::Map<const V>(nnc->trans().data(), numNNC);
113 nnc_cells.resize(0, 2);
119 ngrad.resize(num_connections, nc);
120 caver.resize(num_connections, nc);
121 typedef Eigen::Triplet<double> Tri;
122 std::vector<Tri> ngrad_tri;
123 std::vector<Tri> caver_tri;
124 ngrad_tri.reserve(2*num_connections);
125 caver_tri.reserve(2*num_connections);
126 for (
int i = 0; i < num_internal; ++i) {
127 ngrad_tri.emplace_back(i, nbi(i,0), 1.0);
128 ngrad_tri.emplace_back(i, nbi(i,1), -1.0);
129 caver_tri.emplace_back(i, nbi(i,0), 0.5);
130 caver_tri.emplace_back(i, nbi(i,1), 0.5);
134 for (
int i = 0; i < numNNC; ++i) {
135 ngrad_tri.emplace_back(i+num_internal,
nnc_cells(i,0), 1.0);
136 ngrad_tri.emplace_back(i+num_internal,
nnc_cells(i,1), -1.0);
137 caver_tri.emplace_back(i+num_internal,
nnc_cells(i,0), 0.5);
138 caver_tri.emplace_back(i+num_internal,
nnc_cells(i,1), 0.5);
141 ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end());
142 caver.setFromTriplets(caver_tri.begin(), caver_tri.end());
144 div = ngrad.transpose();
146 std::vector<Tri> fullngrad_tri;
147 fullngrad_tri.reserve(2*(nf+numNNC));
149 for (
int i = 0; i < nf; ++i) {
151 fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
154 fullngrad_tri.emplace_back(i, nb(i,1), -1.0);
159 for (
int i = 0; i < numNNC; ++i) {
160 fullngrad_tri.emplace_back(i+nf,
nnc_cells(i,0), 1.0);
161 fullngrad_tri.emplace_back(i+nf,
nnc_cells(i,1), -1.0);
164 fullngrad.resize(nf+numNNC, nc);
165 fullngrad.setFromTriplets(fullngrad_tri.begin(), fullngrad_tri.end());
166 fulldiv = fullngrad.transpose();
175 template <
typename Scalar>
183 const typename ADB::V& ifaceflux)
185 using namespace AutoDiffGrid;
186 typedef HelperOps::IFaces::Index IFIndex;
188 typename ADFaceCellTraits<Grid>::Type
193 int num_connections = nif + num_nnc;
194 assert(num_connections == ifaceflux.size());
197 typedef typename Eigen::Triplet<Scalar> Triplet;
198 std::vector<Triplet> s; s.reserve(num_connections);
199 for (IFIndex iface = 0; iface < nif; ++iface) {
201 const int c1 = face_cells(f,0);
202 const int c2 = face_cells(f,1);
204 assert ((c1 >= 0) && (c2 >= 0));
207 const int c = (ifaceflux[iface] >= 0) ? c1 : c2;
209 s.push_back(Triplet(iface, c, Scalar(1)));
212 for (
int i = 0; i < num_nnc; ++i) {
214 s.push_back(Triplet(i+nif,c,Scalar(1)));
219 select_.resize(num_connections, numCells(g));
220 select_.setFromTriplets(s.begin(), s.end());
230 std::vector<ADB> xf; xf.reserve(xc.size());
231 for (
typename std::vector<ADB>::const_iterator
232 b = xc.begin(), e = xc.end(); b != e; ++b)
234 xf.push_back(select_ * (*b));
249 return (select_*xc.matrix()).array();
253 Eigen::SparseMatrix<double> select_;
261 template <
typename Scalar,
class IntVec>
262 typename Eigen::SparseMatrix<Scalar>
263 constructSupersetSparseMatrix(
const int full_size,
const IntVec& indices)
265 const int subset_size = indices.size();
267 if (subset_size == 0) {
268 Eigen::SparseMatrix<Scalar> mat(full_size, 0);
272 Eigen::SparseMatrix<Scalar> mat(full_size, subset_size);
273 mat.reserve(Eigen::VectorXi::Constant(subset_size, 1));
274 for (
int i = 0; i < subset_size; ++i) {
275 mat.insert(indices[i], i) = 1;
285 template <
typename Scalar,
class IntVec>
286 Eigen::Array<Scalar, Eigen::Dynamic, 1>
287 subset(
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
288 const IntVec& indices)
290 typedef typename Eigen::Array<Scalar, Eigen::Dynamic, 1>::Index Index;
291 const Index size = indices.size();
292 Eigen::Array<Scalar, Eigen::Dynamic, 1> ret( size );
293 for( Index i=0; i<size; ++i )
294 ret[ i ] = x[ indices[ i ] ];
296 return std::move(ret);
300 template <
typename Scalar,
class IntVec>
301 AutoDiffBlock<Scalar>
303 const IntVec& indices)
305 const Eigen::SparseMatrix<Scalar> sub
306 = constructSupersetSparseMatrix<Scalar>(x.
value().size(), indices).transpose().eval();
312 template <
typename Scalar,
class IntVec>
313 AutoDiffBlock<Scalar>
315 const IntVec& indices,
318 return constructSupersetSparseMatrix<Scalar>(n, indices) * x;
324 template <
typename Scalar,
class IntVec>
325 Eigen::Array<Scalar, Eigen::Dynamic, 1>
326 superset(
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
327 const IntVec& indices,
330 return constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
339 Eigen::SparseMatrix<double>
342 typedef Eigen::SparseMatrix<double>
M;
344 const int n = d.
size();
346 mat.reserve(Eigen::ArrayXi::Ones(n, 1));
347 for (M::Index i = 0; i < n; ++i) {
348 mat.insert(i, i) = d[i];
358 template <
typename Scalar>
369 const int n = selection_basis.size();
371 left_elems_.reserve(n);
372 right_elems_.reserve(n);
373 for (
int i = 0; i < n; ++i) {
374 bool chooseleft =
false;
376 case GreaterEqualZero:
377 chooseleft = selection_basis[i] >= 0.0;
380 chooseleft = selection_basis[i] > 0.0;
383 chooseleft = selection_basis[i] == 0.0;
386 chooseleft = selection_basis[i] != 0.0;
389 chooseleft = selection_basis[i] < 0.0;
392 chooseleft = selection_basis[i] <= 0.0;
395 chooseleft = !isnan(selection_basis[i]);
398 OPM_THROW(std::logic_error,
"No such criterion: " << crit);
401 left_elems_.push_back(i);
403 right_elems_.push_back(i);
409 ADB
select(
const ADB& x1,
const ADB& x2)
const
411 if (right_elems_.empty()) {
413 }
else if (left_elems_.empty()) {
424 if (right_elems_.empty()) {
426 }
else if (left_elems_.empty()) {
435 std::vector<int> left_elems_;
436 std::vector<int> right_elems_;
443 template <
class Matrix>
450 typedef Eigen::Triplet<double> Tri;
452 for (
int block = 0; block < nb; ++block) {
457 int block_col_start = 0;
458 for (
int block = 0; block < nb; ++block) {
460 Eigen::SparseMatrix<double> jac;
462 for (Eigen::SparseMatrix<double>::Index k = 0; k < jac.outerSize(); ++k) {
463 for (Eigen::SparseMatrix<double>::InnerIterator i(jac, k); i ; ++i) {
464 if (i.value() != 0.0) {
465 t.push_back(Tri(i.row(),
466 i.col() + block_col_start,
471 block_col_start += jac.cols();
474 jacobian = Matrix(x.
size(), block_col_start);
475 jacobian.setFromTriplets(t.begin(), t.end());
482 AutoDiffBlock<double>
485 Eigen::SparseMatrix<double> comb_jac;
489 std::vector<ADB::M> jacs(1);
499 AutoDiffBlock<double>
503 const int nx = x.
size();
504 const int ny = y.
size();
505 const int n = nx + ny;
506 std::vector<int> xind(nx);
507 for (
int i = 0; i < nx; ++i) {
510 std::vector<int> yind(ny);
511 for (
int i = 0; i < ny; ++i) {
523 AutoDiffBlock<double>
532 const int nx = x.size();
535 int elem_with_deriv = -1;
537 for (
int elem = 0; elem < nx; ++elem) {
538 size += x[elem].size();
539 if (x[elem].derivative().empty()) {
543 if (elem_with_deriv == -1) {
544 elem_with_deriv = elem;
545 num_blocks = x[elem].numBlocks();
548 if (x[elem].blockPattern() != x[elem_with_deriv].blockPattern()) {
549 OPM_THROW(std::runtime_error,
"vertcatCollapseJacs(): all arguments must have the same block pattern");
551 for (
int block = 0; block < num_blocks; ++block) {
552 nnz += x[elem].derivative()[block].nonZeros();
556 for (
int block = 0; block < num_blocks; ++block) {
557 num_cols += x[elem_with_deriv].derivative()[block].cols();
563 for (
int elem = 0; elem < nx; ++elem) {
564 const int loc_size = x[elem].size();
565 val.segment(pos, loc_size) = x[elem].value();
571 if (num_blocks == 0) {
576 typedef Eigen::SparseMatrix<double>
M;
577 typedef Eigen::Triplet<double> Tri;
581 int block_row_start = 0;
583 for (
int elem = 0; elem < nx; ++elem) {
584 int block_col_start = 0;
585 if (!x[elem].derivative().empty()) {
586 for (
int block = 0; block < num_blocks; ++block) {
587 x[elem].derivative()[block].toSparse(jac);
588 for (M::Index k = 0; k < jac.outerSize(); ++k) {
589 for (M::InnerIterator i(jac, k); i ; ++i) {
590 t.push_back(Tri(i.row() + block_row_start,
591 i.col() + block_col_start,
595 block_col_start += jac.
cols();
598 block_row_start += x[elem].size();
603 M comb_jac =
M(size, num_cols);
604 comb_jac.reserve(nnz);
605 comb_jac.setFromTriplets(t.begin(), t.end());
606 std::vector<ADB::M> jac(1);
607 jac[0] =
ADB::M(std::move(comb_jac));
626 Span(
const int num,
const int stride,
const int start)
634 assert(i >= 0 && i < num_);
635 return start_ + i*stride_;
660 return before_increment;
664 assert(span_ == rhs.span_);
665 return index_ < rhs.index_;
669 assert(span_ == rhs.span_);
670 return index_ == rhs.index_;
674 assert(span_ == rhs.span_);
675 return index_ != rhs.index_;
679 return (*span_)[index_];
701 return num_ == rhs.num_ && start_ == rhs.start_ && stride_ == rhs.stride_;
713 inline Eigen::ArrayXd
sign (
const Eigen::ArrayXd& x)
715 const int n = x.size();
716 Eigen::ArrayXd retval(n);
717 for (
int i = 0; i < n; ++i) {
718 retval[i] = x[i] < 0.0 ? -1.0 : (x[i] > 0.0 ? 1.0 : 0.0);
725 #endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:287
SpanIterator begin() const
Definition: AutoDiffHelpers.hpp:689
SpanIterator end() const
Definition: AutoDiffHelpers.hpp:694
Eigen::Array< int, Eigen::Dynamic, 2, Eigen::RowMajor > TwoColInt
Non-neighboring connections.
Definition: AutoDiffHelpers.hpp:64
M caver
Extract for each face the average of its adjacent cells' values.
Definition: AutoDiffHelpers.hpp:54
AutoDiffBlock< double > ADB
Definition: BlackoilModelBase_impl.hpp:83
Eigen::Array< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:98
Definition: AdditionalObjectDeleter.hpp:22
AutoDiffMatrix M
Underlying type for jacobians.
Definition: AutoDiffBlock.hpp:100
AutoDiffBlock< Scalar > ADB
Definition: AutoDiffHelpers.hpp:178
int cols() const
Definition: AutoDiffMatrix.hpp:563
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:103
CriterionForLeftElement
Definition: AutoDiffHelpers.hpp:363
SpanIterator iterator
Definition: AutoDiffHelpers.hpp:686
int size() const
Definition: AutoDiffHelpers.hpp:637
int numBlocks() const
Number of Jacobian blocks.
Definition: AutoDiffBlock.hpp:419
Definition: AutoDiffHelpers.hpp:617
void toSparse(Eigen::SparseMatrix< Scalar, Options, Index > &s) const
Definition: AutoDiffMatrix.hpp:530
ADFaceCellTraits< UnstructuredGrid >::Type faceCellsToEigen(const UnstructuredGrid &grid)
Get the face to cell mapping of a grid.
Definition: AutoDiffHelpers.hpp:363
Selection. Choose first of two elements if selection basis element is nonnegative.
Definition: AutoDiffHelpers.hpp:359
TwoColInt nnc_cells
Definition: AutoDiffHelpers.hpp:65
Eigen::Array< int, Eigen::Dynamic, 1 > IFaces
A list of internal faces.
Definition: AutoDiffHelpers.hpp:46
M grad
Extract for each face the difference of its adjacent cells' values (second - first).
Definition: AutoDiffHelpers.hpp:52
SpanIterator const_iterator
Definition: AutoDiffHelpers.hpp:687
ADB select(const ADB &x1, const ADB &x2) const
Apply selector to ADB quantities.
Definition: AutoDiffHelpers.hpp:409
bool operator==(const SpanIterator &rhs) const
Definition: AutoDiffHelpers.hpp:667
UpwindSelector(const Grid &g, const HelperOps &h, const typename ADB::V &ifaceflux)
Definition: AutoDiffHelpers.hpp:181
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:413
AutoDiffBlock< double >::V V
Definition: AutoDiffHelpers.hpp:43
Span(const int num)
Definition: AutoDiffHelpers.hpp:620
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:314
Selector(const typename ADB::V &selection_basis, CriterionForLeftElement crit=GreaterEqualZero)
Definition: AutoDiffHelpers.hpp:365
bool operator!=(const SpanIterator &rhs) const
Definition: AutoDiffHelpers.hpp:672
HelperOps(const Grid &grid, Opm::EclipseStateConstPtr eclState=EclipseStateConstPtr())
Constructs all helper vectors and matrices.
Definition: AutoDiffHelpers.hpp:72
static AutoDiffBlock constant(V &&val)
Definition: AutoDiffBlock.hpp:110
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:713
int operator*()
Definition: AutoDiffHelpers.hpp:677
Definition: AutoDiffMatrix.hpp:43
bool operator<(const SpanIterator &rhs) const
Definition: AutoDiffHelpers.hpp:662
M ngrad
Extract for each internal face the difference of its adjacent cells' values (first - second)...
Definition: AutoDiffHelpers.hpp:50
SpanIterator operator++()
Definition: AutoDiffHelpers.hpp:651
Eigen::SparseMatrix< double > spdiag(const AutoDiffBlock< double >::V &d)
Definition: AutoDiffHelpers.hpp:340
V nnc_trans
The NNC transmissibilities.
Definition: AutoDiffHelpers.hpp:68
AutoDiffBlock< Scalar > ADB
Definition: AutoDiffHelpers.hpp:361
M div
Extract for each cell the sum of its adjacent interior faces' (signed) values.
Definition: AutoDiffHelpers.hpp:56
std::vector< ADB > select(const std::vector< ADB > &xc) const
Apply selector to multiple per-cell quantities.
Definition: AutoDiffHelpers.hpp:225
Definition: AutoDiffHelpers.hpp:40
int operator[](const int i) const
Definition: AutoDiffHelpers.hpp:632
ADB::V select(const typename ADB::V &x1, const typename ADB::V &x2) const
Apply selector to ADB quantities.
Definition: AutoDiffHelpers.hpp:422
void collapseJacs(const AutoDiffBlock< double > &x, Matrix &jacobian)
Returns the input expression, but with all Jacobians collapsed to one.
Definition: AutoDiffHelpers.hpp:446
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Definition: AutoDiffBlock.hpp:183
IFaces internal_faces
Definition: AutoDiffHelpers.hpp:47
ADB::V select(const typename ADB::V &xc) const
Apply selector to single per-cell constant quantity.
Definition: AutoDiffHelpers.hpp:247
AutoDiffBlock< double > vertcatCollapseJacs(const std::vector< AutoDiffBlock< double > > &x)
Definition: AutoDiffHelpers.hpp:524
ADB select(const ADB &xc) const
Apply selector to single per-cell ADB quantity.
Definition: AutoDiffHelpers.hpp:241
Definition: AutoDiffHelpers.hpp:176
SpanIterator operator++(int)
Definition: AutoDiffHelpers.hpp:656
Definition: AutoDiffHelpers.hpp:643
Eigen::SparseMatrix< double > M
Definition: AutoDiffHelpers.hpp:42
Definition: AutoDiffBlock.hpp:94
Span(const int num, const int stride, const int start)
Definition: AutoDiffHelpers.hpp:626
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:436
bool operator==(const Span &rhs)
Definition: AutoDiffHelpers.hpp:699
M fullngrad
Definition: AutoDiffHelpers.hpp:59
M fulldiv
Extract for each cell the sum of all its adjacent faces' (signed) values.
Definition: AutoDiffHelpers.hpp:61
AutoDiffBlock< double > vertcat(const AutoDiffBlock< double > &x, const AutoDiffBlock< double > &y)
Returns the vertical concatenation [ x; y ] of the inputs.
Definition: AutoDiffHelpers.hpp:500
SpanIterator(const Span *span, const int index)
Definition: AutoDiffHelpers.hpp:646
void extractInternalFaces(const UnstructuredGrid &grid, Eigen::Array< int, Eigen::Dynamic, 1 > &internal_faces, Eigen::Array< int, Eigen::Dynamic, 2, Eigen::RowMajor > &nbi)
extracts the internal faces of a grid.
const std::vector< M > & derivative() const
Function derivatives.
Definition: AutoDiffBlock.hpp:442