21 #ifndef OPM_AUTODIFFHELPERS_HEADER_INCLUDED
22 #define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
24 #include <opm/autodiff/AutoDiffBlock.hpp>
25 #include <opm/autodiff/GridHelpers.hpp>
26 #include <opm/autodiff/GeoProps.hpp>
27 #include <opm/core/grid.h>
28 #include <opm/core/grid/PinchProcessor.hpp>
29 #include <opm/core/props/rock/RockFromDeck.hpp>
31 #include <opm/common/ErrorMacros.hpp>
32 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
33 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
46 typedef Eigen::SparseMatrix<double> M;
50 typedef Eigen::Array<int, Eigen::Dynamic, 1>
IFaces;
68 typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
TwoColInt;
79 HelperOps(
const Grid& grid,
const NNC& nnc = NNC())
81 using namespace AutoDiffGrid;
82 const int nc = numCells(grid);
83 const int nf = numFaces(grid);
87 extractInternalFaces(grid, internal_faces, nbi);
88 const int num_internal = internal_faces.size();
91 const bool has_nnc = nnc.hasNNC();
92 int numNNC = nnc.numNNC();
95 const int num_connections = num_internal + numNNC;
97 const int *cartDims = AutoDiffGrid::cartDims(grid);
98 const int numCartesianCells = cartDims[0] * cartDims[1] * cartDims[2];
101 std::vector<int> global2localIdx(numCartesianCells,0);
102 for (
int i = 0; i<
nc; ++i) {
103 global2localIdx[ globalCell( grid )[i] ] = i;
105 nnc_cells.resize(numNNC,2);
107 for (
int i = 0; i < numNNC; ++i) {
108 nnc_cells(i,0) = global2localIdx[nnc.nncdata()[i].cell1];
109 nnc_cells(i,1) = global2localIdx[nnc.nncdata()[i].cell2];
118 ngrad.resize(num_connections, nc);
119 caver.resize(num_connections, nc);
120 typedef Eigen::Triplet<double> Tri;
121 std::vector<Tri> ngrad_tri;
122 std::vector<Tri> caver_tri;
123 ngrad_tri.reserve(2*num_connections);
124 caver_tri.reserve(2*num_connections);
125 for (
int i = 0; i < num_internal; ++i) {
126 ngrad_tri.emplace_back(i, nbi(i,0), 1.0);
127 ngrad_tri.emplace_back(i, nbi(i,1), -1.0);
128 caver_tri.emplace_back(i, nbi(i,0), 0.5);
129 caver_tri.emplace_back(i, nbi(i,1), 0.5);
133 for (
int i = 0; i < numNNC; ++i) {
134 ngrad_tri.emplace_back(i+num_internal, nnc_cells(i,0), 1.0);
135 ngrad_tri.emplace_back(i+num_internal, nnc_cells(i,1), -1.0);
136 caver_tri.emplace_back(i+num_internal, nnc_cells(i,0), 0.5);
137 caver_tri.emplace_back(i+num_internal, nnc_cells(i,1), 0.5);
140 ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end());
141 caver.setFromTriplets(caver_tri.begin(), caver_tri.end());
145 std::vector<Tri> fullngrad_tri;
146 fullngrad_tri.reserve(2*(nf+numNNC));
147 typename ADFaceCellTraits<Grid>::Type nb = faceCellsToEigen(grid);
148 for (
int i = 0; i < nf; ++i) {
150 fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
153 fullngrad_tri.emplace_back(i, nb(i,1), -1.0);
158 for (
int i = 0; i < numNNC; ++i) {
159 fullngrad_tri.emplace_back(i+nf, nnc_cells(i,0), 1.0);
160 fullngrad_tri.emplace_back(i+nf, nnc_cells(i,1), -1.0);
164 fullngrad.setFromTriplets(fullngrad_tri.begin(), fullngrad_tri.end());
180 template <
typename Scalar>
188 const typename ADB::V& ifaceflux)
190 using namespace AutoDiffGrid;
191 typedef HelperOps::IFaces::Index IFIndex;
192 const IFIndex nif = h.internal_faces.size();
193 typename ADFaceCellTraits<Grid>::Type
194 face_cells = faceCellsToEigen(g);
198 int num_connections = nif + num_nnc;
199 assert(num_connections == ifaceflux.size());
202 typedef typename Eigen::Triplet<Scalar> Triplet;
203 std::vector<Triplet> s; s.reserve(num_connections);
204 for (IFIndex iface = 0; iface < nif; ++iface) {
205 const int f = h.internal_faces[iface];
206 const int c1 = face_cells(f,0);
207 const int c2 = face_cells(f,1);
209 assert ((c1 >= 0) && (c2 >= 0));
212 const int c = (ifaceflux[iface] >= 0) ? c1 : c2;
214 s.push_back(Triplet(iface, c, Scalar(1)));
217 for (
int i = 0; i < num_nnc; ++i) {
218 const int c = (ifaceflux[i+nif] >= 0) ? h.nnc_cells(i,0) : h.nnc_cells(i,1);
219 s.push_back(Triplet(i+nif,c,Scalar(1)));
224 select_.resize(num_connections, numCells(g));
225 select_.setFromTriplets(s.begin(), s.end());
235 std::vector<ADB> xf; xf.reserve(xc.size());
236 for (
typename std::vector<ADB>::const_iterator
237 b = xc.begin(), e = xc.end(); b != e; ++b)
239 xf.push_back(select_ * (*b));
254 return (select_*xc.matrix()).array();
258 Eigen::SparseMatrix<double> select_;
266 template <
typename Scalar,
class IntVec>
267 typename Eigen::SparseMatrix<Scalar>
268 constructSupersetSparseMatrix(
const int full_size,
const IntVec& indices)
270 const int subset_size = indices.size();
272 if (subset_size == 0) {
273 Eigen::SparseMatrix<Scalar> mat(full_size, 0);
277 Eigen::SparseMatrix<Scalar> mat(full_size, subset_size);
278 mat.reserve(Eigen::VectorXi::Constant(subset_size, 1));
279 for (
int i = 0; i < subset_size; ++i) {
280 mat.insert(indices[i], i) = 1;
290 template <
typename Scalar,
class IntVec>
291 Eigen::Array<Scalar, Eigen::Dynamic, 1>
292 subset(
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
293 const IntVec& indices)
295 typedef typename Eigen::Array<Scalar, Eigen::Dynamic, 1>::Index Index;
296 const Index size = indices.size();
297 Eigen::Array<Scalar, Eigen::Dynamic, 1> ret( size );
298 for( Index i=0; i<size; ++i )
299 ret[ i ] = x[ indices[ i ] ];
301 return std::move(ret);
305 template <
typename Scalar,
class IntVec>
306 AutoDiffBlock<Scalar>
308 const IntVec& indices)
310 const Eigen::SparseMatrix<Scalar> sub
311 = constructSupersetSparseMatrix<Scalar>(x.
value().size(), indices).transpose().eval();
317 template <
typename Scalar,
class IntVec>
318 AutoDiffBlock<Scalar>
320 const IntVec& indices,
323 return constructSupersetSparseMatrix<Scalar>(n, indices) * x;
329 template <
typename Scalar,
class IntVec>
330 Eigen::Array<Scalar, Eigen::Dynamic, 1>
331 superset(
const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
332 const IntVec& indices,
335 return constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
344 Eigen::SparseMatrix<double>
347 typedef Eigen::SparseMatrix<double>
M;
349 const int n = d.
size();
351 mat.reserve(Eigen::ArrayXi::Ones(n, 1));
352 for (M::Index i = 0; i < n; ++i) {
353 mat.insert(i, i) = d[i];
363 template <
typename Scalar>
368 enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN };
371 CriterionForLeftElement crit = GreaterEqualZero)
375 const int n = selection_basis.size();
377 left_elems_.reserve(n);
378 right_elems_.reserve(n);
379 for (
int i = 0; i < n; ++i) {
380 bool chooseleft =
false;
382 case GreaterEqualZero:
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 = selection_basis[i] < 0.0;
398 chooseleft = selection_basis[i] <= 0.0;
401 chooseleft = !isnan(selection_basis[i]);
404 OPM_THROW(std::logic_error,
"No such criterion: " << crit);
407 left_elems_.push_back(i);
409 right_elems_.push_back(i);
417 if (right_elems_.empty()) {
419 }
else if (left_elems_.empty()) {
430 if (right_elems_.empty()) {
432 }
else if (left_elems_.empty()) {
441 std::vector<int> left_elems_;
442 std::vector<int> right_elems_;
449 template <
class Matrix>
456 typedef Eigen::Triplet<double> Tri;
458 for (
int block = 0; block < nb; ++block) {
463 int block_col_start = 0;
464 for (
int block = 0; block < nb; ++block) {
466 Eigen::SparseMatrix<double> jac;
468 for (Eigen::SparseMatrix<double>::Index k = 0; k < jac.outerSize(); ++k) {
469 for (Eigen::SparseMatrix<double>::InnerIterator i(jac, k); i ; ++i) {
470 if (i.value() != 0.0) {
471 t.push_back(Tri(i.row(),
472 i.col() + block_col_start,
477 block_col_start += jac.cols();
480 jacobian = Matrix(x.
size(), block_col_start);
481 jacobian.setFromTriplets(t.begin(), t.end());
488 AutoDiffBlock<double>
491 Eigen::SparseMatrix<double> comb_jac;
495 std::vector<ADB::M> jacs(1);
505 AutoDiffBlock<double>
509 const int nx = x.
size();
510 const int ny = y.
size();
511 const int n = nx + ny;
512 std::vector<int> xind(nx);
513 for (
int i = 0; i < nx; ++i) {
516 std::vector<int> yind(ny);
517 for (
int i = 0; i < ny; ++i) {
529 AutoDiffBlock<double>
538 const int nx = x.size();
541 int elem_with_deriv = -1;
543 for (
int elem = 0; elem < nx; ++elem) {
544 size += x[elem].size();
545 if (x[elem].derivative().empty()) {
549 if (elem_with_deriv == -1) {
550 elem_with_deriv = elem;
551 num_blocks = x[elem].numBlocks();
554 if (x[elem].blockPattern() != x[elem_with_deriv].blockPattern()) {
555 OPM_THROW(std::runtime_error,
"vertcatCollapseJacs(): all arguments must have the same block pattern");
557 for (
int block = 0; block < num_blocks; ++block) {
558 nnz += x[elem].derivative()[block].nonZeros();
562 for (
int block = 0; block < num_blocks; ++block) {
563 num_cols += x[elem_with_deriv].derivative()[block].cols();
569 for (
int elem = 0; elem < nx; ++elem) {
570 const int loc_size = x[elem].size();
571 val.segment(pos, loc_size) = x[elem].value();
577 if (num_blocks == 0) {
582 typedef Eigen::SparseMatrix<double>
M;
583 typedef Eigen::Triplet<double> Tri;
587 int block_row_start = 0;
589 for (
int elem = 0; elem < nx; ++elem) {
590 int block_col_start = 0;
591 if (!x[elem].derivative().empty()) {
592 for (
int block = 0; block < num_blocks; ++block) {
593 x[elem].derivative()[block].
toSparse(jac);
594 for (M::Index k = 0; k < jac.outerSize(); ++k) {
595 for (M::InnerIterator i(jac, k); i ; ++i) {
596 t.push_back(Tri(i.row() + block_row_start,
597 i.col() + block_col_start,
601 block_col_start += jac.
cols();
604 block_row_start += x[elem].size();
609 M comb_jac = M(size, num_cols);
610 comb_jac.reserve(nnz);
611 comb_jac.setFromTriplets(t.begin(), t.end());
612 std::vector<ADB::M> jac(1);
613 jac[0] =
ADB::M(std::move(comb_jac));
626 explicit Span(
const int num)
632 Span(
const int num,
const int stride,
const int start)
638 int operator[](
const int i)
const
640 assert(i >= 0 && i < num_);
641 return start_ + i*stride_;
666 return before_increment;
670 assert(span_ == rhs.span_);
671 return index_ < rhs.index_;
675 assert(span_ == rhs.span_);
676 return index_ == rhs.index_;
680 assert(span_ == rhs.span_);
681 return index_ != rhs.index_;
685 return (*span_)[index_];
700 SpanIterator end()
const
702 return SpanIterator(
this, num_);
705 bool operator==(
const Span& rhs)
707 return num_ == rhs.num_ && start_ == rhs.start_ && stride_ == rhs.stride_;
719 inline Eigen::ArrayXd
sign (
const Eigen::ArrayXd& x)
721 const int n = x.size();
722 Eigen::ArrayXd retval(n);
723 for (
int i = 0; i < n; ++i) {
724 retval[i] = x[i] < 0.0 ? -1.0 : (x[i] > 0.0 ? 1.0 : 0.0);
731 #endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED
Contains vectors and sparse matrices that represent subsets or operations on (AD or regular) vectors ...
Definition: AutoDiffHelpers.hpp:44
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
const std::vector< M > & derivative() const
Function derivatives.
Definition: AutoDiffBlock.hpp:444
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
V nnc_trans
The NNC transmissibilities.
Definition: AutoDiffHelpers.hpp:72
M div
Extract for each cell the sum of its adjacent interior faces' (signed) values.
Definition: AutoDiffHelpers.hpp:60
const int nc
Constructs all helper vectors and matrices.
Definition: AutoDiffHelpers.hpp:82
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
AutoDiffMatrix M
Underlying type for jacobians.
Definition: AutoDiffBlock.hpp:101
ADB select(const ADB &xc) const
Apply selector to single per-cell ADB quantity.
Definition: AutoDiffHelpers.hpp:246
Eigen::SparseMatrix< double > spdiag(const AutoDiffBlock< double >::V &d)
Construct square sparse matrix with the elements of d on the diagonal.
Definition: AutoDiffHelpers.hpp:345
void collapseJacs(const AutoDiffBlock< double > &x, Matrix &jacobian)
Returns the input expression, but with all Jacobians collapsed to one.
Definition: AutoDiffHelpers.hpp:452
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
Eigen::Array< int, Eigen::Dynamic, 2, Eigen::RowMajor > TwoColInt
Non-neighboring connections.
Definition: AutoDiffHelpers.hpp:68
void toSparse(Eigen::SparseMatrix< Scalar, Options, Index > &s) const
Converts the AutoDiffMatrix to an Eigen SparseMatrix.This might be an expensive operation to perform ...
Definition: AutoDiffMatrix.hpp:542
A class for forward-mode automatic differentiation with vector values and sparse jacobian matrices...
Definition: AutoDiffBlock.hpp:95
Definition: AutoDiffHelpers.hpp:649
int numBlocks() const
Number of Jacobian blocks.
Definition: AutoDiffBlock.hpp:421
Selection. Choose first of two elements if selection basis element is nonnegative.
Definition: AutoDiffHelpers.hpp:364
TwoColInt connection_cells
The set of all connections' cells (face or nnc).
Definition: AutoDiffHelpers.hpp:75
M ngrad
Extract for each internal face the difference of its adjacent cells' values (first - second)...
Definition: AutoDiffHelpers.hpp:54
M caver
Extract for each face the average of its adjacent cells' values.
Definition: AutoDiffHelpers.hpp:58
int cols() const
Returns number of columns in the matrix.
Definition: AutoDiffMatrix.hpp:575
std::vector< ADB > select(const std::vector< ADB > &xc) const
Apply selector to multiple per-cell quantities.
Definition: AutoDiffHelpers.hpp:230
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:319
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
ADB select(const ADB &x1, const ADB &x2) const
Apply selector to ADB quantities.
Definition: AutoDiffHelpers.hpp:415
ADB::V select(const typename ADB::V &x1, const typename ADB::V &x2) const
Apply selector to ADB quantities.
Definition: AutoDiffHelpers.hpp:428
Definition: AutoDiffHelpers.hpp:623
AutoDiffBlock< double > vertcatCollapseJacs(const std::vector< AutoDiffBlock< double > > &x)
Returns the vertical concatenation [ x[0]; x[1]; ...; x[n-1] ] of the inputs.
Definition: AutoDiffHelpers.hpp:530
Upwind selection in absence of counter-current flow (i.e., without effects of gravity and/or capillar...
Definition: AutoDiffHelpers.hpp:181
Eigen::Array< int, Eigen::Dynamic, 1 > IFaces
A list of internal faces.
Definition: AutoDiffHelpers.hpp:50
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
M grad
Extract for each face the difference of its adjacent cells' values (second - first).
Definition: AutoDiffHelpers.hpp:56
M fulldiv
Extract for each cell the sum of all its adjacent faces' (signed) values.
Definition: AutoDiffHelpers.hpp:65
Eigen::Array< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
M fullngrad
Extract for each face the difference of its adjacent cells' values (first - second).
Definition: AutoDiffHelpers.hpp:63
AutoDiffBlock< double > vertcat(const AutoDiffBlock< double > &x, const AutoDiffBlock< double > &y)
Returns the vertical concatenation [ x; y ] of the inputs.
Definition: AutoDiffHelpers.hpp:506
ADB::V select(const typename ADB::V &xc) const
Apply selector to single per-cell constant quantity.
Definition: AutoDiffHelpers.hpp:252
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:719