opm-simulators/opm/autodiff/MultisegmentWells_impl.hpp

98 lines
3.8 KiB
C++

/*
Copyright 2016 SINTEF ICT, Applied Mathematics.
Copyright 2016 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED
namespace Opm
{
template <class WellState>
void
MultisegmentWells::
updateWellState(const Vector& dwells,
const int np,
const double dpmaxrel,
WellState& well_state) const
{
if (!wells().empty())
{
const int nw = wells().size();
const int nseg_total = nseg_total_;
// Extract parts of dwells corresponding to each part.
int varstart = 0;
const Vector dsegqs = subset(dwells, Span(np * nseg_total, 1, varstart));
varstart += dsegqs.size();
const Vector dsegp = subset(dwells, Span(nseg_total, 1, varstart));
varstart += dsegp.size();
assert(varstart == dwells.size());
// segment phase rates update
// in dwells, the phase rates are ordered by phase.
// while in WellStateMultiSegment, the phase rates are ordered by segments
const DataBlock wsr = Eigen::Map<const DataBlock>(dsegqs.data(), np, nseg_total).transpose();
const Vector dwsr = Eigen::Map<const Vector>(wsr.data(), nseg_total * np);
const Vector wsr_old = Eigen::Map<const Vector>(&well_state.segPhaseRates()[0], nseg_total * np);
const Vector sr = wsr_old - dwsr;
std::copy(&sr[0], &sr[0] + sr.size(), well_state.segPhaseRates().begin());
// segment pressure updates
const Vector segp_old = Eigen::Map<const Vector>(&well_state.segPress()[0], nseg_total, 1);
// TODO: applying the pressure change limiter to all the segments, not sure if it is the correct thing to do
const Vector dsegp_limited = sign(dsegp) * dsegp.abs().min(segp_old.abs() * dpmaxrel);
const Vector segp = segp_old - dsegp_limited;
std::copy(&segp[0], &segp[0] + segp.size(), well_state.segPress().begin());
// update the well rates and bhps, which are not anymore primary vabriables.
// they are updated directly from the updated segment phase rates and segment pressures.
// Bhp update.
Vector bhp = Vector::Zero(nw);
Vector wr = Vector::Zero(nw * np);
// it is better to use subset
int start_segment = 0;
for (int w = 0; w < nw; ++w) {
bhp[w] = well_state.segPress()[start_segment];
// insert can be faster
for (int p = 0; p < np; ++p) {
wr[p + np * w] = well_state.segPhaseRates()[p + np * start_segment];
}
const int nseg = wells()[w]->numberOfSegments();
start_segment += nseg;
}
assert(start_segment == nseg_total);
std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
// TODO: handling the THP control related.
}
}
}
#endif // OPM_MULTISEGMENTWELLS_IMPL_HEADER_INCLUDED