All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
AutoDiffHelpers.hpp
1 /*
2  Copyright 2013 SINTEF ICT, Applied Mathematics.
3  Copyright 2015 IRIS
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_AUTODIFFHELPERS_HEADER_INCLUDED
22 #define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
23 
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>
30 
31 #include <opm/common/ErrorMacros.hpp>
32 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
33 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
34 #include <iostream>
35 #include <vector>
36 
37 namespace Opm
38 {
39 
40 // -------------------- class HelperOps --------------------
41 
44 struct HelperOps
45 {
46  typedef Eigen::SparseMatrix<double> M;
47  typedef AutoDiffBlock<double>::V V;
48 
50  typedef Eigen::Array<int, Eigen::Dynamic, 1> IFaces;
51  IFaces internal_faces;
52 
54  M ngrad;
56  M grad;
58  M caver;
60  M div;
66 
68  typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
69  TwoColInt nnc_cells;
70 
73 
76 
78  template<class Grid>
79  HelperOps(const Grid& grid, const NNC& nnc = NNC())
80  {
81  using namespace AutoDiffGrid;
82  const int nc = numCells(grid);
83  const int nf = numFaces(grid);
84  // Define some neighbourhood-derived helper arrays.
85 
86  TwoColInt nbi;
87  extractInternalFaces(grid, internal_faces, nbi);
88  const int num_internal = internal_faces.size();
89 
90  // handle non-neighboring connections
91  const bool has_nnc = nnc.hasNNC();
92  int numNNC = nnc.numNNC();
93 
94  // num_connections may also include non-neighboring connections
95  const int num_connections = num_internal + numNNC;
96  if (has_nnc) {
97  const int *cartDims = AutoDiffGrid::cartDims(grid);
98  const int numCartesianCells = cartDims[0] * cartDims[1] * cartDims[2];
99 
100  // the nnc's acts on global indicies and must be mapped to cell indicies
101  std::vector<int> global2localIdx(numCartesianCells,0);
102  for (int i = 0; i< nc; ++i) {
103  global2localIdx[ globalCell( grid )[i] ] = i;
104  }
105  nnc_cells.resize(numNNC,2);
106  nnc_trans.resize(numNNC);
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];
110  // store the nnc transmissibilities for later usage.
111  nnc_trans(i) = nnc.nncdata()[i].trans;
112  }
113  }
114 
115 
116  // std::cout << "nbi = \n" << nbi << std::endl;
117  // Create matrices.
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);
130  }
131  // add contribution from NNC
132  if (has_nnc) {
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);
138  }
139  }
140  ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end());
141  caver.setFromTriplets(caver_tri.begin(), caver_tri.end());
142  grad = -ngrad;
143  div = ngrad.transpose();
144 
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) {
149  if (nb(i,0) >= 0) {
150  fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
151  }
152  if (nb(i,1) >= 0) {
153  fullngrad_tri.emplace_back(i, nb(i,1), -1.0);
154  }
155  }
156  // add contribution from NNC
157  if (has_nnc) {
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);
161  }
162  }
163  fullngrad.resize(nf+numNNC, nc);
164  fullngrad.setFromTriplets(fullngrad_tri.begin(), fullngrad_tri.end());
165  fulldiv = fullngrad.transpose();
166 
167  if (has_nnc) {
168  connection_cells.resize(nbi.rows() + nnc_cells.rows(), 2);
169  connection_cells << nbi, nnc_cells;
170  } else {
171  connection_cells = nbi;
172  }
173  }
174 };
175 // -------------------- upwinding helper class --------------------
176 
177 
180  template <typename Scalar>
182  public:
183  typedef AutoDiffBlock<Scalar> ADB;
184 
185  template<class Grid>
186  UpwindSelector(const Grid& g,
187  const HelperOps& h,
188  const typename ADB::V& ifaceflux)
189  {
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);
195 
196  // num connections may possibly include NNCs
197  int num_nnc = h.nnc_trans.size();
198  int num_connections = nif + num_nnc;
199  assert(num_connections == ifaceflux.size());
200 
201  // Define selector structure.
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);
208 
209  assert ((c1 >= 0) && (c2 >= 0));
210 
211  // Select upwind cell.
212  const int c = (ifaceflux[iface] >= 0) ? c1 : c2;
213 
214  s.push_back(Triplet(iface, c, Scalar(1)));
215  }
216  if (num_nnc > 0) {
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)));
220  }
221  }
222 
223  // Assemble explicit selector operator.
224  select_.resize(num_connections, numCells(g));
225  select_.setFromTriplets(s.begin(), s.end());
226  }
227 
229  std::vector<ADB>
230  select(const std::vector<ADB>& xc) const
231  {
232  // Absence of counter-current flow means that the same
233  // selector applies to all quantities, 'x', defined per
234  // cell.
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)
238  {
239  xf.push_back(select_ * (*b));
240  }
241 
242  return xf;
243  }
244 
246  ADB select(const ADB& xc) const
247  {
248  return select_*xc;
249  }
250 
252  typename ADB::V select(const typename ADB::V& xc) const
253  {
254  return (select_*xc.matrix()).array();
255  }
256 
257  private:
258  Eigen::SparseMatrix<double> select_;
259  };
260 
261 
262 
263 namespace {
264 
265 
266  template <typename Scalar, class IntVec>
267  typename Eigen::SparseMatrix<Scalar>
268  constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
269  {
270  const int subset_size = indices.size();
271 
272  if (subset_size == 0) {
273  Eigen::SparseMatrix<Scalar> mat(full_size, 0);
274  return mat;
275  }
276 
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;
281  }
282  return mat;
283  }
284 
285 } // anon namespace
286 
287 
288 
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)
294 {
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 ] ];
300 
301  return std::move(ret);
302 }
303 
305 template <typename Scalar, class IntVec>
306 AutoDiffBlock<Scalar>
308  const IntVec& indices)
309 {
310  const Eigen::SparseMatrix<Scalar> sub
311  = constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose().eval();
312  return sub * x;
313 }
314 
315 
317 template <typename Scalar, class IntVec>
318 AutoDiffBlock<Scalar>
320  const IntVec& indices,
321  const int n)
322 {
323  return constructSupersetSparseMatrix<Scalar>(n, indices) * x;
324 }
325 
326 
327 
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,
333  const int n)
334 {
335  return constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
336 }
337 
338 
339 
343 inline
344 Eigen::SparseMatrix<double>
346 {
347  typedef Eigen::SparseMatrix<double> M;
348 
349  const int n = d.size();
350  M mat(n, n);
351  mat.reserve(Eigen::ArrayXi::Ones(n, 1));
352  for (M::Index i = 0; i < n; ++i) {
353  mat.insert(i, i) = d[i];
354  }
355 
356  return mat;
357 }
358 
359 
360 
361 
363  template <typename Scalar>
364  class Selector {
365  public:
366  typedef AutoDiffBlock<Scalar> ADB;
367 
368  enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN };
369 
370  Selector(const typename ADB::V& selection_basis,
371  CriterionForLeftElement crit = GreaterEqualZero)
372  {
373  using std::isnan;
374  // Define selector structure.
375  const int n = selection_basis.size();
376  // Over-reserving so we do not have to count.
377  left_elems_.reserve(n);
378  right_elems_.reserve(n);
379  for (int i = 0; i < n; ++i) {
380  bool chooseleft = false;
381  switch (crit) {
382  case GreaterEqualZero:
383  chooseleft = selection_basis[i] >= 0.0;
384  break;
385  case GreaterZero:
386  chooseleft = selection_basis[i] > 0.0;
387  break;
388  case Zero:
389  chooseleft = selection_basis[i] == 0.0;
390  break;
391  case NotEqualZero:
392  chooseleft = selection_basis[i] != 0.0;
393  break;
394  case LessZero:
395  chooseleft = selection_basis[i] < 0.0;
396  break;
397  case LessEqualZero:
398  chooseleft = selection_basis[i] <= 0.0;
399  break;
400  case NotNaN:
401  chooseleft = !isnan(selection_basis[i]);
402  break;
403  default:
404  OPM_THROW(std::logic_error, "No such criterion: " << crit);
405  }
406  if (chooseleft) {
407  left_elems_.push_back(i);
408  } else {
409  right_elems_.push_back(i);
410  }
411  }
412  }
413 
415  ADB select(const ADB& x1, const ADB& x2) const
416  {
417  if (right_elems_.empty()) {
418  return x1;
419  } else if (left_elems_.empty()) {
420  return x2;
421  } else {
422  return superset(subset(x1, left_elems_), left_elems_, x1.size())
423  + superset(subset(x2, right_elems_), right_elems_, x2.size());
424  }
425  }
426 
428  typename ADB::V select(const typename ADB::V& x1, const typename ADB::V& x2) const
429  {
430  if (right_elems_.empty()) {
431  return x1;
432  } else if (left_elems_.empty()) {
433  return x2;
434  } else {
435  return superset(subset(x1, left_elems_), left_elems_, x1.size())
436  + superset(subset(x2, right_elems_), right_elems_, x2.size());
437  }
438  }
439 
440  private:
441  std::vector<int> left_elems_;
442  std::vector<int> right_elems_;
443  };
444 
445 
446 
447 
449 template <class Matrix>
450 inline
451 void
452 collapseJacs(const AutoDiffBlock<double>& x, Matrix& jacobian)
453 {
454  typedef AutoDiffBlock<double> ADB;
455  const int nb = x.numBlocks();
456  typedef Eigen::Triplet<double> Tri;
457  int nnz = 0;
458  for (int block = 0; block < nb; ++block) {
459  nnz += x.derivative()[block].nonZeros();
460  }
461  std::vector<Tri> t;
462  t.reserve(nnz);
463  int block_col_start = 0;
464  for (int block = 0; block < nb; ++block) {
465  const ADB::M& jac1 = x.derivative()[block];
466  Eigen::SparseMatrix<double> jac;
467  jac1.toSparse(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,
473  i.value()));
474  }
475  }
476  }
477  block_col_start += jac.cols();
478  }
479  // Build final jacobian.
480  jacobian = Matrix(x.size(), block_col_start);
481  jacobian.setFromTriplets(t.begin(), t.end());
482 }
483 
484 
485 
487 inline
488 AutoDiffBlock<double>
490 {
491  Eigen::SparseMatrix<double> comb_jac;
492  collapseJacs( x, comb_jac );
493  // Build final jacobian.
494  typedef AutoDiffBlock<double> ADB;
495  std::vector<ADB::M> jacs(1);
496  jacs[0] = AutoDiffMatrix(std::move(comb_jac));
497  ADB::V val = x.value();
498  return ADB::function(std::move(val), std::move(jacs));
499 }
500 
501 
502 
504 inline
505 AutoDiffBlock<double>
507  const AutoDiffBlock<double>& y)
508 {
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) {
514  xind[i] = i;
515  }
516  std::vector<int> yind(ny);
517  for (int i = 0; i < ny; ++i) {
518  yind[i] = nx + i;
519  }
520  return superset(x, xind, n) + superset(y, yind, n);
521 }
522 
523 
524 
525 
528 inline
529 AutoDiffBlock<double>
531 {
532  typedef AutoDiffBlock<double> ADB;
533  if (x.empty()) {
534  return ADB::null();
535  }
536 
537  // Count sizes, nonzeros.
538  const int nx = x.size();
539  int size = 0;
540  int nnz = 0;
541  int elem_with_deriv = -1;
542  int num_blocks = 0;
543  for (int elem = 0; elem < nx; ++elem) {
544  size += x[elem].size();
545  if (x[elem].derivative().empty()) {
546  // No nnz contributions from this element.
547  continue;
548  } else {
549  if (elem_with_deriv == -1) {
550  elem_with_deriv = elem;
551  num_blocks = x[elem].numBlocks();
552  }
553  }
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");
556  }
557  for (int block = 0; block < num_blocks; ++block) {
558  nnz += x[elem].derivative()[block].nonZeros();
559  }
560  }
561  int num_cols = 0;
562  for (int block = 0; block < num_blocks; ++block) {
563  num_cols += x[elem_with_deriv].derivative()[block].cols();
564  }
565 
566  // Build value for result.
567  ADB::V val(size);
568  int pos = 0;
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();
572  pos += loc_size;
573  }
574  assert(pos == size);
575 
576  // Return a constant if we have no derivatives at all.
577  if (num_blocks == 0) {
578  return ADB::constant(std::move(val));
579  }
580 
581  // Set up for batch insertion of all Jacobian elements.
582  typedef Eigen::SparseMatrix<double> M;
583  typedef Eigen::Triplet<double> Tri;
584  std::vector<Tri> t;
585  t.reserve(nnz);
586  {
587  int block_row_start = 0;
588  M jac;
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,
598  i.value()));
599  }
600  }
601  block_col_start += jac.cols();
602  }
603  }
604  block_row_start += x[elem].size();
605  }
606  }
607 
608  // Build final jacobian.
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));
614 
615  // Use move semantics to return result efficiently.
616  return ADB::function(std::move(val), std::move(jac));
617 }
618 
619 
620 
621 
622 
623 class Span
624 {
625 public:
626  explicit Span(const int num)
627  : num_(num),
628  stride_(1),
629  start_(0)
630  {
631  }
632  Span(const int num, const int stride, const int start)
633  : num_(num),
634  stride_(stride),
635  start_(start)
636  {
637  }
638  int operator[](const int i) const
639  {
640  assert(i >= 0 && i < num_);
641  return start_ + i*stride_;
642  }
643  int size() const
644  {
645  return num_;
646  }
647 
648 
650  {
651  public:
652  SpanIterator(const Span* span, const int index)
653  : span_(span),
654  index_(index)
655  {
656  }
657  SpanIterator operator++()
658  {
659  ++index_;
660  return *this;
661  }
662  SpanIterator operator++(int)
663  {
664  SpanIterator before_increment(*this);
665  ++index_;
666  return before_increment;
667  }
668  bool operator<(const SpanIterator& rhs) const
669  {
670  assert(span_ == rhs.span_);
671  return index_ < rhs.index_;
672  }
673  bool operator==(const SpanIterator& rhs) const
674  {
675  assert(span_ == rhs.span_);
676  return index_ == rhs.index_;
677  }
678  bool operator!=(const SpanIterator& rhs) const
679  {
680  assert(span_ == rhs.span_);
681  return index_ != rhs.index_;
682  }
683  int operator*()
684  {
685  return (*span_)[index_];
686  }
687  private:
688  const Span* span_;
689  int index_;
690  };
691 
692  typedef SpanIterator iterator;
694 
695  SpanIterator begin() const
696  {
697  return SpanIterator(this, 0);
698  }
699 
700  SpanIterator end() const
701  {
702  return SpanIterator(this, num_);
703  }
704 
705  bool operator==(const Span& rhs)
706  {
707  return num_ == rhs.num_ && start_ == rhs.start_ && stride_ == rhs.stride_;
708  }
709 
710 private:
711  const int num_;
712  const int stride_;
713  const int start_;
714 };
715 
716 
717 
719 inline Eigen::ArrayXd sign (const Eigen::ArrayXd& x)
720 {
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);
725  }
726  return retval;
727 }
728 
729 } // namespace Opm
730 
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&#39; (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&#39; cells (face or nnc).
Definition: AutoDiffHelpers.hpp:75
M ngrad
Extract for each internal face the difference of its adjacent cells&#39; values (first - second)...
Definition: AutoDiffHelpers.hpp:54
M caver
Extract for each face the average of its adjacent cells&#39; 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&#39; values (second - first).
Definition: AutoDiffHelpers.hpp:56
M fulldiv
Extract for each cell the sum of all its adjacent faces&#39; (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&#39; 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