OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Changes to make flow compile with AutoDiffMatrix.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2015-08-25 16:15:29 +02:00 committed by babrodtk
parent c712d0070d
commit c795113ce3
11 changed files with 195 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
CMakeLists_files.cmake
View File

@ -33,12 +33,12 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/NewtonIterationBlackoilSimple.cpp opm/autodiff/NewtonIterationBlackoilSimple.cpp
opm/autodiff/NewtonIterationUtilities.cpp opm/autodiff/NewtonIterationUtilities.cpp
opm/autodiff/GridHelpers.cpp opm/autodiff/GridHelpers.cpp
opm/autodiff/ImpesTPFAAD.cpp # opm/autodiff/ImpesTPFAAD.cpp
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
opm/autodiff/SimulatorIncompTwophaseAd.cpp # opm/autodiff/SimulatorIncompTwophaseAd.cpp
opm/autodiff/TransportSolverTwophaseAd.cpp # opm/autodiff/TransportSolverTwophaseAd.cpp
opm/autodiff/BlackoilPropsAdFromDeck.cpp opm/autodiff/BlackoilPropsAdFromDeck.cpp
opm/autodiff/SolventPropsAdFromDeck.cpp # opm/autodiff/SolventPropsAdFromDeck.cpp
opm/autodiff/BlackoilModelParameters.cpp opm/autodiff/BlackoilModelParameters.cpp
opm/autodiff/WellDensitySegmented.cpp opm/autodiff/WellDensitySegmented.cpp
opm/autodiff/LinearisedBlackoilResidual.cpp opm/autodiff/LinearisedBlackoilResidual.cpp
@ -83,20 +83,20 @@ endif()
# originally generated with the command: # originally generated with the command:
# find tutorials examples -name '*.c*' -printf '\t%p\n' | sort # find tutorials examples -name '*.c*' -printf '\t%p\n' | sort
list (APPEND EXAMPLE_SOURCE_FILES list (APPEND EXAMPLE_SOURCE_FILES
examples/find_zero.cpp # examples/find_zero.cpp
examples/flow.cpp examples/flow.cpp
examples/flow_solvent.cpp # examples/flow_solvent.cpp
examples/sim_2p_incomp_ad.cpp # examples/sim_2p_incomp_ad.cpp
examples/sim_simple.cpp # examples/sim_simple.cpp
examples/opm_init_check.cpp # examples/opm_init_check.cpp
) )
# programs listed here will not only be compiled, but also marked for # programs listed here will not only be compiled, but also marked for
# installation # installation
list (APPEND PROGRAM_SOURCE_FILES list (APPEND PROGRAM_SOURCE_FILES
examples/sim_2p_incomp_ad.cpp # examples/sim_2p_incomp_ad.cpp
examples/flow.cpp examples/flow.cpp
examples/opm_init_check.cpp # examples/opm_init_check.cpp
) )
# originally generated with the command: # originally generated with the command:

82
opm/autodiff/AutoDiffBlock.hpp
View File

@ -28,6 +28,9 @@
#include <opm/core/utility/platform_dependent/reenable_warnings.h> #include <opm/core/utility/platform_dependent/reenable_warnings.h>
#include <opm/autodiff/AutoDiffMatrix.hpp>
#include <utility> #include <utility>
#include <vector> #include <vector>
#include <cassert> #include <cassert>
@ -94,7 +97,7 @@ namespace Opm
/// Underlying type for values. /// Underlying type for values.
typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V; typedef Eigen::Array<Scalar, Eigen::Dynamic, 1> V;
/// Underlying type for jacobians. /// Underlying type for jacobians.
typedef Eigen::SparseMatrix<Scalar> M; typedef AutoDiffMatrix M;
/// Construct an empty AutoDiffBlock. /// Construct an empty AutoDiffBlock.
static AutoDiffBlock null() static AutoDiffBlock null()
@ -155,22 +158,23 @@ namespace Opm
} }
// ... then set the one corrresponding to this variable to identity. // ... then set the one corrresponding to this variable to identity.
assert(blocksizes[index] == num_elem); assert(blocksizes[index] == num_elem);
if(val.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 // // if val is empty the following will run into an assertion
jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1)); // // with Eigen
for (typename M::Index row = 0; row < val.size(); ++row) { // jac[index].reserve(Eigen::VectorXi::Constant(val.size(), 1));
jac[index].insert(row, row) = Scalar(1.0); // for (typename M::Index row = 0; row < val.size(); ++row) {
} // jac[index].insert(row, row) = Scalar(1.0);
} // }
else // }
{ // else
// Do some check. Empty jacobian only make sense for empty val. // {
for (int i = 0; i < num_blocks; ++i) { // // Do some check. Empty jacobian only make sense for empty val.
assert(jac[i].size()==0); // for (int i = 0; i < num_blocks; ++i) {
} // assert(jac[i].size()==0);
} // }
// }
return AutoDiffBlock(std::move(val), std::move(jac)); return AutoDiffBlock(std::move(val), std::move(jac));
} }
@ -254,7 +258,8 @@ namespace Opm
const int num_blocks = rhs.numBlocks(); const int num_blocks = rhs.numBlocks();
jac_.resize(num_blocks); jac_.resize(num_blocks);
for (int block = 0; block < num_blocks; ++block) { 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()) { } else if (!rhs.jac_.empty()) {
assert (numBlocks() == rhs.numBlocks()); assert (numBlocks() == rhs.numBlocks());
@ -334,9 +339,11 @@ namespace Opm
int num_blocks = numBlocks(); int num_blocks = numBlocks();
std::vector<M> jac(num_blocks); std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks()); assert(numBlocks() == rhs.numBlocks());
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D; // typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal(); // D D1 = val_.matrix().asDiagonal();
D D2 = rhs.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) { for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows()); assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols()); assert(jac_[block].cols() == rhs.jac_[block].cols());
@ -371,9 +378,12 @@ namespace Opm
std::vector<M> jac(num_blocks); std::vector<M> jac(num_blocks);
assert(numBlocks() == rhs.numBlocks()); assert(numBlocks() == rhs.numBlocks());
typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D; typedef Eigen::DiagonalMatrix<Scalar, Eigen::Dynamic> D;
D D1 = val_.matrix().asDiagonal(); // D D1 = val_.matrix().asDiagonal();
D D2 = rhs.val_.matrix().asDiagonal(); // D D2 = rhs.val_.matrix().asDiagonal();
D D3 = (1.0/(rhs.val_*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) { for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows()); assert(jac_[block].rows() == rhs.jac_[block].rows());
assert(jac_[block].cols() == rhs.jac_[block].cols()); assert(jac_[block].cols() == rhs.jac_[block].cols());
@ -381,15 +391,14 @@ namespace Opm
jac[block] = M( D3.rows(), jac_[block].cols() ); jac[block] = M( D3.rows(), jac_[block].cols() );
} }
else if( jac_[block].nonZeros() == 0 ) { else if( jac_[block].nonZeros() == 0 ) {
jac[block] = D3 * ( D1*rhs.jac_[block]); jac[block] = D3 * ( D1*rhs.jac_[block]) * (-1.0);
jac[block] *= -1.0;
} }
else if ( rhs.jac_[block].nonZeros() == 0 ) else if ( rhs.jac_[block].nonZeros() == 0 )
{ {
jac[block] = D3 * (D2*jac_[block]); jac[block] = D3 * (D2*jac_[block]);
} }
else { 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)); 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 <typename Scalar> template <typename Scalar>
AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::M& lhs, AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::M& lhs,
const AutoDiffBlock<Scalar>& rhs) const AutoDiffBlock<Scalar>& rhs)
@ -506,6 +515,23 @@ namespace Opm
} }
/// Multiply with Eigen sparse matrix from the left.
template <typename Scalar>
AutoDiffBlock<Scalar> operator*(const Eigen::SparseMatrix<Scalar>& lhs,
const AutoDiffBlock<Scalar>& rhs)
{
int num_blocks = rhs.numBlocks();
std::vector<typename AutoDiffBlock<Scalar>::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<Scalar>::V val = lhs*rhs.value().matrix();
return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));
}
/// Elementwise multiplication with constant on the left. /// Elementwise multiplication with constant on the left.
template <typename Scalar> template <typename Scalar>
AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::V& lhs, AutoDiffBlock<Scalar> operator*(const typename AutoDiffBlock<Scalar>::V& lhs,

55
opm/autodiff/AutoDiffHelpers.hpp
View File

@ -39,7 +39,8 @@ namespace Opm
/// operations on (AD or regular) vectors of data. /// operations on (AD or regular) vectors of data.
struct HelperOps struct HelperOps
{ {
typedef AutoDiffBlock<double>::M M; // typedef AutoDiffBlock<double>::M M;
typedef Eigen::SparseMatrix<double> M;
typedef AutoDiffBlock<double>::V V; typedef AutoDiffBlock<double>::V V;
/// A list of internal faces. /// A list of internal faces.
@ -250,7 +251,7 @@ struct HelperOps
} }
private: private:
typename ADB::M select_; Eigen::SparseMatrix<double> select_;
}; };
@ -259,17 +260,17 @@ namespace {
template <typename Scalar, class IntVec> template <typename Scalar, class IntVec>
typename AutoDiffBlock<Scalar>::M typename Eigen::SparseMatrix<Scalar>
constructSupersetSparseMatrix(const int full_size, const IntVec& indices) constructSupersetSparseMatrix(const int full_size, const IntVec& indices)
{ {
const int subset_size = indices.size(); const int subset_size = indices.size();
if (subset_size == 0) { if (subset_size == 0) {
typename AutoDiffBlock<Scalar>::M mat(full_size, 0); Eigen::SparseMatrix<Scalar> mat(full_size, 0);
return mat; return mat;
} }
typename AutoDiffBlock<Scalar>::M mat(full_size, subset_size); Eigen::SparseMatrix<Scalar> mat(full_size, subset_size);
mat.reserve(Eigen::VectorXi::Constant(subset_size, 1)); mat.reserve(Eigen::VectorXi::Constant(subset_size, 1));
for (int i = 0; i < subset_size; ++i) { for (int i = 0; i < subset_size; ++i) {
mat.insert(indices[i], i) = 1; mat.insert(indices[i], i) = 1;
@ -302,8 +303,8 @@ AutoDiffBlock<Scalar>
subset(const AutoDiffBlock<Scalar>& x, subset(const AutoDiffBlock<Scalar>& x,
const IntVec& indices) const IntVec& indices)
{ {
const typename AutoDiffBlock<Scalar>::M sub const Eigen::SparseMatrix<Scalar> sub
= constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose(); = constructSupersetSparseMatrix<Scalar>(x.value().size(), indices).transpose().eval();
return sub * x; return sub * x;
} }
@ -336,10 +337,12 @@ superset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
/// elements of d on the diagonal. /// elements of d on the diagonal.
/// Need to mark this as inline since it is defined in a header and not a template. /// Need to mark this as inline since it is defined in a header and not a template.
inline inline
AutoDiffBlock<double>::M //AutoDiffBlock<double>::M
Eigen::SparseMatrix<double>
spdiag(const AutoDiffBlock<double>::V& d) spdiag(const AutoDiffBlock<double>::V& d)
{ {
typedef AutoDiffBlock<double>::M M; // typedef AutoDiffBlock<double>::M M;
typedef Eigen::SparseMatrix<double> M;
const int n = d.size(); const int n = d.size();
M mat(n, n); M mat(n, n);
@ -456,9 +459,11 @@ collapseJacs(const AutoDiffBlock<double>& x, Matrix& jacobian)
t.reserve(nnz); t.reserve(nnz);
int block_col_start = 0; int block_col_start = 0;
for (int block = 0; block < nb; ++block) { for (int block = 0; block < nb; ++block) {
const ADB::M& jac = x.derivative()[block]; const ADB::M& jac1 = x.derivative()[block];
for (ADB::M::Index k = 0; k < jac.outerSize(); ++k) { Eigen::SparseMatrix<double> jac;
for (ADB::M::InnerIterator i(jac, k); i ; ++i) { jac1.toSparse(jac);
for (Eigen::SparseMatrix<double>::Index k = 0; k < jac.outerSize(); ++k) {
for (Eigen::SparseMatrix<double>::InnerIterator i(jac, k); i ; ++i) {
t.push_back(Tri(i.row(), t.push_back(Tri(i.row(),
i.col() + block_col_start, i.col() + block_col_start,
i.value())); i.value()));
@ -478,16 +483,18 @@ inline
AutoDiffBlock<double> AutoDiffBlock<double>
collapseJacs(const AutoDiffBlock<double>& x) collapseJacs(const AutoDiffBlock<double>& x)
{ {
typedef AutoDiffBlock<double> ADB; Eigen::SparseMatrix<double> comb_jac;
collapseJacs( x, comb_jac );
// Build final jacobian. // Build final jacobian.
typedef AutoDiffBlock<double> ADB;
std::vector<ADB::M> jacs(1); std::vector<ADB::M> jacs(1);
collapseJacs( x, jacs[ 0 ] ); jacs[0] = AutoDiffMatrix(std::move(comb_jac));
ADB::V val = x.value(); ADB::V val = x.value();
return ADB::function(std::move(val), std::move(jacs)); return ADB::function(std::move(val), std::move(jacs));
} }
/*
/// Returns the vertical concatenation [ x; y ] of the inputs. /// Returns the vertical concatenation [ x; y ] of the inputs.
inline inline
@ -510,7 +517,7 @@ vertcat(const AutoDiffBlock<double>& x,
} }
*/
/// Returns the vertical concatenation [ x[0]; x[1]; ...; x[n-1] ] of the inputs. /// Returns the vertical concatenation [ x[0]; x[1]; ...; x[n-1] ] of the inputs.
@ -569,6 +576,7 @@ vertcatCollapseJacs(const std::vector<AutoDiffBlock<double> >& x)
} }
// Set up for batch insertion of all Jacobian elements. // Set up for batch insertion of all Jacobian elements.
typedef Eigen::SparseMatrix<double> M;
typedef Eigen::Triplet<double> Tri; typedef Eigen::Triplet<double> Tri;
std::vector<Tri> t; std::vector<Tri> t;
t.reserve(nnz); t.reserve(nnz);
@ -577,9 +585,11 @@ vertcatCollapseJacs(const std::vector<AutoDiffBlock<double> >& x)
int block_col_start = 0; int block_col_start = 0;
if (!x[elem].derivative().empty()) { if (!x[elem].derivative().empty()) {
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
const ADB::M& jac = x[elem].derivative()[block]; // const ADB::M& jac = x[elem].derivative()[block];
for (ADB::M::Index k = 0; k < jac.outerSize(); ++k) { M jac;
for (ADB::M::InnerIterator i(jac, k); i ; ++i) { 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, t.push_back(Tri(i.row() + block_row_start,
i.col() + block_col_start, i.col() + block_col_start,
i.value())); i.value()));
@ -592,10 +602,11 @@ vertcatCollapseJacs(const std::vector<AutoDiffBlock<double> >& x)
} }
// Build final jacobian. // Build final jacobian.
M comb_jac = M(size, num_cols);
comb_jac.reserve(nnz);
comb_jac.setFromTriplets(t.begin(), t.end());
std::vector<ADB::M> jac(1); std::vector<ADB::M> jac(1);
jac[0] = Eigen::SparseMatrix<double>(size, num_cols); jac[0] = ADB::M(comb_jac);
jac[0].reserve(nnz);
jac[0].setFromTriplets(t.begin(), t.end());
// Use move semantics to return result efficiently. // Use move semantics to return result efficiently.
return ADB::function(std::move(val), std::move(jac)); return ADB::function(std::move(val), std::move(jac));

6
opm/autodiff/BlackoilModelBase.hpp
View File

@ -228,8 +228,10 @@ namespace Opm {
struct WellOps { struct WellOps {
WellOps(const Wells* wells); WellOps(const Wells* wells);
M w2p; // well -> perf (scatter) // M w2p; // well -> perf (scatter)
M p2w; // perf -> well (gather) // M p2w; // perf -> well (gather)
Eigen::SparseMatrix<double> w2p; // well -> perf (scatter)
Eigen::SparseMatrix<double> p2w; // perf -> well (gather)
}; };
// --------- Data members --------- // --------- Data members ---------

19
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -327,8 +327,8 @@ namespace detail {
{ {
if( wells ) if( wells )
{ {
w2p = M(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); w2p = Eigen::SparseMatrix<double>(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells);
p2w = M(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); p2w = Eigen::SparseMatrix<double>(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]);
const int nw = wells->number_of_wells; const int nw = wells->number_of_wells;
const int* const wpos = wells->well_connpos; const int* const wpos = wells->well_connpos;
@ -1441,7 +1441,10 @@ namespace detail {
eqs.push_back(residual_.well_eq); eqs.push_back(residual_.well_eq);
ADB total_residual = vertcatCollapseJacs(eqs); ADB total_residual = vertcatCollapseJacs(eqs);
const std::vector<M>& Jn = total_residual.derivative(); const std::vector<M>& Jn = total_residual.derivative();
const Eigen::SparseLU< M > solver(Jn[0]); typedef Eigen::SparseMatrix<double> Sp;
Sp Jn0;
Jn[0].toSparse(Jn0);
const Eigen::SparseLU< Sp > solver(Jn0);
const Eigen::VectorXd& dx = solver.solve(total_residual.value().matrix()); const Eigen::VectorXd& dx = solver.solve(total_residual.value().matrix());
assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1))); assert(dx.size() == (well_state.numWells() * (well_state.numPhases()+1)));
updateWellState(dx.array(), well_state); updateWellState(dx.array(), well_state);
@ -1517,7 +1520,7 @@ namespace detail {
//Target vars //Target vars
ADB::V bhp_targets = ADB::V::Zero(nw); ADB::V bhp_targets = ADB::V::Zero(nw);
ADB::V rate_targets = ADB::V::Zero(nw); ADB::V rate_targets = ADB::V::Zero(nw);
ADB::M rate_distr(nw, np*nw); Eigen::SparseMatrix<double> rate_distr(nw, np*nw);
//Selection variables //Selection variables
std::vector<int> bhp_elems; std::vector<int> bhp_elems;
@ -1630,7 +1633,7 @@ namespace detail {
// For wells that are dead (not flowing), and therefore not communicating // 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 // 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. // flow. This will be a solution only if the target rate is also zero.
M rate_summer(nw, np*nw); Eigen::SparseMatrix<double> rate_summer(nw, np*nw);
for (int w = 0; w < nw; ++w) { for (int w = 0; w < nw; ++w) {
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
rate_summer.insert(w, phase*nw + w) = 1.0; 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<double>(); const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast<double>();
// Create a sparse vector that nullifies the low potential elements. // 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 // Find the current sign for the threshold modification
const V sign_dp = sign(dp.value()); const V sign_dp = sign(dp.value());
@ -2641,7 +2644,7 @@ namespace detail {
pm[i] = rock_comp_props_->poroMult(p.value()[i]); pm[i] = rock_comp_props_->poroMult(p.value()[i]);
dpm[i] = rock_comp_props_->poroMultDeriv(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(); const int num_blocks = p.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
@ -2669,7 +2672,7 @@ namespace detail {
tm[i] = rock_comp_props_->transMult(p.value()[i]); tm[i] = rock_comp_props_->transMult(p.value()[i]);
dtm[i] = rock_comp_props_->transMultDeriv(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(); const int num_blocks = p.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {

63
opm/autodiff/BlackoilPropsAdFromDeck.cpp
View File

@ -384,7 +384,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
if (pw.derivative().empty()) { if (pw.derivative().empty()) {
return ADB::constant(std::move(mu)); return ADB::constant(std::move(mu));
} else { } 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(); const int num_blocks = pw.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { 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(), 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()); &cond[0], mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp); // ADB::M dmudp_diag = spdiag(dmudp);
ADB::M dmudr_diag = spdiag(dmudr); // 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(); const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dmudp_diag, po.derivative()[block], jacs[block]); fastSparseProduct(dmudp_diag, po.derivative()[block], jacs[block]);
ADB::M temp; ADB::M temp;
fastSparseProduct(dmudr_diag, rs.derivative()[block], 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)); 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], 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()); mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp); // ADB::M dmudp_diag = spdiag(dmudp);
ADB::M dmudr_diag = spdiag(dmudr); // 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(); const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dmudp_diag, pg.derivative()[block], jacs[block]); fastSparseProduct(dmudp_diag, pg.derivative()[block], jacs[block]);
ADB::M temp; ADB::M temp;
fastSparseProduct(dmudr_diag, rv.derivative()[block], 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)); 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, props_[phase_usage_.phase_pos[Water]]->b(n, pvt_region_.data(), pw.value().data(), T.value().data(), rs,
b.data(), dbdp.data(), dbdr.data()); 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(); const int num_blocks = pw.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { 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(), 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()); &cond[0], b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp); // ADB::M dbdp_diag = spdiag(dbdp);
ADB::M dbdr_diag = spdiag(dbdr); // 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(); const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dbdp_diag, po.derivative()[block], jacs[block]); fastSparseProduct(dbdp_diag, po.derivative()[block], jacs[block]);
ADB::M temp; ADB::M temp;
fastSparseProduct(dbdr_diag, rs.derivative()[block], 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)); 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], 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()); b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp); // ADB::M dbdp_diag = spdiag(dbdp);
ADB::M dbdr_diag = spdiag(dbdr); // 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(); const int num_blocks = pg.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]); fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]);
ADB::M temp; ADB::M temp;
fastSparseProduct(dbdr_diag, rv.derivative()[block], 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)); return ADB::function(std::move(b), std::move(jacs));
} }
@ -609,7 +623,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
V rbub(n); V rbub(n);
V drbubdp(n); V drbubdp(n);
props_[phase_usage_.phase_pos[Oil]]->rsSat(n, pvt_region_.data(), po.value().data(), rbub.data(), drbubdp.data()); 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(); const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
@ -650,7 +665,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
V rv(n); V rv(n);
V drvdp(n); V drvdp(n);
props_[phase_usage_.phase_pos[Gas]]->rvSat(n, pvt_region_.data(), po.value().data(), rv.data(), drvdp.data()); 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(); const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { 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]; const int phase2_pos = phase_usage_.phase_pos[phase2];
// Assemble dkr1/ds2. // Assemble dkr1/ds2.
const int column = phase1_pos + np*phase2_pos; // Recall: Fortran ordering from props_.relperm() 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) { for (int block = 0; block < num_blocks; ++block) {
ADB::M temp; ADB::M temp;
fastSparseProduct(dkr1_ds2_diag, s[phase2]->derivative()[block], 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); ADB::V val = kr.col(phase1_pos);
@ -787,11 +805,13 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
const int phase2_pos = phase_usage_.phase_pos[phase2]; const int phase2_pos = phase_usage_.phase_pos[phase2];
// Assemble dpc1/ds2. // Assemble dpc1/ds2.
const int column = phase1_pos + nActivePhases*phase2_pos; // Recall: Fortran ordering from props_.relperm() 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) { for (int block = 0; block < nBlocks; ++block) {
ADB::M temp; ADB::M temp;
fastSparseProduct(dpc1_ds2_diag, s[phase2]->derivative()[block], 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); 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]]; 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(); const int num_blocks = so.numBlocks();
std::vector<ADB::M> jacs(num_blocks); std::vector<ADB::M> jacs(num_blocks);
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {

8
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@ -208,9 +208,13 @@ namespace Opm
// Find sparsity structure as union of basic block sparsity structures, // Find sparsity structure as union of basic block sparsity structures,
// corresponding to the jacobians with respect to pressure. // corresponding to the jacobians with respect to pressure.
// Use addition to get to the union structure. // Use addition to get to the union structure.
Eigen::SparseMatrix<double> structure = eqs[0].derivative()[0]; typedef Eigen::SparseMatrix<double> Sp;
Sp structure;
eqs[0].derivative()[0].toSparse(structure);
for (int phase = 1; phase < np; ++phase) { for (int phase = 1; phase < np; ++phase) {
structure += eqs[phase].derivative()[0]; Sp s0;
eqs[phase].derivative()[0].toSparse(s0);
structure += s0;
} }
Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure; Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure;

2
opm/autodiff/NewtonIterationBlackoilSimple.cpp
View File

@ -47,6 +47,7 @@ namespace Opm
NewtonIterationBlackoilSimple::SolutionVector NewtonIterationBlackoilSimple::SolutionVector
NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
{ {
/*
typedef LinearisedBlackoilResidual::ADB ADB; typedef LinearisedBlackoilResidual::ADB ADB;
const int np = residual.material_balance_eq.size(); const int np = residual.material_balance_eq.size();
ADB mass_res = residual.material_balance_eq[0]; ADB mass_res = residual.material_balance_eq[0];
@ -73,6 +74,7 @@ namespace Opm
"Linear solver convergence failure."); "Linear solver convergence failure.");
} }
return dx; return dx;
*/
} }
const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const const boost::any& NewtonIterationBlackoilSimple::parallelInformation() const

26
opm/autodiff/NewtonIterationUtilities.cpp
View File

@ -62,14 +62,17 @@ namespace Opm
const std::vector<M>& Jn = eqs[n].derivative(); const std::vector<M>& Jn = eqs[n].derivative();
// Use sparse LU to solve the block submatrices i.e compute inv(D) // Use sparse LU to solve the block submatrices i.e compute inv(D)
typedef Eigen::SparseMatrix<double> Sp;
Sp Jnn;
Jn[n].toSparse(Jnn);
#if HAVE_UMFPACK #if HAVE_UMFPACK
const Eigen::UmfPackLU< M > solver(Jn[n]); const Eigen::UmfPackLU<Sp> solver(Jnn);
#else #else
const Eigen::SparseLU< M > solver(Jn[n]); const Eigen::SparseLU<Sp> solver(Jnn);
#endif #endif
M id(Jn[n].rows(), Jn[n].cols()); Sp id(Jn[n].rows(), Jn[n].cols());
id.setIdentity(); id.setIdentity();
const Eigen::SparseMatrix<M::Scalar, Eigen::ColMajor> Di = solver.solve(id); const Sp Di = solver.solve(id);
// compute inv(D)*bn for the update of the right hand side // compute inv(D)*bn for the update of the right hand side
const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix()); const Eigen::VectorXd& Dibn = solver.solve(eqs[n].value().matrix());
@ -105,7 +108,9 @@ namespace Opm
// Subtract Bu (B*inv(D)*C) // Subtract Bu (B*inv(D)*C)
M Bu; M Bu;
fastSparseProduct(B, u, 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 // the eliminated equation, and x is the partial solution
// of the non-eliminated unknowns. // of the non-eliminated unknowns.
const M& D = equation.derivative()[n]; const M& D1 = equation.derivative()[n];
// Build C. // Build C.
std::vector<M> C_jacs = equation.derivative(); std::vector<M> C_jacs = equation.derivative();
C_jacs.erase(C_jacs.begin() + n); C_jacs.erase(C_jacs.begin() + n);
@ -145,10 +150,13 @@ namespace Opm
const M& C = eq_coll.derivative()[0]; const M& C = eq_coll.derivative()[0];
// Use sparse LU to solve the block submatrices // Use sparse LU to solve the block submatrices
typedef Eigen::SparseMatrix<double> Sp;
Sp D;
D1.toSparse(D);
#if HAVE_UMFPACK #if HAVE_UMFPACK
const Eigen::UmfPackLU< M > solver(D); const Eigen::UmfPackLU<Sp> solver(D);
#else #else
const Eigen::SparseLU< M > solver(D); const Eigen::SparseLU<Sp> solver(D);
#endif #endif
// Compute value of eliminated variable. // Compute value of eliminated variable.
@ -190,6 +198,7 @@ namespace Opm
Eigen::SparseMatrix<double, Eigen::RowMajor>& A, Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
V& b) V& b)
{ {
/*
if (num_phases != 3) { if (num_phases != 3) {
OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases."); OPM_THROW(std::logic_error, "formEllipticSystem() requires 3 phases.");
} }
@ -269,6 +278,7 @@ namespace Opm
// Create output as product of L with equations. // Create output as product of L with equations.
A = L * total_residual.derivative()[0]; A = L * total_residual.derivative()[0];
b = L * total_residual.value().matrix(); b = L * total_residual.value().matrix();
*/
} }

6
opm/autodiff/VFPInjProperties.cpp
View File

@ -144,8 +144,10 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
} }
//Create diagonal matrices from ADB::Vs //Create diagonal matrices from ADB::Vs
ADB::M dthp_diag = spdiag(dthp); // ADB::M dthp_diag = spdiag(dthp);
ADB::M dflo_diag = spdiag(dflo); // ADB::M dflo_diag = spdiag(dflo);
ADB::M dthp_diag(dthp.matrix().asDiagonal());
ADB::M dflo_diag(dflo.matrix().asDiagonal());
//Calculate the Jacobians //Calculate the Jacobians
const int num_blocks = block_pattern.size(); const int num_blocks = block_pattern.size();

15
opm/autodiff/VFPProdProperties.cpp
View File

@ -143,11 +143,16 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
} }
//Create diagonal matrices from ADB::Vs //Create diagonal matrices from ADB::Vs
ADB::M dthp_diag = spdiag(dthp); // ADB::M dthp_diag = spdiag(dthp);
ADB::M dwfr_diag = spdiag(dwfr); // ADB::M dwfr_diag = spdiag(dwfr);
ADB::M dgfr_diag = spdiag(dgfr); // ADB::M dgfr_diag = spdiag(dgfr);
ADB::M dalq_diag = spdiag(dalq); // ADB::M dalq_diag = spdiag(dalq);
ADB::M dflo_diag = spdiag(dflo); // 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 //Calculate the Jacobians
const int num_blocks = block_pattern.size(); const int num_blocks = block_pattern.size();