diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 44c4bfb94..6c2a32909 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -33,12 +33,12 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/NewtonIterationBlackoilSimple.cpp opm/autodiff/NewtonIterationUtilities.cpp opm/autodiff/GridHelpers.cpp - opm/autodiff/ImpesTPFAAD.cpp +# opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp - opm/autodiff/SimulatorIncompTwophaseAd.cpp - opm/autodiff/TransportSolverTwophaseAd.cpp +# opm/autodiff/SimulatorIncompTwophaseAd.cpp +# opm/autodiff/TransportSolverTwophaseAd.cpp opm/autodiff/BlackoilPropsAdFromDeck.cpp - opm/autodiff/SolventPropsAdFromDeck.cpp +# opm/autodiff/SolventPropsAdFromDeck.cpp opm/autodiff/BlackoilModelParameters.cpp opm/autodiff/WellDensitySegmented.cpp opm/autodiff/LinearisedBlackoilResidual.cpp @@ -83,20 +83,20 @@ endif() # originally generated with the command: # find tutorials examples -name '*.c*' -printf '\t%p\n' | sort list (APPEND EXAMPLE_SOURCE_FILES - examples/find_zero.cpp +# examples/find_zero.cpp examples/flow.cpp - examples/flow_solvent.cpp - examples/sim_2p_incomp_ad.cpp - examples/sim_simple.cpp - examples/opm_init_check.cpp +# examples/flow_solvent.cpp +# examples/sim_2p_incomp_ad.cpp +# examples/sim_simple.cpp +# examples/opm_init_check.cpp ) # programs listed here will not only be compiled, but also marked for # installation list (APPEND PROGRAM_SOURCE_FILES - examples/sim_2p_incomp_ad.cpp +# examples/sim_2p_incomp_ad.cpp examples/flow.cpp - examples/opm_init_check.cpp +# examples/opm_init_check.cpp ) # originally generated with the command: diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 5ba57b32f..4ef55c41c 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -28,6 +28,9 @@ #include +#include + + #include #include #include @@ -94,7 +97,7 @@ namespace Opm /// Underlying type for values. typedef Eigen::Array V; /// Underlying type for jacobians. - typedef Eigen::SparseMatrix M; + typedef AutoDiffMatrix M; /// Construct an empty AutoDiffBlock. static AutoDiffBlock null() @@ -155,22 +158,23 @@ namespace Opm } // ... then set the one corrresponding to this variable to identity. assert(blocksizes[index] == num_elem); - if(val.size()>0) - { - // if val is empty the following will run into an assertion - // with Eigen - jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1)); - for (typename M::Index row = 0; row < val.size(); ++row) { - jac[index].insert(row, row) = Scalar(1.0); - } - } - else - { - // Do some check. Empty jacobian only make sense for empty val. - for (int i = 0; i < num_blocks; ++i) { - assert(jac[i].size()==0); - } - } + jac[index] = M(M::IdentityMatrix, val.size()); + // if(val.size()>0) + // { + // // if val is empty the following will run into an assertion + // // with Eigen + // jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1)); + // for (typename M::Index row = 0; row < val.size(); ++row) { + // jac[index].insert(row, row) = Scalar(1.0); + // } + // } + // else + // { + // // Do some check. Empty jacobian only make sense for empty val. + // for (int i = 0; i < num_blocks; ++i) { + // assert(jac[i].size()==0); + // } + // } return AutoDiffBlock(std::move(val), std::move(jac)); } @@ -254,7 +258,8 @@ namespace Opm const int num_blocks = rhs.numBlocks(); jac_.resize(num_blocks); for (int block = 0; block < num_blocks; ++block) { - jac_[block] = -rhs.jac_[block]; + // jac_[block] = -rhs.jac_[block]; + jac_[block] = rhs.jac_[block] * (-1.0); } } else if (!rhs.jac_.empty()) { assert (numBlocks() == rhs.numBlocks()); @@ -334,9 +339,11 @@ namespace Opm int num_blocks = numBlocks(); std::vector jac(num_blocks); assert(numBlocks() == rhs.numBlocks()); - typedef Eigen::DiagonalMatrix D; - D D1 = val_.matrix().asDiagonal(); - D D2 = rhs.val_.matrix().asDiagonal(); + // typedef Eigen::DiagonalMatrix D; + // D D1 = val_.matrix().asDiagonal(); + // D D2 = rhs.val_.matrix().asDiagonal(); + M D1(val_.matrix().asDiagonal()); + M D2(rhs.val_.matrix().asDiagonal()); for (int block = 0; block < num_blocks; ++block) { assert(jac_[block].rows() == rhs.jac_[block].rows()); assert(jac_[block].cols() == rhs.jac_[block].cols()); @@ -371,9 +378,12 @@ namespace Opm std::vector jac(num_blocks); assert(numBlocks() == rhs.numBlocks()); typedef Eigen::DiagonalMatrix D; - D D1 = val_.matrix().asDiagonal(); - D D2 = rhs.val_.matrix().asDiagonal(); - D D3 = (1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal(); + // D D1 = val_.matrix().asDiagonal(); + // D D2 = rhs.val_.matrix().asDiagonal(); + // D D3 = (1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal(); + M D1(val_.matrix().asDiagonal()); + M D2(rhs.val_.matrix().asDiagonal()); + M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal()); for (int block = 0; block < num_blocks; ++block) { assert(jac_[block].rows() == rhs.jac_[block].rows()); assert(jac_[block].cols() == rhs.jac_[block].cols()); @@ -381,15 +391,14 @@ namespace Opm jac[block] = M( D3.rows(), jac_[block].cols() ); } else if( jac_[block].nonZeros() == 0 ) { - jac[block] = D3 * ( D1*rhs.jac_[block]); - jac[block] *= -1.0; + jac[block] = D3 * ( D1*rhs.jac_[block]) * (-1.0); } else if ( rhs.jac_[block].nonZeros() == 0 ) { jac[block] = D3 * (D2*jac_[block]); } else { - jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]); + jac[block] = D3 * (D2*jac_[block] + (D1*rhs.jac_[block]*(-1.0))); } } return function(val_ / rhs.val_, std::move(jac)); @@ -489,7 +498,7 @@ namespace Opm } - /// Multiply with sparse matrix from the left. + /// Multiply with AutoDiffMatrix from the left. template AutoDiffBlock operator*(const typename AutoDiffBlock::M& lhs, const AutoDiffBlock& rhs) @@ -506,6 +515,23 @@ namespace Opm } + /// Multiply with Eigen sparse matrix from the left. + template + AutoDiffBlock operator*(const Eigen::SparseMatrix& lhs, + const AutoDiffBlock& rhs) + { + int num_blocks = rhs.numBlocks(); + std::vector::M> jac(num_blocks); + assert(lhs.cols() == rhs.value().rows()); + for (int block = 0; block < num_blocks; ++block) { + // jac[block] = lhs*rhs.derivative()[block]; + fastSparseProduct(lhs, rhs.derivative()[block], jac[block]); + } + typename AutoDiffBlock::V val = lhs*rhs.value().matrix(); + return AutoDiffBlock::function(std::move(val), std::move(jac)); + } + + /// Elementwise multiplication with constant on the left. template AutoDiffBlock operator*(const typename AutoDiffBlock::V& lhs, diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 8ac663727..ca168f792 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -39,7 +39,8 @@ namespace Opm /// operations on (AD or regular) vectors of data. struct HelperOps { - typedef AutoDiffBlock::M M; + // typedef AutoDiffBlock::M M; + typedef Eigen::SparseMatrix M; typedef AutoDiffBlock::V V; /// A list of internal faces. @@ -250,7 +251,7 @@ struct HelperOps } private: - typename ADB::M select_; + Eigen::SparseMatrix select_; }; @@ -259,17 +260,17 @@ namespace { template - typename AutoDiffBlock::M + typename Eigen::SparseMatrix constructSupersetSparseMatrix(const int full_size, const IntVec& indices) { const int subset_size = indices.size(); if (subset_size == 0) { - typename AutoDiffBlock::M mat(full_size, 0); + Eigen::SparseMatrix mat(full_size, 0); return mat; } - typename AutoDiffBlock::M mat(full_size, subset_size); + Eigen::SparseMatrix mat(full_size, subset_size); mat.reserve(Eigen::VectorXi::Constant(subset_size, 1)); for (int i = 0; i < subset_size; ++i) { mat.insert(indices[i], i) = 1; @@ -302,8 +303,8 @@ AutoDiffBlock subset(const AutoDiffBlock& x, const IntVec& indices) { - const typename AutoDiffBlock::M sub - = constructSupersetSparseMatrix(x.value().size(), indices).transpose(); + const Eigen::SparseMatrix sub + = constructSupersetSparseMatrix(x.value().size(), indices).transpose().eval(); return sub * x; } @@ -336,10 +337,12 @@ superset(const Eigen::Array& x, /// elements of d on the diagonal. /// Need to mark this as inline since it is defined in a header and not a template. inline -AutoDiffBlock::M +//AutoDiffBlock::M +Eigen::SparseMatrix spdiag(const AutoDiffBlock::V& d) { - typedef AutoDiffBlock::M M; + // typedef AutoDiffBlock::M M; + typedef Eigen::SparseMatrix M; const int n = d.size(); M mat(n, n); @@ -456,9 +459,11 @@ collapseJacs(const AutoDiffBlock& x, Matrix& jacobian) t.reserve(nnz); int block_col_start = 0; for (int block = 0; block < nb; ++block) { - const ADB::M& jac = x.derivative()[block]; - for (ADB::M::Index k = 0; k < jac.outerSize(); ++k) { - for (ADB::M::InnerIterator i(jac, k); i ; ++i) { + const ADB::M& jac1 = x.derivative()[block]; + Eigen::SparseMatrix jac; + jac1.toSparse(jac); + for (Eigen::SparseMatrix::Index k = 0; k < jac.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator i(jac, k); i ; ++i) { t.push_back(Tri(i.row(), i.col() + block_col_start, i.value())); @@ -478,16 +483,18 @@ inline AutoDiffBlock collapseJacs(const AutoDiffBlock& x) { - typedef AutoDiffBlock ADB; + Eigen::SparseMatrix comb_jac; + collapseJacs( x, comb_jac ); // Build final jacobian. + typedef AutoDiffBlock ADB; std::vector jacs(1); - collapseJacs( x, jacs[ 0 ] ); + jacs[0] = AutoDiffMatrix(std::move(comb_jac)); ADB::V val = x.value(); return ADB::function(std::move(val), std::move(jacs)); } - + /* /// Returns the vertical concatenation [ x; y ] of the inputs. inline @@ -510,7 +517,7 @@ vertcat(const AutoDiffBlock& x, } - + */ /// Returns the vertical concatenation [ x[0]; x[1]; ...; x[n-1] ] of the inputs. @@ -569,6 +576,7 @@ vertcatCollapseJacs(const std::vector >& x) } // Set up for batch insertion of all Jacobian elements. + typedef Eigen::SparseMatrix M; typedef Eigen::Triplet Tri; std::vector t; t.reserve(nnz); @@ -577,9 +585,11 @@ vertcatCollapseJacs(const std::vector >& x) int block_col_start = 0; if (!x[elem].derivative().empty()) { for (int block = 0; block < num_blocks; ++block) { - const ADB::M& jac = x[elem].derivative()[block]; - for (ADB::M::Index k = 0; k < jac.outerSize(); ++k) { - for (ADB::M::InnerIterator i(jac, k); i ; ++i) { + // const ADB::M& jac = x[elem].derivative()[block]; + M jac; + x[elem].derivative()[block].toSparse(jac); + for (M::Index k = 0; k < jac.outerSize(); ++k) { + for (M::InnerIterator i(jac, k); i ; ++i) { t.push_back(Tri(i.row() + block_row_start, i.col() + block_col_start, i.value())); @@ -592,10 +602,11 @@ vertcatCollapseJacs(const std::vector >& x) } // Build final jacobian. + M comb_jac = M(size, num_cols); + comb_jac.reserve(nnz); + comb_jac.setFromTriplets(t.begin(), t.end()); std::vector jac(1); - jac[0] = Eigen::SparseMatrix(size, num_cols); - jac[0].reserve(nnz); - jac[0].setFromTriplets(t.begin(), t.end()); + jac[0] = ADB::M(comb_jac); // Use move semantics to return result efficiently. return ADB::function(std::move(val), std::move(jac)); diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 5468a4620..1d870fd8f 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -228,8 +228,10 @@ namespace Opm { struct WellOps { WellOps(const Wells* wells); - M w2p; // well -> perf (scatter) - M p2w; // perf -> well (gather) + // M w2p; // well -> perf (scatter) + // M p2w; // perf -> well (gather) + Eigen::SparseMatrix w2p; // well -> perf (scatter) + Eigen::SparseMatrix p2w; // perf -> well (gather) }; // --------- Data members --------- diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index bbbf2592c..7120ef777 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -327,8 +327,8 @@ namespace detail { { if( wells ) { - w2p = M(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); - p2w = M(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); + w2p = Eigen::SparseMatrix(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); + p2w = Eigen::SparseMatrix(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); const int nw = wells->number_of_wells; const int* const wpos = wells->well_connpos; @@ -1441,7 +1441,10 @@ namespace detail { eqs.push_back(residual_.well_eq); ADB total_residual = vertcatCollapseJacs(eqs); const std::vector& Jn = total_residual.derivative(); - const Eigen::SparseLU< M > solver(Jn[0]); + typedef Eigen::SparseMatrix Sp; + Sp Jn0; + Jn[0].toSparse(Jn0); + const Eigen::SparseLU< Sp > solver(Jn0); const Eigen::VectorXd& dx = solver.solve(total_residual.value().matrix()); assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1))); updateWellState(dx.array(), well_state); @@ -1517,7 +1520,7 @@ namespace detail { //Target vars ADB::V bhp_targets = ADB::V::Zero(nw); ADB::V rate_targets = ADB::V::Zero(nw); - ADB::M rate_distr(nw, np*nw); + Eigen::SparseMatrix rate_distr(nw, np*nw); //Selection variables std::vector bhp_elems; @@ -1630,7 +1633,7 @@ namespace detail { // For wells that are dead (not flowing), and therefore not communicating // with the reservoir, we set the equation to be equal to the well's total // flow. This will be a solution only if the target rate is also zero. - M rate_summer(nw, np*nw); + Eigen::SparseMatrix rate_summer(nw, np*nw); for (int w = 0; w < nw; ++w) { for (int phase = 0; phase < np; ++phase) { rate_summer.insert(w, phase*nw + w) = 1.0; @@ -2175,7 +2178,7 @@ namespace detail { const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast(); // Create a sparse vector that nullifies the low potential elements. - const M keep_high_potential = spdiag(high_potential); + const M keep_high_potential(high_potential.matrix().asDiagonal()); // Find the current sign for the threshold modification const V sign_dp = sign(dp.value()); @@ -2641,7 +2644,7 @@ namespace detail { pm[i] = rock_comp_props_->poroMult(p.value()[i]); dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]); } - ADB::M dpm_diag = spdiag(dpm); + ADB::M dpm_diag(dpm.matrix().asDiagonal()); const int num_blocks = p.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { @@ -2669,7 +2672,7 @@ namespace detail { tm[i] = rock_comp_props_->transMult(p.value()[i]); dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]); } - ADB::M dtm_diag = spdiag(dtm); + ADB::M dtm_diag(dtm.matrix().asDiagonal()); const int num_blocks = p.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index 57b3fc5af..7c18577d4 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -384,7 +384,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& if (pw.derivative().empty()) { return ADB::constant(std::move(mu)); } else { - ADB::M dmudp_diag = spdiag(dmudp); + // ADB::M dmudp_diag = spdiag(dmudp); + ADB::M dmudp_diag(dmudp.matrix().asDiagonal()); const int num_blocks = pw.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { @@ -420,15 +421,18 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& props_[phase_usage_.phase_pos[Oil]]->mu(n, pvt_region_.data(), po.value().data(), T.value().data(), rs.value().data(), &cond[0], mu.data(), dmudp.data(), dmudr.data()); - ADB::M dmudp_diag = spdiag(dmudp); - ADB::M dmudr_diag = spdiag(dmudr); + // ADB::M dmudp_diag = spdiag(dmudp); + // ADB::M dmudr_diag = spdiag(dmudr); + ADB::M dmudp_diag(dmudp.matrix().asDiagonal()); + ADB::M dmudr_diag(dmudr.matrix().asDiagonal()); const int num_blocks = po.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { fastSparseProduct(dmudp_diag, po.derivative()[block], jacs[block]); ADB::M temp; fastSparseProduct(dmudr_diag, rs.derivative()[block], temp); - jacs[block] += temp; + // jacs[block] += temp; + jacs[block] = jacs[block] + temp; } return ADB::function(std::move(mu), std::move(jacs)); } @@ -459,15 +463,18 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& props_[phase_usage_.phase_pos[Gas]]->mu(n, pvt_region_.data(), pg.value().data(), T.value().data(), rv.value().data(),&cond[0], mu.data(), dmudp.data(), dmudr.data()); - ADB::M dmudp_diag = spdiag(dmudp); - ADB::M dmudr_diag = spdiag(dmudr); + // ADB::M dmudp_diag = spdiag(dmudp); + // ADB::M dmudr_diag = spdiag(dmudr); + ADB::M dmudp_diag(dmudp.matrix().asDiagonal()); + ADB::M dmudr_diag(dmudr.matrix().asDiagonal()); const int num_blocks = pg.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { fastSparseProduct(dmudp_diag, pg.derivative()[block], jacs[block]); ADB::M temp; fastSparseProduct(dmudr_diag, rv.derivative()[block], temp); - jacs[block] += temp; + // jacs[block] += temp; + jacs[block] = jacs[block] + temp; } return ADB::function(std::move(mu), std::move(jacs)); } @@ -500,7 +507,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& props_[phase_usage_.phase_pos[Water]]->b(n, pvt_region_.data(), pw.value().data(), T.value().data(), rs, b.data(), dbdp.data(), dbdr.data()); - ADB::M dbdp_diag = spdiag(dbdp); + // ADB::M dbdp_diag = spdiag(dbdp); + ADB::M dbdp_diag(dbdp.matrix().asDiagonal()); const int num_blocks = pw.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { @@ -536,15 +544,18 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& props_[phase_usage_.phase_pos[Oil]]->b(n, pvt_region_.data(), po.value().data(), T.value().data(), rs.value().data(), &cond[0], b.data(), dbdp.data(), dbdr.data()); - ADB::M dbdp_diag = spdiag(dbdp); - ADB::M dbdr_diag = spdiag(dbdr); + // ADB::M dbdp_diag = spdiag(dbdp); + // ADB::M dbdr_diag = spdiag(dbdr); + ADB::M dbdp_diag(dbdp.matrix().asDiagonal()); + ADB::M dbdr_diag(dbdr.matrix().asDiagonal()); const int num_blocks = po.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { fastSparseProduct(dbdp_diag, po.derivative()[block], jacs[block]); ADB::M temp; fastSparseProduct(dbdr_diag, rs.derivative()[block], temp); - jacs[block] += temp; + // jacs[block] += temp; + jacs[block] = jacs[block] + temp; } return ADB::function(std::move(b), std::move(jacs)); } @@ -576,15 +587,18 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& props_[phase_usage_.phase_pos[Gas]]->b(n, pvt_region_.data(), pg.value().data(), T.value().data(), rv.value().data(), &cond[0], b.data(), dbdp.data(), dbdr.data()); - ADB::M dbdp_diag = spdiag(dbdp); - ADB::M dbdr_diag = spdiag(dbdr); + // ADB::M dbdp_diag = spdiag(dbdp); + // ADB::M dbdr_diag = spdiag(dbdr); + ADB::M dbdp_diag(dbdp.matrix().asDiagonal()); + ADB::M dbdr_diag(dbdr.matrix().asDiagonal()); const int num_blocks = pg.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]); ADB::M temp; fastSparseProduct(dbdr_diag, rv.derivative()[block], temp); - jacs[block] += temp; + // jacs[block] += temp; + jacs[block] = jacs[block] + temp; } return ADB::function(std::move(b), std::move(jacs)); } @@ -609,7 +623,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& V rbub(n); V drbubdp(n); props_[phase_usage_.phase_pos[Oil]]->rsSat(n, pvt_region_.data(), po.value().data(), rbub.data(), drbubdp.data()); - ADB::M drbubdp_diag = spdiag(drbubdp); + // ADB::M drbubdp_diag = spdiag(drbubdp); + ADB::M drbubdp_diag(drbubdp.matrix().asDiagonal()); const int num_blocks = po.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { @@ -650,7 +665,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& V rv(n); V drvdp(n); props_[phase_usage_.phase_pos[Gas]]->rvSat(n, pvt_region_.data(), po.value().data(), rv.data(), drvdp.data()); - ADB::M drvdp_diag = spdiag(drvdp); + // ADB::M drvdp_diag = spdiag(drvdp); + ADB::M drvdp_diag(drvdp.matrix().asDiagonal()); const int num_blocks = po.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { @@ -726,11 +742,13 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& const int phase2_pos = phase_usage_.phase_pos[phase2]; // Assemble dkr1/ds2. const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm() - ADB::M dkr1_ds2_diag = spdiag(dkr.col(column)); + // ADB::M dkr1_ds2_diag = spdiag(dkr.col(column)); + ADB::M dkr1_ds2_diag(dkr.col(column).matrix().asDiagonal()); for (int block = 0; block < num_blocks; ++block) { ADB::M temp; fastSparseProduct(dkr1_ds2_diag, s[phase2]->derivative()[block], temp); - jacs[block] += temp; + // jacs[block] += temp; + jacs[block] = jacs[block] + temp; } } ADB::V val = kr.col(phase1_pos); @@ -787,11 +805,13 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& const int phase2_pos = phase_usage_.phase_pos[phase2]; // Assemble dpc1/ds2. const int column = phase1_pos + nActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm() - ADB::M dpc1_ds2_diag = spdiag(dpc.col(column)); + // ADB::M dpc1_ds2_diag = spdiag(dpc.col(column)); + ADB::M dpc1_ds2_diag(dpc.col(column).matrix().asDiagonal()); for (int block = 0; block < nBlocks; ++block) { ADB::M temp; fastSparseProduct(dpc1_ds2_diag, s[phase2]->derivative()[block], temp); - jacs[block] += temp; + // jacs[block] += temp; + jacs[block] = jacs[block] + temp; } } ADB::V val = pc.col(phase1_pos); @@ -892,7 +912,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck& dfactor_dso[i] = vap*std::pow(so_i/satOilMax_[cells[i]], vap-1.0)/satOilMax_[cells[i]]; } } - ADB::M dfactor_dso_diag = spdiag(dfactor_dso); + // ADB::M dfactor_dso_diag = spdiag(dfactor_dso); + ADB::M dfactor_dso_diag(dfactor_dso.matrix().asDiagonal()); const int num_blocks = so.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp index c4906abca..163069757 100644 --- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp +++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp @@ -208,9 +208,13 @@ namespace Opm // Find sparsity structure as union of basic block sparsity structures, // corresponding to the jacobians with respect to pressure. // Use addition to get to the union structure. - Eigen::SparseMatrix structure = eqs[0].derivative()[0]; + typedef Eigen::SparseMatrix Sp; + Sp structure; + eqs[0].derivative()[0].toSparse(structure); for (int phase = 1; phase < np; ++phase) { - structure += eqs[phase].derivative()[0]; + Sp s0; + eqs[phase].derivative()[0].toSparse(s0); + structure += s0; } Eigen::SparseMatrix s = structure; diff --git a/opm/autodiff/NewtonIterationBlackoilSimple.cpp b/opm/autodiff/NewtonIterationBlackoilSimple.cpp index b7a0d89f4..0adb17679 100644 --- a/opm/autodiff/NewtonIterationBlackoilSimple.cpp +++ b/opm/autodiff/NewtonIterationBlackoilSimple.cpp @@ -47,6 +47,7 @@ namespace Opm NewtonIterationBlackoilSimple::SolutionVector NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const { + /* typedef LinearisedBlackoilResidual::ADB ADB; const int np = residual.material_balance_eq.size(); ADB mass_res = residual.material_balance_eq[0]; @@ -73,6 +74,7 @@ namespace Opm "Linear solver convergence failure."); } return dx; + */ } const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const diff --git a/opm/autodiff/NewtonIterationUtilities.cpp b/opm/autodiff/NewtonIterationUtilities.cpp index 71c389280..0d813ef7f 100644 --- a/opm/autodiff/NewtonIterationUtilities.cpp +++ b/opm/autodiff/NewtonIterationUtilities.cpp @@ -62,14 +62,17 @@ namespace Opm const std::vector& Jn = eqs[n].derivative(); // Use sparse LU to solve the block submatrices i.e compute inv(D) + typedef Eigen::SparseMatrix Sp; + Sp Jnn; + Jn[n].toSparse(Jnn); #if HAVE_UMFPACK - const Eigen::UmfPackLU< M > solver(Jn[n]); + const Eigen::UmfPackLU solver(Jnn); #else - const Eigen::SparseLU< M > solver(Jn[n]); + const Eigen::SparseLU solver(Jnn); #endif - M id(Jn[n].rows(), Jn[n].cols()); + Sp id(Jn[n].rows(), Jn[n].cols()); id.setIdentity(); - const Eigen::SparseMatrix Di = solver.solve(id); + const Sp Di = solver.solve(id); // compute inv(D)*bn for the update of the right hand side const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix()); @@ -105,7 +108,9 @@ namespace Opm // Subtract Bu (B*inv(D)*C) M Bu; fastSparseProduct(B, u, Bu); - J -= Bu; + // J -= Bu; + Bu = Bu * -1.0; + J = J + Bu; } } @@ -136,7 +141,7 @@ namespace Opm // the eliminated equation, and x is the partial solution // of the non-eliminated unknowns. - const M& D = equation.derivative()[n]; + const M& D1 = equation.derivative()[n]; // Build C. std::vector C_jacs = equation.derivative(); C_jacs.erase(C_jacs.begin() + n); @@ -145,10 +150,13 @@ namespace Opm const M& C = eq_coll.derivative()[0]; // Use sparse LU to solve the block submatrices + typedef Eigen::SparseMatrix Sp; + Sp D; + D1.toSparse(D); #if HAVE_UMFPACK - const Eigen::UmfPackLU< M > solver(D); + const Eigen::UmfPackLU solver(D); #else - const Eigen::SparseLU< M > solver(D); + const Eigen::SparseLU solver(D); #endif // Compute value of eliminated variable. @@ -190,6 +198,7 @@ namespace Opm Eigen::SparseMatrix& A, V& b) { + /* if (num_phases != 3) { OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases."); } @@ -269,6 +278,7 @@ namespace Opm // Create output as product of L with equations. A = L * total_residual.derivative()[0]; b = L * total_residual.value().matrix(); + */ } diff --git a/opm/autodiff/VFPInjProperties.cpp b/opm/autodiff/VFPInjProperties.cpp index fb1345ce8..5300a0199 100644 --- a/opm/autodiff/VFPInjProperties.cpp +++ b/opm/autodiff/VFPInjProperties.cpp @@ -144,8 +144,10 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, } //Create diagonal matrices from ADB::Vs - ADB::M dthp_diag = spdiag(dthp); - ADB::M dflo_diag = spdiag(dflo); + // ADB::M dthp_diag = spdiag(dthp); + // ADB::M dflo_diag = spdiag(dflo); + ADB::M dthp_diag(dthp.matrix().asDiagonal()); + ADB::M dflo_diag(dflo.matrix().asDiagonal()); //Calculate the Jacobians const int num_blocks = block_pattern.size(); diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index f69280b29..d2f667500 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -143,11 +143,16 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, } //Create diagonal matrices from ADB::Vs - ADB::M dthp_diag = spdiag(dthp); - ADB::M dwfr_diag = spdiag(dwfr); - ADB::M dgfr_diag = spdiag(dgfr); - ADB::M dalq_diag = spdiag(dalq); - ADB::M dflo_diag = spdiag(dflo); + // ADB::M dthp_diag = spdiag(dthp); + // ADB::M dwfr_diag = spdiag(dwfr); + // ADB::M dgfr_diag = spdiag(dgfr); + // ADB::M dalq_diag = spdiag(dalq); + // ADB::M dflo_diag = spdiag(dflo); + ADB::M dthp_diag(dthp.matrix().asDiagonal()); + ADB::M dwfr_diag(dwfr.matrix().asDiagonal()); + ADB::M dgfr_diag(dgfr.matrix().asDiagonal()); + ADB::M dalq_diag(dalq.matrix().asDiagonal()); + ADB::M dflo_diag(dflo.matrix().asDiagonal()); //Calculate the Jacobians const int num_blocks = block_pattern.size();