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

Fix issues for iterated sequential implicit.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2017-05-03 14:20:44 +02:00
parent e01eca28ca
commit e0570fc980
1 changed files with 21 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

64
opm/autodiff/BlackoilTransportModel.hpp
View File

@ -74,16 +74,25 @@ namespace Opm {
{ {
} }
void prepareStep(const SimulatorTimerInterface& timer, void prepareStep(const SimulatorTimerInterface& timer,
const ReservoirState& reservoir_state, const ReservoirState& reservoir_state,
const WellState& well_state) const WellState& well_state)
{ {
Base::prepareStep(timer, reservoir_state, well_state); Base::prepareStep(timer, reservoir_state, well_state);
Base::param_.solve_welleq_initially_ = false; Base::param_.solve_welleq_initially_ = false;
state0_ = variableState(reservoir_state, well_state); SolutionState state0 = variableState(reservoir_state, well_state);
asImpl().makeConstantState(state0_); asImpl().makeConstantState(state0);
asImpl().computeAccum(state0, 0);
} }
SimulatorReport SimulatorReport
assemble(const ReservoirState& reservoir_state, assemble(const ReservoirState& reservoir_state,
WellState& well_state, WellState& well_state,
@ -107,11 +116,7 @@ namespace Opm {
// get reasonable initial conditions for the wells // get reasonable initial conditions for the wells
asImpl().wellModel().updateWellControls(well_state); asImpl().wellModel().updateWellControls(well_state);
// Create the primary variables.
SolutionState state = asImpl().variableState(reservoir_state, well_state);
if (initial_assembly) { if (initial_assembly) {
// HACK // HACK
const_cast<V&>(total_flux_) const_cast<V&>(total_flux_)
= Eigen::Map<const V>(reservoir_state.faceflux().data(), reservoir_state.faceflux().size()); = Eigen::Map<const V>(reservoir_state.faceflux().data(), reservoir_state.faceflux().size());
@ -120,19 +125,19 @@ namespace Opm {
const_cast<DataBlock&>(comp_wellperf_flux_) const_cast<DataBlock&>(comp_wellperf_flux_)
= Eigen::Map<const DataBlock>(well_state.perfPhaseRates().data(), well_state.perfRates().size(), numPhases()); = Eigen::Map<const DataBlock>(well_state.perfPhaseRates().data(), well_state.perfRates().size(), numPhases());
assert(numPhases() * well_state.perfRates().size() == well_state.perfPhaseRates().size()); assert(numPhases() * well_state.perfRates().size() == well_state.perfPhaseRates().size());
asImpl().updatePhaseCondFromPrimalVariable(reservoir_state);
}
is_first_iter_ = true; // Create the primary variables.
// Create the (constant, derivativeless) initial state. SolutionState state = asImpl().variableState(reservoir_state, well_state);
if (initial_assembly) {
SolutionState state0 = state; SolutionState state0 = state;
asImpl().makeConstantState(state0); asImpl().makeConstantState(state0);
// Compute initial accumulation contributions
// and well connection pressures.
asImpl().computeAccum(state0, 0);
asImpl().wellModel().computeWellConnectionPressures(state0, well_state); asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
} else {
is_first_iter_ = false;
} }
// -------- Mass balance equations -------- // -------- Mass balance equations --------
asImpl().assembleMassBalanceEq(state); asImpl().assembleMassBalanceEq(state);
@ -229,11 +234,12 @@ namespace Opm {
V total_flux_; // HACK, should be part of a revised (transport-specific) SolutionState. V total_flux_; // HACK, should be part of a revised (transport-specific) SolutionState.
V total_wellperf_flux_; V total_wellperf_flux_;
DataBlock comp_wellperf_flux_; DataBlock comp_wellperf_flux_;
SolutionState state0_ = SolutionState(3);
bool is_first_iter_ = false;
Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind_flags_; Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> upwind_flags_;
SolutionState SolutionState
variableState(const ReservoirState& x, variableState(const ReservoirState& x,
const WellState& xw) const const WellState& xw) const
@ -252,26 +258,6 @@ namespace Opm {
void computeAccum(const SolutionState& state,
const int aix )
{
if (aix == 0) {
// The pressure passed in state is from after
// the pressure solver, but we need to use the original
// b factors etc. to get the initial accumulation term
// correct.
Base::computeAccum(state0_, aix);
} else {
Base::computeAccum(state, aix);
}
}
void assembleMassBalanceEq(const SolutionState& state) void assembleMassBalanceEq(const SolutionState& state)
@ -293,10 +279,6 @@ namespace Opm {
const ADB tr_mult = asImpl().transMult(state.pressure); const ADB tr_mult = asImpl().transMult(state.pressure);
const V gdz = geo_.gravity()[2] * (ops_.grad * geo_.z().matrix()); const V gdz = geo_.gravity()[2] * (ops_.grad * geo_.z().matrix());
if (is_first_iter_) {
upwind_flags_.resize(gdz.size(), numPhases());
}
// Compute mobilities and heads // Compute mobilities and heads
const std::vector<PhasePresence>& cond = asImpl().phaseCondition(); const std::vector<PhasePresence>& cond = asImpl().phaseCondition();
const std::vector<ADB> kr = asImpl().computeRelPerm(state); const std::vector<ADB> kr = asImpl().computeRelPerm(state);
@ -316,10 +298,6 @@ namespace Opm {
const ADB rhoavg = ops_.caver * sd_.rq[phase_idx].rho; const ADB rhoavg = ops_.caver * sd_.rq[phase_idx].rho;
sd_.rq[ phase_idx ].dh = ops_.grad * phase_pressure - rhoavg * gdz; sd_.rq[ phase_idx ].dh = ops_.grad * phase_pressure - rhoavg * gdz;
if (is_first_iter_) {
upwind_flags_.col(phase_idx) = -sd_.rq[phase_idx].dh.value();
}
if (use_threshold_pressure_) { if (use_threshold_pressure_) {
asImpl().applyThresholdPressures(sd_.rq[ phase_idx ].dh); asImpl().applyThresholdPressures(sd_.rq[ phase_idx ].dh);
} }