Merge pull request #334 from atgeirr/use-move-semantics

Use move semantics to gain performance
This commit is contained in:
dr-robertk 2015-03-20 10:47:07 +01:00
commit f288719e07
10 changed files with 131 additions and 117 deletions

View File

@ -90,8 +90,8 @@ phaseMobility(const Opm::IncompPropertiesInterface& props,
std::vector<M> dmw = { krwjac/mu[0] }; std::vector<M> dmw = { krwjac/mu[0] };
std::vector<M> dmo = { krojac/mu[1] }; std::vector<M> dmo = { krojac/mu[1] };
std::vector<ADB> pmobc = { ADB::function(krw / mu[0], dmw) , std::vector<ADB> pmobc = { ADB::function(krw / mu[0], std::move(dmw)) ,
ADB::function(kro / mu[1], dmo) }; ADB::function(kro / mu[1], std::move(dmo)) };
return pmobc; return pmobc;
} }

View File

@ -28,6 +28,7 @@
#include <opm/core/utility/platform_dependent/reenable_warnings.h> #include <opm/core/utility/platform_dependent/reenable_warnings.h>
#include <utility>
#include <vector> #include <vector>
#include <cassert> #include <cassert>
#include <iostream> #include <iostream>
@ -98,9 +99,14 @@ namespace Opm
/// Construct an empty AutoDiffBlock. /// Construct an empty AutoDiffBlock.
static AutoDiffBlock null() static AutoDiffBlock null()
{ {
V val; return AutoDiffBlock(V(), {});
std::vector<M> jac; }
return AutoDiffBlock(val, jac);
/// Create an AutoDiffBlock representing a constant.
/// \param[in] val values
static AutoDiffBlock constant(V&& val)
{
return AutoDiffBlock(std::move(val));
} }
/// Create an AutoDiffBlock representing a constant. /// Create an AutoDiffBlock representing a constant.
@ -126,7 +132,8 @@ namespace Opm
for (int i = 0; i < num_blocks; ++i) { for (int i = 0; i < num_blocks; ++i) {
jac[i] = M(num_elem, blocksizes[i]); jac[i] = M(num_elem, blocksizes[i]);
} }
return AutoDiffBlock(val, jac); V val_copy(val);
return AutoDiffBlock(std::move(val_copy), std::move(jac));
} }
/// Create an AutoDiffBlock representing a single variable block. /// Create an AutoDiffBlock representing a single variable block.
@ -136,7 +143,7 @@ namespace Opm
/// The resulting object will have size() equal to block_pattern[index]. /// The resulting object will have size() equal to block_pattern[index].
/// Its jacobians will all be zero, except for derivative()[index], which /// Its jacobians will all be zero, except for derivative()[index], which
/// will be an identity matrix. /// will be an identity matrix.
static AutoDiffBlock variable(const int index, const V& val, const std::vector<int>& blocksizes) static AutoDiffBlock variable(const int index, V&& val, const std::vector<int>& blocksizes)
{ {
std::vector<M> jac; std::vector<M> jac;
const int num_elem = val.size(); const int num_elem = val.size();
@ -164,15 +171,42 @@ namespace Opm
assert(jac[i].size()==0); assert(jac[i].size()==0);
} }
} }
return AutoDiffBlock(val, jac); return AutoDiffBlock(std::move(val), std::move(jac));
}
/// Create an AutoDiffBlock representing a single variable block.
/// \param[in] index index of the variable you are constructing
/// \param[in] val values
/// \param[in] blocksizes block pattern
/// The resulting object will have size() equal to block_pattern[index].
/// Its jacobians will all be zero, except for derivative()[index], which
/// will be an identity matrix.
static AutoDiffBlock variable(const int index, const V& val, const std::vector<int>& blocksizes)
{
V value = val;
return variable(index, std::move(value), blocksizes);
} }
/// Create an AutoDiffBlock by directly specifying values and jacobians. /// Create an AutoDiffBlock by directly specifying values and jacobians.
/// This version of function() moves its arguments and is therefore
/// quite efficient, but leaves the argument variables empty (but valid).
/// \param[in] val values
/// \param[in] jac vector of jacobians
static AutoDiffBlock function(V&& val, std::vector<M>&& jac)
{
return AutoDiffBlock(std::move(val), std::move(jac));
}
/// Create an AutoDiffBlock by directly specifying values and jacobians.
/// This version of function() copies its arguments and is therefore
/// less efficient than the other (moving) overload.
/// \param[in] val values /// \param[in] val values
/// \param[in] jac vector of jacobians /// \param[in] jac vector of jacobians
static AutoDiffBlock function(const V& val, const std::vector<M>& jac) static AutoDiffBlock function(const V& val, const std::vector<M>& jac)
{ {
return AutoDiffBlock(val, jac); V val_copy(val);
std::vector<M> jac_copy(jac);
return AutoDiffBlock(std::move(val_copy), std::move(jac_copy));
} }
/// Construct a set of primary variables, each initialized to /// Construct a set of primary variables, each initialized to
@ -259,7 +293,7 @@ namespace Opm
assert(jac[block].cols() == rhs.jac_[block].cols()); assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] += rhs.jac_[block]; jac[block] += rhs.jac_[block];
} }
return function(val_ + rhs.val_, jac); return function(val_ + rhs.val_, std::move(jac));
} }
/// Elementwise operator - /// Elementwise operator -
@ -282,7 +316,7 @@ namespace Opm
assert(jac[block].cols() == rhs.jac_[block].cols()); assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] -= rhs.jac_[block]; jac[block] -= rhs.jac_[block];
} }
return function(val_ - rhs.val_, jac); return function(val_ - rhs.val_, std::move(jac));
} }
/// Elementwise operator * /// Elementwise operator *
@ -318,7 +352,7 @@ namespace Opm
jac[block] = D2*jac_[block] + D1*rhs.jac_[block]; jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
} }
} }
return function(val_ * rhs.val_, jac); return function(val_ * rhs.val_, std::move(jac));
} }
/// Elementwise operator / /// Elementwise operator /
@ -358,7 +392,7 @@ namespace Opm
jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]); jac[block] = D3 * (D2*jac_[block] - D1*rhs.jac_[block]);
} }
} }
return function(val_ / rhs.val_, jac); return function(val_ / rhs.val_, std::move(jac));
} }
/// I/O. /// I/O.
@ -374,6 +408,13 @@ namespace Opm
return os; return os;
} }
/// Efficient swap function.
void swap(AutoDiffBlock& other)
{
val_.swap(other.val_);
jac_.swap(other.jac_);
}
/// Number of elements /// Number of elements
int size() const int size() const
{ {
@ -411,13 +452,17 @@ namespace Opm
private: private:
AutoDiffBlock(const V& val) AutoDiffBlock(const V& val)
: val_(val) : val_(val)
{ {
} }
AutoDiffBlock(const V& val, AutoDiffBlock(V&& val)
const std::vector<M>& jac) : val_(std::move(val))
: val_(val), jac_(jac) {
}
AutoDiffBlock(V&& val, std::vector<M>&& jac)
: val_(std::move(val)), jac_(std::move(jac))
{ {
#ifndef NDEBUG #ifndef NDEBUG
const int num_elem = val_.size(); const int num_elem = val_.size();
@ -457,7 +502,7 @@ namespace Opm
fastSparseProduct(lhs, rhs.derivative()[block], jac[block]); fastSparseProduct(lhs, rhs.derivative()[block], jac[block]);
} }
typename AutoDiffBlock<Scalar>::V val = lhs*rhs.value().matrix(); typename AutoDiffBlock<Scalar>::V val = lhs*rhs.value().matrix();
return AutoDiffBlock<Scalar>::function(val, jac); return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));
} }
@ -549,7 +594,7 @@ namespace Opm
for (int block=0; block<lhs.numBlocks(); block++) { for (int block=0; block<lhs.numBlocks(); block++) {
jac.emplace_back( lhs.derivative()[block] * rhs ); jac.emplace_back( lhs.derivative()[block] * rhs );
} }
return AutoDiffBlock<Scalar>::function( lhs.value() * rhs, jac ); return AutoDiffBlock<Scalar>::function( lhs.value() * rhs, std::move(jac) );
} }

View File

@ -403,7 +403,8 @@ collapseJacs(const AutoDiffBlock<double>& x)
// Build final jacobian. // Build final jacobian.
std::vector<ADB::M> jacs(1); std::vector<ADB::M> jacs(1);
collapseJacs( x, jacs[ 0 ] ); collapseJacs( x, jacs[ 0 ] );
return ADB::function(x.value(), jacs); ADB::V val = x.value();
return ADB::function(std::move(val), std::move(jacs));
} }

View File

@ -335,7 +335,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
props_[phase_usage_.phase_pos[Water]]->mu(n, pvt_region_.data(), pw.value().data(), T.value().data(), rs, props_[phase_usage_.phase_pos[Water]]->mu(n, pvt_region_.data(), pw.value().data(), T.value().data(), rs,
mu.data(), dmudp.data(), dmudr.data()); mu.data(), dmudp.data(), dmudr.data());
if (pw.derivative().empty()) { if (pw.derivative().empty()) {
return ADB::constant(mu); return ADB::constant(std::move(mu));
} else { } else {
ADB::M dmudp_diag = spdiag(dmudp); ADB::M dmudp_diag = spdiag(dmudp);
const int num_blocks = pw.numBlocks(); const int num_blocks = pw.numBlocks();
@ -343,7 +343,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dmudp_diag, pw.derivative()[block], jacs[block]); fastSparseProduct(dmudp_diag, pw.derivative()[block], jacs[block]);
} }
return ADB::function(mu, jacs); return ADB::function(std::move(mu), std::move(jacs));
} }
} }
@ -383,7 +383,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dmudr_diag, rs.derivative()[block], temp); fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
jacs[block] += temp; jacs[block] += temp;
} }
return ADB::function(mu, jacs); return ADB::function(std::move(mu), std::move(jacs));
} }
/// Gas viscosity. /// Gas viscosity.
@ -422,7 +422,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dmudr_diag, rv.derivative()[block], temp); fastSparseProduct(dmudr_diag, rv.derivative()[block], temp);
jacs[block] += temp; jacs[block] += temp;
} }
return ADB::function(mu, jacs); return ADB::function(std::move(mu), std::move(jacs));
} }
@ -459,7 +459,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dbdp_diag, pw.derivative()[block], jacs[block]); fastSparseProduct(dbdp_diag, pw.derivative()[block], jacs[block]);
} }
return ADB::function(b, jacs); return ADB::function(std::move(b), std::move(jacs));
} }
/// Oil formation volume factor. /// Oil formation volume factor.
@ -499,7 +499,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dbdr_diag, rs.derivative()[block], temp); fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
jacs[block] += temp; jacs[block] += temp;
} }
return ADB::function(b, jacs); return ADB::function(std::move(b), std::move(jacs));
} }
/// Gas formation volume factor. /// Gas formation volume factor.
@ -539,7 +539,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dbdr_diag, rv.derivative()[block], temp); fastSparseProduct(dbdr_diag, rv.derivative()[block], temp);
jacs[block] += temp; jacs[block] += temp;
} }
return ADB::function(b, jacs); return ADB::function(std::move(b), std::move(jacs));
} }
@ -568,7 +568,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(drbubdp_diag, po.derivative()[block], jacs[block]); fastSparseProduct(drbubdp_diag, po.derivative()[block], jacs[block]);
} }
return ADB::function(rbub, jacs); return ADB::function(std::move(rbub), std::move(jacs));
} }
/// Bubble point curve for Rs as function of oil pressure. /// Bubble point curve for Rs as function of oil pressure.
@ -609,7 +609,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(drvdp_diag, po.derivative()[block], jacs[block]); fastSparseProduct(drvdp_diag, po.derivative()[block], jacs[block]);
} }
return ADB::function(rv, jacs); return ADB::function(std::move(rv), std::move(jacs));
} }
/// Condensation curve for Rv as function of oil pressure. /// Condensation curve for Rv as function of oil pressure.
@ -686,7 +686,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
jacs[block] += temp; jacs[block] += temp;
} }
} }
relperms.emplace_back(ADB::function(kr.col(phase1_pos), jacs)); ADB::V val = kr.col(phase1_pos);
relperms.emplace_back(ADB::function(std::move(val), std::move(jacs)));
} else { } else {
relperms.emplace_back(ADB::null()); relperms.emplace_back(ADB::null());
} }
@ -746,7 +747,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
jacs[block] += temp; jacs[block] += temp;
} }
} }
adbCapPressures.emplace_back(ADB::function(pc.col(phase1_pos), jacs)); ADB::V val = pc.col(phase1_pos);
adbCapPressures.emplace_back(ADB::function(std::move(val), std::move(jacs)));
} else { } else {
adbCapPressures.emplace_back(ADB::null()); adbCapPressures.emplace_back(ADB::null());
} }
@ -849,7 +851,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
jacs[block] = dfactor_dso_diag * so.derivative()[block]; jacs[block] = dfactor_dso_diag * so.derivative()[block];
} }
r = ADB::function(factor, jacs)*r; r = ADB::function(std::move(factor), std::move(jacs))*r;
} }
} }

View File

@ -281,11 +281,6 @@ namespace Opm {
std::vector<ADB> std::vector<ADB>
computeRelPerm(const SolutionState& state) const; computeRelPerm(const SolutionState& state) const;
std::vector<ADB>
computeRelPermWells(const SolutionState& state,
const DataBlock& well_s,
const std::vector<int>& well_cells) const;
void void
computeMassFlux(const int actph , computeMassFlux(const int actph ,
const V& transi, const V& transi,

View File

@ -542,21 +542,20 @@ namespace detail {
// Pressure. // Pressure.
int nextvar = 0; int nextvar = 0;
state.pressure = vars[ nextvar++ ]; state.pressure = std::move(vars[ nextvar++ ]);
// temperature // temperature
const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size()); const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
state.temperature = ADB::constant(temp); state.temperature = ADB::constant(temp);
// Saturations // Saturations
const std::vector<int>& bpat = vars[0].blockPattern();
{ {
ADB so = ADB::constant(V::Ones(nc, 1), bpat); ADB so = ADB::constant(V::Ones(nc, 1));
if (active_[ Water ]) { if (active_[ Water ]) {
ADB& sw = vars[ nextvar++ ]; state.saturation[pu.phase_pos[ Water ]] = std::move(vars[ nextvar++ ]);
state.saturation[pu.phase_pos[ Water ]] = sw; const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
so -= sw; so -= sw;
} }
@ -572,7 +571,7 @@ namespace detail {
// RS and RV is only defined if both oil and gas phase are active. // RS and RV is only defined if both oil and gas phase are active.
const ADB& sw = (active_[ Water ] const ADB& sw = (active_[ Water ]
? state.saturation[ pu.phase_pos[ Water ] ] ? state.saturation[ pu.phase_pos[ Water ] ]
: ADB::constant(V::Zero(nc, 1), bpat)); : ADB::constant(V::Zero(nc, 1)));
state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg); state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_); const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
if (has_disgas_) { if (has_disgas_) {
@ -591,15 +590,15 @@ namespace detail {
if (active_[ Oil ]) { if (active_[ Oil ]) {
// Note that so is never a primary variable. // Note that so is never a primary variable.
state.saturation[pu.phase_pos[ Oil ]] = so; state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
} }
} }
// Qs. // Qs.
state.qs = vars[ nextvar++ ]; state.qs = std::move(vars[ nextvar++ ]);
// Bhp. // Bhp.
state.bhp = vars[ nextvar++ ]; state.bhp = std::move(vars[ nextvar++ ]);
assert(nextvar == int(vars.size())); assert(nextvar == int(vars.size()));
@ -946,7 +945,7 @@ namespace detail {
} }
// total rates at std // total rates at std
ADB qt_s = ADB::constant(V::Zero(nw), state.bhp.blockPattern()); ADB qt_s = ADB::constant(V::Zero(nw));
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
qt_s += subset(state.qs, Span(nw, 1, phase*nw)); qt_s += subset(state.qs, Span(nw, 1, phase*nw));
} }
@ -954,7 +953,7 @@ namespace detail {
// compute avg. and total wellbore phase volumetric rates at std. conds // compute avg. and total wellbore phase volumetric rates at std. conds
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np); const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
std::vector<ADB> wbq(np, ADB::null()); std::vector<ADB> wbq(np, ADB::null());
ADB wbqt = ADB::constant(V::Zero(nw), state.pressure.blockPattern()); ADB wbqt = ADB::constant(V::Zero(nw));
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
const int pos = pu.phase_pos[phase]; const int pos = pu.phase_pos[phase];
wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase]; wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase];
@ -995,7 +994,7 @@ namespace detail {
ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown); ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown);
// compute volume ratio between connection at standard conditions // compute volume ratio between connection at standard conditions
ADB volRat = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); ADB volRat = ADB::constant(V::Zero(nperf));
std::vector<ADB> cmix_s(np, ADB::null()); std::vector<ADB> cmix_s(np, ADB::null());
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
cmix_s[phase] = wops_.w2p * mix_s[phase]; cmix_s[phase] = wops_.w2p * mix_s[phase];
@ -1223,16 +1222,22 @@ namespace detail {
// Update primary variables, if necessary. // Update primary variables, if necessary.
if (bhp_changed) { if (bhp_changed) {
// We will set the bhp primary variable to the new ones,
// but we do not change the derivatives here.
ADB::V new_bhp = Eigen::Map<ADB::V>(xw.bhp().data(), nw); ADB::V new_bhp = Eigen::Map<ADB::V>(xw.bhp().data(), nw);
bhp = ADB::function(new_bhp, bhp.derivative()); // Avoiding the copy below would require a value setter method
// in AutoDiffBlock.
std::vector<ADB::M> old_derivs = bhp.derivative();
bhp = ADB::function(std::move(new_bhp), std::move(old_derivs));
} }
if (rates_changed) { if (rates_changed) {
// Need to reshuffle well rates, from phase running fastest // Need to reshuffle well rates, from phase running fastest
// to wells running fastest. // to wells running fastest.
// The transpose() below switches the ordering. // The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(xw.wellRates().data(), nw, np).transpose(); const DataBlock wrates = Eigen::Map<const DataBlock>(xw.wellRates().data(), nw, np).transpose();
const ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np); ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np);
well_phase_flow_rate = ADB::function(new_qs, well_phase_flow_rate.derivative()); std::vector<ADB::M> old_derivs = well_phase_flow_rate.derivative();
well_phase_flow_rate = ADB::function(std::move(new_qs), std::move(old_derivs));
} }
} }
@ -1607,22 +1612,21 @@ namespace detail {
{ {
using namespace Opm::AutoDiffGrid; using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_); const int nc = numCells(grid_);
const std::vector<int>& bpat = state.pressure.blockPattern();
const ADB null = ADB::constant(V::Zero(nc, 1), bpat); const ADB zero = ADB::constant(V::Zero(nc));
const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB sw = (active_[ Water ] const ADB& sw = (active_[ Water ]
? state.saturation[ pu.phase_pos[ Water ] ] ? state.saturation[ pu.phase_pos[ Water ] ]
: null); : zero);
const ADB so = (active_[ Oil ] const ADB& so = (active_[ Oil ]
? state.saturation[ pu.phase_pos[ Oil ] ] ? state.saturation[ pu.phase_pos[ Oil ] ]
: null); : zero);
const ADB sg = (active_[ Gas ] const ADB& sg = (active_[ Gas ]
? state.saturation[ pu.phase_pos[ Gas ] ] ? state.saturation[ pu.phase_pos[ Gas ] ]
: null); : zero);
return fluid_.relperm(sw, so, sg, cells_); return fluid_.relperm(sw, so, sg, cells_);
} }
@ -1637,9 +1641,8 @@ namespace detail {
{ {
using namespace Opm::AutoDiffGrid; using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_); const int nc = numCells(grid_);
const std::vector<int>& bpat = state.pressure.blockPattern();
const ADB null = ADB::constant(V::Zero(nc, 1), bpat); const ADB null = ADB::constant(V::Zero(nc));
const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB& sw = (active_[ Water ] const ADB& sw = (active_[ Water ]
@ -1718,38 +1721,6 @@ namespace detail {
template<class T>
std::vector<ADB>
FullyImplicitBlackoilSolver<T>::computeRelPermWells(const SolutionState& state,
const DataBlock& well_s,
const std::vector<int>& well_cells) const
{
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const std::vector<int>& bpat = state.pressure.blockPattern();
const ADB null = ADB::constant(V::Zero(nperf), bpat);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB sw = (active_[ Water ]
? ADB::constant(well_s.col(pu.phase_pos[ Water ]), bpat)
: null);
const ADB so = (active_[ Oil ]
? ADB::constant(well_s.col(pu.phase_pos[ Oil ]), bpat)
: null);
const ADB sg = (active_[ Gas ]
? ADB::constant(well_s.col(pu.phase_pos[ Gas ]), bpat)
: null);
return fluid_.relperm(sw, so, sg, well_cells);
}
template<class T> template<class T>
void void
FullyImplicitBlackoilSolver<T>::computeMassFlux(const int actph , FullyImplicitBlackoilSolver<T>::computeMassFlux(const int actph ,
@ -2263,9 +2234,9 @@ namespace detail {
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]); fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]);
} }
return ADB::function(pm, jacs); return ADB::function(std::move(pm), std::move(jacs));
} else { } else {
return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); return ADB::constant(V::Constant(n, 1.0));
} }
} }
@ -2291,9 +2262,9 @@ namespace detail {
for (int block = 0; block < num_blocks; ++block) { for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]); fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]);
} }
return ADB::function(tm, jacs); return ADB::function(std::move(tm), std::move(jacs));
} else { } else {
return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); return ADB::constant(V::Constant(n, 1.0));
} }
} }

View File

@ -326,7 +326,6 @@ namespace {
const ADB& p = vars[0]; const ADB& p = vars[0];
const ADB T = ADB::constant(T0); const ADB T = ADB::constant(T0);
const ADB& bhp = vars[1]; const ADB& bhp = vars[1];
std::vector<int> bpat = p.blockPattern();
// Compute T_ij * (p_i - p_j). // Compute T_ij * (p_i - p_j).
const ADB nkgradp = transi * (ops_.ngrad * p); const ADB nkgradp = transi * (ops_.ngrad * p);
@ -351,10 +350,10 @@ namespace {
const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); const ADB nkgradp_well = transw * (p_perfcell - p_perfwell);
const Selector<double> cell_to_well_selector(nkgradp_well.value()); const Selector<double> cell_to_well_selector(nkgradp_well.value());
cell_residual_ = ADB::constant(pv, bpat); cell_residual_ = ADB::constant(pv);
well_residual_ = ADB::constant(V::Zero(nw,1), bpat); well_residual_ = ADB::constant(V::Zero(nw,1));
ADB divcontrib_sum = ADB::constant(V::Zero(nc,1), bpat); ADB divcontrib_sum = ADB::constant(V::Zero(nc,1));
qs_ = ADB::constant(V::Zero(nw*np, 1), bpat); qs_ = ADB::constant(V::Zero(nw*np, 1));
for (int phase = 0; phase < np; ++phase) { for (int phase = 0; phase < np; ++phase) {
const ADB cell_b = fluidFvf(phase, p, T, cells); const ADB cell_b = fluidFvf(phase, p, T, cells);
const ADB cell_rho = fluidRho(phase, p, T, cells); const ADB cell_rho = fluidRho(phase, p, T, cells);
@ -381,7 +380,7 @@ namespace {
const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc); const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
const ADB divcontrib = delta_t * (ops_.div * (flux * face_b) + well_contrib); const ADB divcontrib = delta_t * (ops_.div * (flux * face_b) + well_contrib);
const V qcontrib = delta_t * q; const V qcontrib = delta_t * q;
const ADB pvcontrib = ADB::constant(pv*z0, bpat); const ADB pvcontrib = ADB::constant(pv*z0);
const ADB component_contrib = pvcontrib + qcontrib; const ADB component_contrib = pvcontrib + qcontrib;
divcontrib_sum = divcontrib_sum - divcontrib/cell_b; divcontrib_sum = divcontrib_sum - divcontrib/cell_b;
cell_residual_ = cell_residual_ - (component_contrib/cell_b); cell_residual_ = cell_residual_ - (component_contrib/cell_b);

View File

@ -314,7 +314,7 @@ namespace Opm
if (eq == n) { if (eq == n) {
continue; continue;
} }
retval.push_back(ADB::function(vals[eq], jacs[eq])); retval.push_back(ADB::function(std::move(vals[eq]), std::move(jacs[eq])));
} }
return retval; return retval;
} }
@ -338,7 +338,8 @@ namespace Opm
// 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);
ADB eq_coll = collapseJacs(ADB::function(equation.value(), C_jacs)); V equation_value = equation.value();
ADB eq_coll = collapseJacs(ADB::function(std::move(equation_value), std::move(C_jacs)));
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
@ -395,7 +396,7 @@ namespace Opm
// A concession to MRST, to obtain more similar behaviour: // A concession to MRST, to obtain more similar behaviour:
// swap the first two equations, so that oil is first, then water. // swap the first two equations, so that oil is first, then water.
auto eqs = eqs_in; auto eqs = eqs_in;
std::swap(eqs[0], eqs[1]); eqs[0].swap(eqs[1]);
// Characterize the material balance equations. // Characterize the material balance equations.
const int n = eqs[0].size(); const int n = eqs[0].size();
@ -462,7 +463,7 @@ namespace Opm
L.setFromTriplets(t.begin(), t.end()); L.setFromTriplets(t.begin(), t.end());
// Combine in single block. // Combine in single block.
ADB total_residual = eqs[0]; ADB total_residual = std::move(eqs[0]);
for (int phase = 1; phase < num_phases; ++phase) { for (int phase = 1; phase < num_phases; ++phase) {
total_residual = vertcat(total_residual, eqs[phase]); total_residual = vertcat(total_residual, eqs[phase]);
} }

View File

@ -129,8 +129,8 @@ namespace Opm
std::vector<M> dmw = { krwjac/mu[0] }; std::vector<M> dmw = { krwjac/mu[0] };
std::vector<M> dmo = { krojac/mu[1] }; std::vector<M> dmo = { krojac/mu[1] };
std::vector<ADB> pmobc = { ADB::function(krw / mu[0], dmw) , std::vector<ADB> pmobc = { ADB::function(krw / mu[0], std::move(dmw)) ,
ADB::function(kro / mu[1], dmo) }; ADB::function(kro / mu[1], std::move(dmo)) };
return pmobc; return pmobc;
} }

View File

@ -77,12 +77,10 @@ BOOST_AUTO_TEST_CASE(ConstantInitialisation)
{ {
typedef AutoDiffBlock<double> ADB; typedef AutoDiffBlock<double> ADB;
std::vector<int> blocksizes = { 3, 1, 2 };
ADB::V v(3); ADB::V v(3);
v << 0.2, 1.2, 13.4; v << 0.2, 1.2, 13.4;
ADB a = ADB::constant(v, blocksizes); ADB a = ADB::constant(v);
BOOST_REQUIRE(a.value().matrix() == v.matrix()); BOOST_REQUIRE(a.value().matrix() == v.matrix());
const std::vector<ADB::M>& J = a.derivative(); const std::vector<ADB::M>& J = a.derivative();
@ -137,7 +135,9 @@ BOOST_AUTO_TEST_CASE(FunctionInitialisation)
jacs[j].insert(0,0) = -1.0; jacs[j].insert(0,0) = -1.0;
} }
ADB f = ADB::function(v, jacs); ADB::V v_copy(v);
std::vector<ADB::M> jacs_copy(jacs);
ADB f = ADB::function(std::move(v_copy), std::move(jacs_copy));
BOOST_REQUIRE(f.value().matrix() == v.matrix()); BOOST_REQUIRE(f.value().matrix() == v.matrix());