00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_AUTODIFFHELPERS_HEADER_INCLUDED
00022 #define OPM_AUTODIFFHELPERS_HEADER_INCLUDED
00023
00024 #include <opm/autodiff/AutoDiffBlock.hpp>
00025 #include <opm/autodiff/GridHelpers.hpp>
00026 #include <opm/autodiff/GeoProps.hpp>
00027 #include <opm/core/grid.h>
00028 #include <opm/core/grid/PinchProcessor.hpp>
00029 #include <opm/core/props/rock/RockFromDeck.hpp>
00030
00031 #include <opm/common/ErrorMacros.hpp>
00032 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
00033 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00034 #include <iostream>
00035 #include <vector>
00036
00037 namespace Opm
00038 {
00039
00040
00041
00044 struct HelperOps
00045 {
00046 typedef Eigen::SparseMatrix<double> M;
00047 typedef AutoDiffBlock<double>::V V;
00048
00050 typedef Eigen::Array<int, Eigen::Dynamic, 1> IFaces;
00051 IFaces internal_faces;
00052
00054 M ngrad;
00056 M grad;
00058 M caver;
00060 M div;
00063 M fullngrad;
00065 M fulldiv;
00066
00068 typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
00069 TwoColInt nnc_cells;
00070
00072 V nnc_trans;
00073
00075 TwoColInt connection_cells;
00076
00078 template<class Grid>
00079 HelperOps(const Grid& grid, const NNC& nnc = NNC())
00080 {
00081 using namespace AutoDiffGrid;
00082 const int nc = numCells(grid);
00083 const int nf = numFaces(grid);
00084
00085
00086 TwoColInt nbi;
00087 extractInternalFaces(grid, internal_faces, nbi);
00088 const int num_internal = internal_faces.size();
00089
00090
00091 const bool has_nnc = nnc.hasNNC();
00092 int numNNC = nnc.numNNC();
00093
00094
00095 const int num_connections = num_internal + numNNC;
00096 if (has_nnc) {
00097 const int *cartDims = AutoDiffGrid::cartDims(grid);
00098 const int numCartesianCells = cartDims[0] * cartDims[1] * cartDims[2];
00099
00100
00101 std::vector<int> global2localIdx(numCartesianCells,0);
00102 for (int i = 0; i< nc; ++i) {
00103 global2localIdx[ globalCell( grid )[i] ] = i;
00104 }
00105 nnc_cells.resize(numNNC,2);
00106 nnc_trans.resize(numNNC);
00107 for (int i = 0; i < numNNC; ++i) {
00108 nnc_cells(i,0) = global2localIdx[nnc.nncdata()[i].cell1];
00109 nnc_cells(i,1) = global2localIdx[nnc.nncdata()[i].cell2];
00110
00111 nnc_trans(i) = nnc.nncdata()[i].trans;
00112 }
00113 }
00114
00115
00116
00117
00118 ngrad.resize(num_connections, nc);
00119 caver.resize(num_connections, nc);
00120 typedef Eigen::Triplet<double> Tri;
00121 std::vector<Tri> ngrad_tri;
00122 std::vector<Tri> caver_tri;
00123 ngrad_tri.reserve(2*num_connections);
00124 caver_tri.reserve(2*num_connections);
00125 for (int i = 0; i < num_internal; ++i) {
00126 ngrad_tri.emplace_back(i, nbi(i,0), 1.0);
00127 ngrad_tri.emplace_back(i, nbi(i,1), -1.0);
00128 caver_tri.emplace_back(i, nbi(i,0), 0.5);
00129 caver_tri.emplace_back(i, nbi(i,1), 0.5);
00130 }
00131
00132 if (has_nnc) {
00133 for (int i = 0; i < numNNC; ++i) {
00134 ngrad_tri.emplace_back(i+num_internal, nnc_cells(i,0), 1.0);
00135 ngrad_tri.emplace_back(i+num_internal, nnc_cells(i,1), -1.0);
00136 caver_tri.emplace_back(i+num_internal, nnc_cells(i,0), 0.5);
00137 caver_tri.emplace_back(i+num_internal, nnc_cells(i,1), 0.5);
00138 }
00139 }
00140 ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end());
00141 caver.setFromTriplets(caver_tri.begin(), caver_tri.end());
00142 grad = -ngrad;
00143 div = ngrad.transpose();
00144
00145 std::vector<Tri> fullngrad_tri;
00146 fullngrad_tri.reserve(2*(nf+numNNC));
00147 typename ADFaceCellTraits<Grid>::Type nb = faceCellsToEigen(grid);
00148 for (int i = 0; i < nf; ++i) {
00149 if (nb(i,0) >= 0) {
00150 fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
00151 }
00152 if (nb(i,1) >= 0) {
00153 fullngrad_tri.emplace_back(i, nb(i,1), -1.0);
00154 }
00155 }
00156
00157 if (has_nnc) {
00158 for (int i = 0; i < numNNC; ++i) {
00159 fullngrad_tri.emplace_back(i+nf, nnc_cells(i,0), 1.0);
00160 fullngrad_tri.emplace_back(i+nf, nnc_cells(i,1), -1.0);
00161 }
00162 }
00163 fullngrad.resize(nf+numNNC, nc);
00164 fullngrad.setFromTriplets(fullngrad_tri.begin(), fullngrad_tri.end());
00165 fulldiv = fullngrad.transpose();
00166
00167 if (has_nnc) {
00168 connection_cells.resize(nbi.rows() + nnc_cells.rows(), 2);
00169 connection_cells << nbi, nnc_cells;
00170 } else {
00171 connection_cells = nbi;
00172 }
00173 }
00174 };
00175
00176
00177
00180 template <typename Scalar>
00181 class UpwindSelector {
00182 public:
00183 typedef AutoDiffBlock<Scalar> ADB;
00184
00185 template<class Grid>
00186 UpwindSelector(const Grid& g,
00187 const HelperOps& h,
00188 const typename ADB::V& ifaceflux)
00189 {
00190 using namespace AutoDiffGrid;
00191 typedef HelperOps::IFaces::Index IFIndex;
00192 const IFIndex nif = h.internal_faces.size();
00193 typename ADFaceCellTraits<Grid>::Type
00194 face_cells = faceCellsToEigen(g);
00195
00196
00197 int num_nnc = h.nnc_trans.size();
00198 int num_connections = nif + num_nnc;
00199 assert(num_connections == ifaceflux.size());
00200
00201
00202 typedef typename Eigen::Triplet<Scalar> Triplet;
00203 std::vector<Triplet> s; s.reserve(num_connections);
00204 for (IFIndex iface = 0; iface < nif; ++iface) {
00205 const int f = h.internal_faces[iface];
00206 const int c1 = face_cells(f,0);
00207 const int c2 = face_cells(f,1);
00208
00209 assert ((c1 >= 0) && (c2 >= 0));
00210
00211
00212 const int c = (ifaceflux[iface] >= 0) ? c1 : c2;
00213
00214 s.push_back(Triplet(iface, c, Scalar(1)));
00215 }
00216 if (num_nnc > 0) {
00217 for (int i = 0; i < num_nnc; ++i) {
00218 const int c = (ifaceflux[i+nif] >= 0) ? h.nnc_cells(i,0) : h.nnc_cells(i,1);
00219 s.push_back(Triplet(i+nif,c,Scalar(1)));
00220 }
00221 }
00222
00223
00224 select_.resize(num_connections, numCells(g));
00225 select_.setFromTriplets(s.begin(), s.end());
00226 }
00227
00229 std::vector<ADB>
00230 select(const std::vector<ADB>& xc) const
00231 {
00232
00233
00234
00235 std::vector<ADB> xf; xf.reserve(xc.size());
00236 for (typename std::vector<ADB>::const_iterator
00237 b = xc.begin(), e = xc.end(); b != e; ++b)
00238 {
00239 xf.push_back(select_ * (*b));
00240 }
00241
00242 return xf;
00243 }
00244
00246 ADB select(const ADB& xc) const
00247 {
00248 return select_*xc;
00249 }
00250
00252 typename ADB::V select(const typename ADB::V& xc) const
00253 {
00254 return (select_*xc.matrix()).array();
00255 }
00256
00257 private:
00258 Eigen::SparseMatrix<double> select_;
00259 };
00260
00261
00262
00263 namespace {
00264
00265
00266 template <typename Scalar, class IntVec>
00267 typename Eigen::SparseMatrix<Scalar>
00268 constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
00269 {
00270 const int subset_size = indices.size();
00271
00272 if (subset_size == 0) {
00273 Eigen::SparseMatrix<Scalar> mat(full_size, 0);
00274 return mat;
00275 }
00276
00277 Eigen::SparseMatrix<Scalar> mat(full_size, subset_size);
00278 mat.reserve(Eigen::VectorXi::Constant(subset_size, 1));
00279 for (int i = 0; i < subset_size; ++i) {
00280 mat.insert(indices[i], i) = 1;
00281 }
00282 return mat;
00283 }
00284
00285 }
00286
00287
00288
00290 template <typename Scalar, class IntVec>
00291 Eigen::Array<Scalar, Eigen::Dynamic, 1>
00292 subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
00293 const IntVec& indices)
00294 {
00295 typedef typename Eigen::Array<Scalar, Eigen::Dynamic, 1>::Index Index;
00296 const Index size = indices.size();
00297 Eigen::Array<Scalar, Eigen::Dynamic, 1> ret( size );
00298 for( Index i=0; i<size; ++i )
00299 ret[ i ] = x[ indices[ i ] ];
00300
00301 return std::move(ret);
00302 }
00303
00305 template <typename Scalar, class IntVec>
00306 AutoDiffBlock<Scalar>
00307 subset(const AutoDiffBlock<Scalar>& x,
00308 const IntVec& indices)
00309 {
00310 const Eigen::SparseMatrix<Scalar> sub
00311 = constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose().eval();
00312 return sub * x;
00313 }
00314
00315
00317 template <typename Scalar, class IntVec>
00318 AutoDiffBlock<Scalar>
00319 superset(const AutoDiffBlock<Scalar>& x,
00320 const IntVec& indices,
00321 const int n)
00322 {
00323 return constructSupersetSparseMatrix<Scalar>(n, indices) * x;
00324 }
00325
00326
00327
00329 template <typename Scalar, class IntVec>
00330 Eigen::Array<Scalar, Eigen::Dynamic, 1>
00331 superset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
00332 const IntVec& indices,
00333 const int n)
00334 {
00335 return constructSupersetSparseMatrix<Scalar>(n, indices) * x.matrix();
00336 }
00337
00338
00339
00343 inline
00344 Eigen::SparseMatrix<double>
00345 spdiag(const AutoDiffBlock<double>::V& d)
00346 {
00347 typedef Eigen::SparseMatrix<double> M;
00348
00349 const int n = d.size();
00350 M mat(n, n);
00351 mat.reserve(Eigen::ArrayXi::Ones(n, 1));
00352 for (M::Index i = 0; i < n; ++i) {
00353 mat.insert(i, i) = d[i];
00354 }
00355
00356 return mat;
00357 }
00358
00359
00360
00361
00363 template <typename Scalar>
00364 class Selector {
00365 public:
00366 typedef AutoDiffBlock<Scalar> ADB;
00367
00368 enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN };
00369
00370 Selector(const typename ADB::V& selection_basis,
00371 CriterionForLeftElement crit = GreaterEqualZero)
00372 {
00373 using std::isnan;
00374
00375 const int n = selection_basis.size();
00376
00377 left_elems_.reserve(n);
00378 right_elems_.reserve(n);
00379 for (int i = 0; i < n; ++i) {
00380 bool chooseleft = false;
00381 switch (crit) {
00382 case GreaterEqualZero:
00383 chooseleft = selection_basis[i] >= 0.0;
00384 break;
00385 case GreaterZero:
00386 chooseleft = selection_basis[i] > 0.0;
00387 break;
00388 case Zero:
00389 chooseleft = selection_basis[i] == 0.0;
00390 break;
00391 case NotEqualZero:
00392 chooseleft = selection_basis[i] != 0.0;
00393 break;
00394 case LessZero:
00395 chooseleft = selection_basis[i] < 0.0;
00396 break;
00397 case LessEqualZero:
00398 chooseleft = selection_basis[i] <= 0.0;
00399 break;
00400 case NotNaN:
00401 chooseleft = !isnan(selection_basis[i]);
00402 break;
00403 default:
00404 OPM_THROW(std::logic_error, "No such criterion: " << crit);
00405 }
00406 if (chooseleft) {
00407 left_elems_.push_back(i);
00408 } else {
00409 right_elems_.push_back(i);
00410 }
00411 }
00412 }
00413
00415 ADB select(const ADB& x1, const ADB& x2) const
00416 {
00417 if (right_elems_.empty()) {
00418 return x1;
00419 } else if (left_elems_.empty()) {
00420 return x2;
00421 } else {
00422 return superset(subset(x1, left_elems_), left_elems_, x1.size())
00423 + superset(subset(x2, right_elems_), right_elems_, x2.size());
00424 }
00425 }
00426
00428 typename ADB::V select(const typename ADB::V& x1, const typename ADB::V& x2) const
00429 {
00430 if (right_elems_.empty()) {
00431 return x1;
00432 } else if (left_elems_.empty()) {
00433 return x2;
00434 } else {
00435 return superset(subset(x1, left_elems_), left_elems_, x1.size())
00436 + superset(subset(x2, right_elems_), right_elems_, x2.size());
00437 }
00438 }
00439
00440 private:
00441 std::vector<int> left_elems_;
00442 std::vector<int> right_elems_;
00443 };
00444
00445
00446
00447
00449 template <class Matrix>
00450 inline
00451 void
00452 collapseJacs(const AutoDiffBlock<double>& x, Matrix& jacobian)
00453 {
00454 typedef AutoDiffBlock<double> ADB;
00455 const int nb = x.numBlocks();
00456 typedef Eigen::Triplet<double> Tri;
00457 int nnz = 0;
00458 for (int block = 0; block < nb; ++block) {
00459 nnz += x.derivative()[block].nonZeros();
00460 }
00461 std::vector<Tri> t;
00462 t.reserve(nnz);
00463 int block_col_start = 0;
00464 for (int block = 0; block < nb; ++block) {
00465 const ADB::M& jac1 = x.derivative()[block];
00466 Eigen::SparseMatrix<double> jac;
00467 jac1.toSparse(jac);
00468 for (Eigen::SparseMatrix<double>::Index k = 0; k < jac.outerSize(); ++k) {
00469 for (Eigen::SparseMatrix<double>::InnerIterator i(jac, k); i ; ++i) {
00470 if (i.value() != 0.0) {
00471 t.push_back(Tri(i.row(),
00472 i.col() + block_col_start,
00473 i.value()));
00474 }
00475 }
00476 }
00477 block_col_start += jac.cols();
00478 }
00479
00480 jacobian = Matrix(x.size(), block_col_start);
00481 jacobian.setFromTriplets(t.begin(), t.end());
00482 }
00483
00484
00485
00487 inline
00488 AutoDiffBlock<double>
00489 collapseJacs(const AutoDiffBlock<double>& x)
00490 {
00491 Eigen::SparseMatrix<double> comb_jac;
00492 collapseJacs( x, comb_jac );
00493
00494 typedef AutoDiffBlock<double> ADB;
00495 std::vector<ADB::M> jacs(1);
00496 jacs[0] = AutoDiffMatrix(std::move(comb_jac));
00497 ADB::V val = x.value();
00498 return ADB::function(std::move(val), std::move(jacs));
00499 }
00500
00501
00502
00504 inline
00505 AutoDiffBlock<double>
00506 vertcat(const AutoDiffBlock<double>& x,
00507 const AutoDiffBlock<double>& y)
00508 {
00509 const int nx = x.size();
00510 const int ny = y.size();
00511 const int n = nx + ny;
00512 std::vector<int> xind(nx);
00513 for (int i = 0; i < nx; ++i) {
00514 xind[i] = i;
00515 }
00516 std::vector<int> yind(ny);
00517 for (int i = 0; i < ny; ++i) {
00518 yind[i] = nx + i;
00519 }
00520 return superset(x, xind, n) + superset(y, yind, n);
00521 }
00522
00523
00524
00525
00528 inline
00529 AutoDiffBlock<double>
00530 vertcatCollapseJacs(const std::vector<AutoDiffBlock<double> >& x)
00531 {
00532 typedef AutoDiffBlock<double> ADB;
00533 if (x.empty()) {
00534 return ADB::null();
00535 }
00536
00537
00538 const int nx = x.size();
00539 int size = 0;
00540 int nnz = 0;
00541 int elem_with_deriv = -1;
00542 int num_blocks = 0;
00543 for (int elem = 0; elem < nx; ++elem) {
00544 size += x[elem].size();
00545 if (x[elem].derivative().empty()) {
00546
00547 continue;
00548 } else {
00549 if (elem_with_deriv == -1) {
00550 elem_with_deriv = elem;
00551 num_blocks = x[elem].numBlocks();
00552 }
00553 }
00554 if (x[elem].blockPattern() != x[elem_with_deriv].blockPattern()) {
00555 OPM_THROW(std::runtime_error, "vertcatCollapseJacs(): all arguments must have the same block pattern");
00556 }
00557 for (int block = 0; block < num_blocks; ++block) {
00558 nnz += x[elem].derivative()[block].nonZeros();
00559 }
00560 }
00561 int num_cols = 0;
00562 for (int block = 0; block < num_blocks; ++block) {
00563 num_cols += x[elem_with_deriv].derivative()[block].cols();
00564 }
00565
00566
00567 ADB::V val(size);
00568 int pos = 0;
00569 for (int elem = 0; elem < nx; ++elem) {
00570 const int loc_size = x[elem].size();
00571 val.segment(pos, loc_size) = x[elem].value();
00572 pos += loc_size;
00573 }
00574 assert(pos == size);
00575
00576
00577 if (num_blocks == 0) {
00578 return ADB::constant(std::move(val));
00579 }
00580
00581
00582 typedef Eigen::SparseMatrix<double> M;
00583 typedef Eigen::Triplet<double> Tri;
00584 std::vector<Tri> t;
00585 t.reserve(nnz);
00586 {
00587 int block_row_start = 0;
00588 M jac;
00589 for (int elem = 0; elem < nx; ++elem) {
00590 int block_col_start = 0;
00591 if (!x[elem].derivative().empty()) {
00592 for (int block = 0; block < num_blocks; ++block) {
00593 x[elem].derivative()[block].toSparse(jac);
00594 for (M::Index k = 0; k < jac.outerSize(); ++k) {
00595 for (M::InnerIterator i(jac, k); i ; ++i) {
00596 t.push_back(Tri(i.row() + block_row_start,
00597 i.col() + block_col_start,
00598 i.value()));
00599 }
00600 }
00601 block_col_start += jac.cols();
00602 }
00603 }
00604 block_row_start += x[elem].size();
00605 }
00606 }
00607
00608
00609 M comb_jac = M(size, num_cols);
00610 comb_jac.reserve(nnz);
00611 comb_jac.setFromTriplets(t.begin(), t.end());
00612 std::vector<ADB::M> jac(1);
00613 jac[0] = ADB::M(std::move(comb_jac));
00614
00615
00616 return ADB::function(std::move(val), std::move(jac));
00617 }
00618
00619
00620
00621
00622
00623 class Span
00624 {
00625 public:
00626 explicit Span(const int num)
00627 : num_(num),
00628 stride_(1),
00629 start_(0)
00630 {
00631 }
00632 Span(const int num, const int stride, const int start)
00633 : num_(num),
00634 stride_(stride),
00635 start_(start)
00636 {
00637 }
00638 int operator[](const int i) const
00639 {
00640 assert(i >= 0 && i < num_);
00641 return start_ + i*stride_;
00642 }
00643 int size() const
00644 {
00645 return num_;
00646 }
00647
00648
00649 class SpanIterator
00650 {
00651 public:
00652 SpanIterator(const Span* span, const int index)
00653 : span_(span),
00654 index_(index)
00655 {
00656 }
00657 SpanIterator operator++()
00658 {
00659 ++index_;
00660 return *this;
00661 }
00662 SpanIterator operator++(int)
00663 {
00664 SpanIterator before_increment(*this);
00665 ++index_;
00666 return before_increment;
00667 }
00668 bool operator<(const SpanIterator& rhs) const
00669 {
00670 assert(span_ == rhs.span_);
00671 return index_ < rhs.index_;
00672 }
00673 bool operator==(const SpanIterator& rhs) const
00674 {
00675 assert(span_ == rhs.span_);
00676 return index_ == rhs.index_;
00677 }
00678 bool operator!=(const SpanIterator& rhs) const
00679 {
00680 assert(span_ == rhs.span_);
00681 return index_ != rhs.index_;
00682 }
00683 int operator*()
00684 {
00685 return (*span_)[index_];
00686 }
00687 private:
00688 const Span* span_;
00689 int index_;
00690 };
00691
00692 typedef SpanIterator iterator;
00693 typedef SpanIterator const_iterator;
00694
00695 SpanIterator begin() const
00696 {
00697 return SpanIterator(this, 0);
00698 }
00699
00700 SpanIterator end() const
00701 {
00702 return SpanIterator(this, num_);
00703 }
00704
00705 bool operator==(const Span& rhs)
00706 {
00707 return num_ == rhs.num_ && start_ == rhs.start_ && stride_ == rhs.stride_;
00708 }
00709
00710 private:
00711 const int num_;
00712 const int stride_;
00713 const int start_;
00714 };
00715
00716
00717
00719 inline Eigen::ArrayXd sign (const Eigen::ArrayXd& x)
00720 {
00721 const int n = x.size();
00722 Eigen::ArrayXd retval(n);
00723 for (int i = 0; i < n; ++i) {
00724 retval[i] = x[i] < 0.0 ? -1.0 : (x[i] > 0.0 ? 1.0 : 0.0);
00725 }
00726 return retval;
00727 }
00728
00729 }
00730
00731 #endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED