Make more use of move semantics in AD code.

This makes some API changes to AutoDiffBlock.
 - Add overload for the constant() constructor taking rvalue ref.
 - Add overload for the variable() constructor taking rvalue ref.
 - Make the function() constructor *require* rvalue refs.
 - Add a swap() function.

The remaining changes in this commit are follow-ups especially
to the third change (adding std::move in many places), and
some removal of unnecessary block pattern arguments from calls to
the constant() static method.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-03-10 12:36:14 +01:00
parent 635f3db814
commit 04b255a03f
9 changed files with 100 additions and 64 deletions

View File

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

View File

@ -98,9 +98,14 @@ namespace Opm
/// Construct an empty AutoDiffBlock.
static AutoDiffBlock null()
{
V val;
std::vector<M> jac;
return AutoDiffBlock(val, jac);
return AutoDiffBlock(V(), {});
}
/// 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.
@ -126,7 +131,8 @@ namespace Opm
for (int i = 0; i < num_blocks; ++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.
@ -136,7 +142,7 @@ namespace Opm
/// 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)
static AutoDiffBlock variable(const int index, V&& val, const std::vector<int>& blocksizes)
{
std::vector<M> jac;
const int num_elem = val.size();
@ -164,15 +170,28 @@ namespace Opm
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.
/// \param[in] val values
/// \param[in] jac vector of jacobians
static AutoDiffBlock function(const V& val, const std::vector<M>& jac)
static AutoDiffBlock function(V&& val, std::vector<M>&& jac)
{
return AutoDiffBlock(val, jac);
return AutoDiffBlock(std::move(val), std::move(jac));
}
/// Construct a set of primary variables, each initialized to
@ -259,7 +278,7 @@ namespace Opm
assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] += rhs.jac_[block];
}
return function(val_ + rhs.val_, jac);
return function(val_ + rhs.val_, std::move(jac));
}
/// Elementwise operator -
@ -282,7 +301,7 @@ namespace Opm
assert(jac[block].cols() == rhs.jac_[block].cols());
jac[block] -= rhs.jac_[block];
}
return function(val_ - rhs.val_, jac);
return function(val_ - rhs.val_, std::move(jac));
}
/// Elementwise operator *
@ -318,7 +337,7 @@ namespace Opm
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 /
@ -358,7 +377,7 @@ namespace Opm
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.
@ -374,6 +393,13 @@ namespace Opm
return os;
}
/// Efficient swap function.
void swap(AutoDiffBlock& other)
{
val_.swap(other.val_);
jac_.swap(other.jac_);
}
/// Number of elements
int size() const
{
@ -411,13 +437,17 @@ namespace Opm
private:
AutoDiffBlock(const V& val)
: val_(val)
: val_(val)
{
}
AutoDiffBlock(const V& val,
const std::vector<M>& jac)
: val_(val), jac_(jac)
AutoDiffBlock(V&& val)
: val_(std::move(val))
{
}
AutoDiffBlock(V&& val, std::vector<M>&& jac)
: val_(std::move(val)), jac_(std::move(jac))
{
#ifndef NDEBUG
const int num_elem = val_.size();
@ -457,7 +487,7 @@ namespace Opm
fastSparseProduct(lhs, rhs.derivative()[block], jac[block]);
}
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 +579,7 @@ namespace Opm
for (int block=0; block<lhs.numBlocks(); block++) {
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.
std::vector<ADB::M> jacs(1);
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,
mu.data(), dmudp.data(), dmudr.data());
if (pw.derivative().empty()) {
return ADB::constant(mu);
return ADB::constant(std::move(mu));
} else {
ADB::M dmudp_diag = spdiag(dmudp);
const int num_blocks = pw.numBlocks();
@ -343,7 +343,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++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);
jacs[block] += temp;
}
return ADB::function(mu, jacs);
return ADB::function(std::move(mu), std::move(jacs));
}
/// Gas viscosity.
@ -422,7 +422,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dmudr_diag, rv.derivative()[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) {
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.
@ -499,7 +499,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
jacs[block] += temp;
}
return ADB::function(b, jacs);
return ADB::function(std::move(b), std::move(jacs));
}
/// Gas formation volume factor.
@ -539,7 +539,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dbdr_diag, rv.derivative()[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) {
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.
@ -609,7 +609,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++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.
@ -686,7 +686,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
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 {
relperms.emplace_back(ADB::null());
}
@ -746,7 +747,8 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
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 {
adbCapPressures.emplace_back(ADB::null());
}
@ -849,7 +851,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++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

@ -945,7 +945,7 @@ namespace detail {
}
// 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) {
qt_s += subset(state.qs, Span(nw, 1, phase*nw));
}
@ -953,7 +953,7 @@ namespace detail {
// compute avg. and total wellbore phase volumetric rates at std. conds
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
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) {
const int pos = pu.phase_pos[phase];
wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase];
@ -994,7 +994,7 @@ namespace detail {
ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown);
// 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());
for (int phase = 0; phase < np; ++phase) {
cmix_s[phase] = wops_.w2p * mix_s[phase];
@ -1222,16 +1222,22 @@ namespace detail {
// Update primary variables, if necessary.
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);
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) {
// Need to reshuffle well rates, from phase running fastest
// to wells running fastest.
// The transpose() below switches the ordering.
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);
well_phase_flow_rate = ADB::function(new_qs, well_phase_flow_rate.derivative());
ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np);
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));
}
}
@ -1635,9 +1641,8 @@ namespace detail {
{
using namespace Opm::AutoDiffGrid;
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 ADB& sw = (active_[ Water ]
@ -1724,21 +1729,20 @@ namespace detail {
{
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 ADB null = ADB::constant(V::Zero(nperf));
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const ADB sw = (active_[ Water ]
? ADB::constant(well_s.col(pu.phase_pos[ Water ]), bpat)
? ADB::constant(well_s.col(pu.phase_pos[ Water ]))
: null);
const ADB so = (active_[ Oil ]
? ADB::constant(well_s.col(pu.phase_pos[ Oil ]), bpat)
? ADB::constant(well_s.col(pu.phase_pos[ Oil ]))
: null);
const ADB sg = (active_[ Gas ]
? ADB::constant(well_s.col(pu.phase_pos[ Gas ]), bpat)
? ADB::constant(well_s.col(pu.phase_pos[ Gas ]))
: null);
return fluid_.relperm(sw, so, sg, well_cells);
@ -2261,9 +2265,9 @@ namespace detail {
for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]);
}
return ADB::function(pm, jacs);
return ADB::function(std::move(pm), std::move(jacs));
} else {
return ADB::constant(V::Constant(n, 1.0), p.blockPattern());
return ADB::constant(V::Constant(n, 1.0));
}
}
@ -2289,9 +2293,9 @@ namespace detail {
for (int block = 0; block < num_blocks; ++block) {
fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]);
}
return ADB::function(tm, jacs);
return ADB::function(std::move(tm), std::move(jacs));
} 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 T = ADB::constant(T0);
const ADB& bhp = vars[1];
std::vector<int> bpat = p.blockPattern();
// Compute T_ij * (p_i - p_j).
const ADB nkgradp = transi * (ops_.ngrad * p);
@ -351,10 +350,10 @@ namespace {
const ADB nkgradp_well = transw * (p_perfcell - p_perfwell);
const Selector<double> cell_to_well_selector(nkgradp_well.value());
cell_residual_ = ADB::constant(pv, bpat);
well_residual_ = ADB::constant(V::Zero(nw,1), bpat);
ADB divcontrib_sum = ADB::constant(V::Zero(nc,1), bpat);
qs_ = ADB::constant(V::Zero(nw*np, 1), bpat);
cell_residual_ = ADB::constant(pv);
well_residual_ = ADB::constant(V::Zero(nw,1));
ADB divcontrib_sum = ADB::constant(V::Zero(nc,1));
qs_ = ADB::constant(V::Zero(nw*np, 1));
for (int phase = 0; phase < np; ++phase) {
const ADB cell_b = fluidFvf(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 divcontrib = delta_t * (ops_.div * (flux * face_b) + well_contrib);
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;
divcontrib_sum = divcontrib_sum - divcontrib/cell_b;
cell_residual_ = cell_residual_ - (component_contrib/cell_b);

View File

@ -314,7 +314,7 @@ namespace Opm
if (eq == n) {
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;
}
@ -338,7 +338,8 @@ namespace Opm
// Build C.
std::vector<M> C_jacs = equation.derivative();
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];
// Use sparse LU to solve the block submatrices

View File

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

View File

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