OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-27 21:36:26 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added rs and rv output capabilities

Browse Source
This commit is contained in:
babrodtk 2016-09-02 14:15:10 +02:00
parent 40b2b95d87
commit 739976535f
14 changed files with 219 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
opm/autodiff/BlackoilModelBase.hpp
View File

@ -129,6 +129,13 @@ namespace Opm {
ADB mob; // Phase mobility (per cell)
};
struct SimulatorData {
SimulatorData(int num_phases);
std::vector<ReservoirResidualQuant> rq;
ADB rs;
ADB rv;
};
typedef typename ModelTraits<Implementation>::ReservoirState ReservoirState;
typedef typename ModelTraits<Implementation>::WellState WellState;
typedef typename ModelTraits<Implementation>::ModelParameters ModelParameters;
@ -272,6 +279,11 @@ namespace Opm {
WellModel& wellModel() { return well_model_; }
const WellModel& wellModel() const { return well_model_; }
/// Return reservoir simulation data (for output functionality)
const SimulatorData& getSimulatorData() const {
return sd_;
}
/// Compute fluid in place.
/// \param[in] ReservoirState
/// \param[in] FIPNUM for active cells not global cells.
@ -280,11 +292,6 @@ namespace Opm {
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum);
/// Return reservoir residual quantitites (in particular for output functionality)
const std::vector<ReservoirResidualQuant>& getReservoirResidualQuantities() const {
return rq_;
}
protected:
// --------- Types and enums ---------
@ -316,7 +323,7 @@ namespace Opm {
bool use_threshold_pressure_;
V threshold_pressures_by_connection_;
std::vector<ReservoirResidualQuant> rq_;
SimulatorData sd_;
std::vector<PhasePresence> phaseCondition_;
// Well Model

75
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -189,7 +189,7 @@ namespace detail {
, has_vapoil_(has_vapoil)
, param_( param )
, use_threshold_pressure_(false)
, rq_ (fluid.numPhases())
, sd_ (fluid.numPhases())
, phaseCondition_(AutoDiffGrid::numCells(grid))
, well_model_ (well_model)
, isRs_(V::Zero(AutoDiffGrid::numCells(grid)))
@ -465,6 +465,15 @@ namespace detail {
{
}
template <class Grid, class WellModel, class Implementation>
BlackoilModelBase<Grid, WellModel, Implementation>::
SimulatorData::SimulatorData(int num_phases)
: rq(num_phases)
, rv(ADB::null())
, rs(ADB::null())
{
}
@ -707,10 +716,10 @@ namespace detail {
for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
const int pos = pu.phase_pos[ phase ];
rq_[pos].b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos];
// OPM_AD_DUMP(rq_[pos].b);
// OPM_AD_DUMP(rq_[pos].accum[aix]);
sd_.rq[pos].b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
sd_.rq[pos].accum[aix] = pv_mult * sd_.rq[pos].b * sat[pos];
// OPM_AD_DUMP(sd_.rq[pos].b);
// OPM_AD_DUMP(sd_.rq[pos].accum[aix]);
}
}
@ -721,11 +730,11 @@ namespace detail {
// Temporary copy to avoid contribution of dissolved gas in the vaporized oil
// when both dissolved gas and vaporized oil are present.
const ADB accum_gas_copy =rq_[pg].accum[aix];
const ADB accum_gas_copy =sd_.rq[pg].accum[aix];
rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix];
rq_[po].accum[aix] += state.rv * accum_gas_copy;
// OPM_AD_DUMP(rq_[pg].accum[aix]);
sd_.rq[pg].accum[aix] += state.rs * sd_.rq[po].accum[aix];
sd_.rq[po].accum[aix] += state.rv * accum_gas_copy;
// OPM_AD_DUMP(sd_.rq[pg].accum[aix]);
}
}
@ -789,7 +798,7 @@ namespace detail {
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
asImpl().wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
@ -820,10 +829,10 @@ namespace detail {
{
// Compute b_p and the accumulation term b_p*s_p for each phase,
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
// These quantities are stored in rq_[phase].accum[1].
// These quantities are stored in sd_.rq[phase].accum[1].
// The corresponding accumulation terms from the start of
// the timestep (b^0_p*s^0_p etc.) were already computed
// on the initial call to assemble() and stored in rq_[phase].accum[0].
// on the initial call to assemble() and stored in sd_.rq[phase].accum[0].
asImpl().computeAccum(state, 1);
// Set up the common parts of the mass balance equations
@ -837,19 +846,19 @@ namespace detail {
{
const std::vector<ADB> kr = asImpl().computeRelPerm(state);
for (int phaseIdx=0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
rq_[phaseIdx].kr = kr[canph_[phaseIdx]];
sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
}
}
#pragma omp parallel for schedule(static)
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
const std::vector<PhasePresence>& cond = phaseCondition();
rq_[phaseIdx].mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
rq_[phaseIdx].rho = asImpl().fluidDensity(canph_[phaseIdx], rq_[phaseIdx].b, state.rs, state.rv);
asImpl().computeMassFlux(phaseIdx, trans_all, rq_[phaseIdx].kr, rq_[phaseIdx].mu, rq_[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
sd_.rq[phaseIdx].mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
sd_.rq[phaseIdx].rho = asImpl().fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
asImpl().computeMassFlux(phaseIdx, trans_all, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
residual_.material_balance_eq[ phaseIdx ] =
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ ops_.div*rq_[phaseIdx].mflux;
pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
+ ops_.div*sd_.rq[phaseIdx].mflux;
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
@ -862,15 +871,15 @@ namespace detail {
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
const UpwindSelector<double> upwindOil(grid_, ops_,
rq_[po].dh.value());
sd_.rq[po].dh.value());
const ADB rs_face = upwindOil.select(state.rs);
const UpwindSelector<double> upwindGas(grid_, ops_,
rq_[pg].dh.value());
sd_.rq[pg].dh.value());
const ADB rv_face = upwindGas.select(state.rv);
residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux);
residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
// OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
@ -897,7 +906,7 @@ namespace detail {
{
if (active_[idx]) {
const int pos = pu.phase_pos[idx];
const ADB& temp_b = rq_[pos].b;
const ADB& temp_b = sd_.rq[pos].b;
B = 1. / temp_b.value();
#if HAVE_MPI
if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
@ -1570,21 +1579,21 @@ namespace detail {
{
// Compute and store mobilities.
const ADB tr_mult = transMult(state.pressure);
rq_[ actph ].mob = tr_mult * kr / mu;
sd_.rq[ actph ].mob = tr_mult * kr / mu;
// Compute head differentials. Gravity potential is done using the face average as in eclipse and MRST.
const ADB rhoavg = ops_.caver * rho;
rq_[ actph ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
sd_.rq[ actph ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
if (use_threshold_pressure_) {
applyThresholdPressures(rq_[ actph ].dh);
applyThresholdPressures(sd_.rq[ actph ].dh);
}
// Compute phase fluxes with upwinding of formation value factor and mobility.
const ADB& b = rq_[ actph ].b;
const ADB& mob = rq_[ actph ].mob;
const ADB& dh = rq_[ actph ].dh;
const ADB& b = sd_.rq[ actph ].b;
const ADB& mob = sd_.rq[ actph ].mob;
const ADB& dh = sd_.rq[ actph ].dh;
UpwindSelector<double> upwind(grid_, ops_, dh.value());
rq_[ actph ].mflux = upwind.select(b * mob) * (transi * dh);
sd_.rq[ actph ].mflux = upwind.select(b * mob) * (transi * dh);
}
@ -1806,7 +1815,7 @@ namespace detail {
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = asImpl().numPhases();
const int nm = asImpl().numMaterials();
assert(int(rq_.size()) == nm);
assert(int(sd_.rq.size()) == nm);
const V& pv = geo_.poreVolume();
@ -1820,7 +1829,7 @@ namespace detail {
for ( int idx = 0; idx < nm; ++idx )
{
const ADB& tempB = rq_[idx].b;
const ADB& tempB = sd_.rq[idx].b;
B.col(idx) = 1./tempB.value();
R.col(idx) = residual_.material_balance_eq[idx].value();
tempV.col(idx) = R.col(idx).abs()/pv;
@ -1941,7 +1950,7 @@ namespace detail {
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
for ( int idx = 0; idx < nm; ++idx )
{
const ADB& tempB = rq_[idx].b;
const ADB& tempB = sd_.rq[idx].b;
B.col(idx) = 1./tempB.value();
R.col(idx) = residual_.material_balance_eq[idx].value();
tempV.col(idx) = R.col(idx).abs()/pv;

2
opm/autodiff/BlackoilMultiSegmentModel.hpp
View File

@ -123,7 +123,7 @@ namespace Opm {
using Base::pvdt_;
using Base::geo_;
using Base::active_;
using Base::rq_;
using Base::sd_;
using Base::fluid_;
using Base::terminal_output_;
using Base::grid_;

4
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
View File

@ -195,7 +195,7 @@ namespace Opm {
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
@ -273,7 +273,7 @@ namespace Opm {
// TODO: make sure the order of the density and the order of the kr are the same.
for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) {
const int canonicalPhaseIdx = canph_[phaseIdx];
fluid_density[phaseIdx] = fluidDensity(canonicalPhaseIdx, rq_[phaseIdx].b, state.rs, state.rv);
fluid_density[phaseIdx] = fluidDensity(canonicalPhaseIdx, sd_.rq[phaseIdx].b, state.rs, state.rv);
}
wellModel().computeWellConnectionPressures(state, well_state, kr_adb, fluid_density);
}

14
opm/autodiff/BlackoilPressureModel.hpp
View File

@ -136,7 +136,7 @@ namespace Opm {
using Base::asImpl;
using Base::linsolver_;
using Base::residual_;
using Base::rq_;
using Base::sd_;
using Base::grid_;
using Base::ops_;
using Base::has_vapoil_;
@ -167,11 +167,11 @@ namespace Opm {
residual_.material_balance_eq[0] = pressure_residual; // HACK
// Compute total reservoir volume flux.
const int n = rq_[0].mflux.size();
const int n = sd_.rq[0].mflux.size();
V flux = V::Zero(n);
for (int phase = 0; phase < numPhases(); ++phase) {
UpwindSelector<double> upwind(grid_, ops_, rq_[phase].dh.value());
flux += rq_[phase].mflux.value() / upwind.select(rq_[phase].b.value());
UpwindSelector<double> upwind(grid_, ops_, sd_.rq[phase].dh.value());
flux += sd_.rq[phase].mflux.value() / upwind.select(sd_.rq[phase].b.value());
}
// Storing the fluxes in the assemble() method is a bit of
@ -235,9 +235,9 @@ namespace Opm {
assert(numPhases() == 3);
assert(numMaterials() == 3);
V one = V::Constant(state.pressure.size(), 1.0);
scaling_[Water] = one / rq_[Water].b;
scaling_[Oil] = one / rq_[Oil].b - state.rs / rq_[Gas].b;
scaling_[Gas] = one / rq_[Gas].b - state.rv / rq_[Oil].b;
scaling_[Water] = one / sd_.rq[Water].b;
scaling_[Oil] = one / sd_.rq[Oil].b - state.rs / sd_.rq[Gas].b;
scaling_[Gas] = one / sd_.rq[Gas].b - state.rv / sd_.rq[Oil].b;
if (has_disgas_ && has_vapoil_) {
ADB r_factor = one / (one - state.rs * state.rv);
scaling_[Oil] = r_factor * scaling_[Oil];

10
opm/autodiff/BlackoilSequentialModel.hpp
View File

@ -53,8 +53,9 @@ namespace Opm {
typedef BlackoilSequentialModelParameters ModelParameters;
typedef DefaultBlackoilSolutionState SolutionState;
typedef BlackoilPressureModel<Grid, WellModel> PressureModel;
typedef BlackoilTransportModel<Grid, WellModel> TransportModel;
typedef typename TransportModel::ReservoirResidualQuant ReservoirResidualQuant;
typedef typename TransportModel::SimulatorData SimulatorData;
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
@ -263,14 +264,13 @@ namespace Opm {
}
/// Return reservoir residual quantitites (in particular for output functionality)
const std::vector<ReservoirResidualQuant>& getReservoirResidualQuantities() const {
return transport_model_->getReservoirResidualQuantities();
/// Return reservoir simulation data (for output functionality)
const SimulatorData& getSimulatorData() const {
return transport_model_->getSimulatorData();
}
protected:
typedef BlackoilPressureModel<Grid, WellModel> PressureModel;
typedef NonlinearSolver<PressureModel> PressureSolver;
typedef NonlinearSolver<TransportModel> TransportSolver;

2
opm/autodiff/BlackoilSolventModel.hpp
View File

@ -122,7 +122,7 @@ namespace Opm {
using Base::param_;
using Base::use_threshold_pressure_;
using Base::threshold_pressures_by_connection_;
using Base::rq_;
using Base::sd_;
using Base::phaseCondition_;
using Base::residual_;
using Base::terminal_output_;

14
opm/autodiff/BlackoilSolventModel_impl.hpp
View File

@ -95,7 +95,7 @@ namespace Opm {
if (has_solvent_) {
// If deck has solvent, residual_ should contain solvent equation.
rq_.resize(fluid_.numPhases() + 1);
sd_.rq.resize(fluid_.numPhases() + 1);
residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
Base::material_name_.push_back("Solvent");
assert(solvent_pos_ == fluid_.numPhases());
@ -270,8 +270,8 @@ namespace Opm {
const ADB& pg = state.canonical_phase_pressures[pu.phase_pos[Gas]];
const std::vector<PhasePresence>& cond = phaseCondition();
rq_[solvent_pos_].b = fluidReciprocFVF(Solvent, pg, state.temperature, state.rs, state.rv,cond);
rq_[solvent_pos_].accum[aix] = pv_mult * rq_[solvent_pos_].b * ss;
sd_.rq[solvent_pos_].b = fluidReciprocFVF(Solvent, pg, state.temperature, state.rs, state.rv,cond);
sd_.rq[solvent_pos_].accum[aix] = pv_mult * sd_.rq[solvent_pos_].b * ss;
}
}
@ -289,8 +289,8 @@ namespace Opm {
if (has_solvent_) {
residual_.material_balance_eq[ solvent_pos_ ] =
pvdt_ * (rq_[solvent_pos_].accum[1] - rq_[solvent_pos_].accum[0])
+ ops_.div*rq_[solvent_pos_].mflux;
pvdt_ * (sd_.rq[solvent_pos_].accum[1] - sd_.rq[solvent_pos_].accum[0])
+ ops_.div*sd_.rq[solvent_pos_].mflux;
}
}
@ -302,7 +302,7 @@ namespace Opm {
Base::updateEquationsScaling();
assert(MaxNumPhases + 1 == residual_.matbalscale.size());
if (has_solvent_) {
const ADB& temp_b = rq_[solvent_pos_].b;
const ADB& temp_b = sd_.rq[solvent_pos_].b;
ADB::V B = 1. / temp_b.value();
#if HAVE_MPI
if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
@ -660,7 +660,7 @@ namespace Opm {
// Compute solvent properties
const std::vector<PhasePresence>& cond = phaseCondition();
ADB mu_s = fluidViscosity(Solvent, phasePressure,state.temperature, state.rs, state.rv, cond);
ADB rho_s = fluidDensity(Solvent,rq_[solvent_pos_].b, state.rs, state.rv);
ADB rho_s = fluidDensity(Solvent,sd_.rq[solvent_pos_].b, state.rs, state.rv);
// Compute solvent relperm and mass flux
ADB krs = solvent_props_.solventRelPermMultiplier(F_solvent, cells_) * kr_mod;

44
opm/autodiff/BlackoilTransportModel.hpp
View File

@ -132,7 +132,7 @@ namespace Opm {
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
asImpl().wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
iter_report = asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state);
@ -203,7 +203,7 @@ namespace Opm {
using Base::linsolver_;
using Base::residual_;
using Base::rq_;
using Base::sd_;
using Base::geo_;
using Base::ops_;
using Base::grid_;
@ -281,10 +281,10 @@ namespace Opm {
{
// Compute b_p and the accumulation term b_p*s_p for each phase,
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
// These quantities are stored in rq_[phase].accum[1].
// These quantities are stored in sd_.rq[phase].accum[1].
// The corresponding accumulation terms from the start of
// the timestep (b^0_p*s^0_p etc.) were already computed
// on the initial call to assemble() and stored in rq_[phase].accum[0].
// on the initial call to assemble() and stored in sd_.rq[phase].accum[0].
asImpl().computeAccum(state, 1);
// Set up the common parts of the mass balance equations
@ -311,19 +311,19 @@ namespace Opm {
const ADB mu = asImpl().fluidViscosity(canonical_phase_idx, phase_pressure, state.temperature, state.rs, state.rv, cond);
// Note that the pressure-dependent transmissibility multipliers are considered
// part of the mobility here.
rq_[ phase_idx ].mob = tr_mult * kr[phase_idx] / mu;
sd_.rq[ phase_idx ].mob = tr_mult * kr[phase_idx] / mu;
// Compute head differentials. Gravity potential is done using the face average as in eclipse and MRST.
const ADB rho = asImpl().fluidDensity(canonical_phase_idx, rq_[phase_idx].b, state.rs, state.rv);
const ADB rho = asImpl().fluidDensity(canonical_phase_idx, sd_.rq[phase_idx].b, state.rs, state.rv);
const ADB rhoavg = ops_.caver * rho;
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) = -rq_[phase_idx].dh.value();
upwind_flags_.col(phase_idx) = -sd_.rq[phase_idx].dh.value();
}
if (use_threshold_pressure_) {
asImpl().applyThresholdPressures(rq_[ phase_idx ].dh);
asImpl().applyThresholdPressures(sd_.rq[ phase_idx ].dh);
}
}
@ -331,7 +331,7 @@ namespace Opm {
const ADB gradp = ops_.grad * state.pressure;
std::vector<ADB> dh_sat(numPhases(), ADB::null());
for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
dh_sat[phase_idx] = gradp - rq_[phase_idx].dh;
dh_sat[phase_idx] = gradp - sd_.rq[phase_idx].dh;
}
// Find upstream directions for each phase.
@ -346,9 +346,9 @@ namespace Opm {
ADB tot_mob = ADB::constant(V::Zero(gdz.size()));
for (int phase_idx = 0; phase_idx < numPhases(); ++phase_idx) {
UpwindSelector<double> upwind(grid_, ops_, upwind_flags_.col(phase_idx));
mob[phase_idx] = upwind.select(rq_[phase_idx].mob);
mob[phase_idx] = upwind.select(sd_.rq[phase_idx].mob);
tot_mob += mob[phase_idx];
b[phase_idx] = upwind.select(rq_[phase_idx].b);
b[phase_idx] = upwind.select(sd_.rq[phase_idx].b);
if (canph_[phase_idx] == Oil) {
rs = upwind.select(state.rs);
}
@ -365,7 +365,7 @@ namespace Opm {
gflux += mob[other_phase] * (dh_sat[phase_idx] - dh_sat[other_phase]);
}
}
rq_[phase_idx].mflux = b[phase_idx] * (mob[phase_idx] / tot_mob) * (total_flux_ + trans_all * gflux);
sd_.rq[phase_idx].mflux = b[phase_idx] * (mob[phase_idx] / tot_mob) * (total_flux_ + trans_all * gflux);
}
#pragma omp parallel for schedule(static)
@ -376,8 +376,8 @@ namespace Opm {
// Material balance equation for this phase.
residual_.material_balance_eq[ phase_idx ] =
pvdt_ * (rq_[phase_idx].accum[1] - rq_[phase_idx].accum[0])
+ ops_.div*rq_[phase_idx].mflux;
pvdt_ * (sd_.rq[phase_idx].accum[1] - sd_.rq[phase_idx].accum[0])
+ ops_.div*sd_.rq[phase_idx].mflux;
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
@ -388,8 +388,8 @@ namespace Opm {
if (active_[ Oil ] && active_[ Gas ]) {
const int po = fluid_.phaseUsage().phase_pos[ Oil ];
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
residual_.material_balance_eq[ pg ] += ops_.div * (rs * rq_[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv * rq_[pg].mflux);
residual_.material_balance_eq[ pg ] += ops_.div * (rs * sd_.rq[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv * sd_.rq[pg].mflux);
}
if (param_.update_equations_scaling_) {
@ -412,7 +412,7 @@ namespace Opm {
// Using the data members:
// total_flux_
// rq_[].mob
// sd_.rq[].mob
// Notation based on paper cited above.
const int num_connections = head_diff[0].size();
@ -441,10 +441,10 @@ namespace Opm {
theta[ell] = q;
for (int j = 0; j < num_phases; ++j) {
if (j < ell) {
theta[ell] += t * (g[ell].first - g[j].first) * rq_[g[j].second].mob.value()[b];
theta[ell] += t * (g[ell].first - g[j].first) * sd_.rq[g[j].second].mob.value()[b];
}
if (j > ell) {
theta[ell] += t * (g[ell].first - g[j].first) * rq_[g[j].second].mob.value()[a];
theta[ell] += t * (g[ell].first - g[j].first) * sd_.rq[g[j].second].mob.value()[a];
}
}
if (theta[ell] <= 0.0) {
@ -542,7 +542,7 @@ namespace Opm {
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = asImpl().numPhases();
const int nm = asImpl().numMaterials();
assert(int(rq_.size()) == nm);
assert(int(sd_.rq.size()) == nm);
const V& pv = geo_.poreVolume();
@ -556,7 +556,7 @@ namespace Opm {
for ( int idx = 0; idx < nm; ++idx )
{
const ADB& tempB = rq_[idx].b;
const ADB& tempB = sd_.rq[idx].b;
B.col(idx) = 1./tempB.value();
R.col(idx) = residual_.material_balance_eq[idx].value();
tempV.col(idx) = R.col(idx).abs()/pv;

45
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@ -25,6 +25,7 @@
#include <opm/core/utility/Compat.hpp>
#include <opm/core/utility/DataMap.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/output/eclipse/EclipseReader.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
@ -478,7 +479,7 @@ namespace Opm
outKeywords["VISC"] = std::max(outKeywords["VISC"], 1);
}
const std::vector<typename Model::ReservoirResidualQuant>& rq = model.getReservoirResidualQuantities();
const typename Model::SimulatorData& sd = model.getSimulatorData();
//Get shorthands for water, oil, gas
@ -498,19 +499,19 @@ namespace Opm
simProps.emplace_back(
"1OVERBW",
Opm::UnitSystem::measure::volume,
adbToDoubleVector(rq[aqua_idx].b));
adbToDoubleVector(sd.rq[aqua_idx].b));
}
if (liquid_active && outKeywords["BO"] > 0) {
simProps.emplace_back(
"1OVERBO",
Opm::UnitSystem::measure::volume,
adbToDoubleVector(rq[liquid_idx].b));
adbToDoubleVector(sd.rq[liquid_idx].b));
}
if (vapour_active && outKeywords["BG"] > 0) {
simProps.emplace_back(
"1OVERBG",
Opm::UnitSystem::measure::volume,
adbToDoubleVector(rq[vapour_idx].b));
adbToDoubleVector(sd.rq[vapour_idx].b));
}
/**
@ -521,19 +522,19 @@ namespace Opm
simProps.emplace_back(
"WAT_DEN",
Opm::UnitSystem::measure::density,
adbToDoubleVector(rq[aqua_idx].rho));
adbToDoubleVector(sd.rq[aqua_idx].rho));
}
if (liquid_active) {
simProps.emplace_back(
"OIL_DEN",
Opm::UnitSystem::measure::density,
adbToDoubleVector(rq[liquid_idx].rho));
adbToDoubleVector(sd.rq[liquid_idx].rho));
}
if (vapour_active) {
simProps.emplace_back(
"GAS_DEN",
Opm::UnitSystem::measure::density,
adbToDoubleVector(rq[vapour_idx].rho));
adbToDoubleVector(sd.rq[vapour_idx].rho));
}
}
@ -545,19 +546,19 @@ namespace Opm
simProps.emplace_back(
"WAT_VISC",
Opm::UnitSystem::measure::viscosity,
adbToDoubleVector(rq[aqua_idx].mu));
adbToDoubleVector(sd.rq[aqua_idx].mu));
}
if (liquid_active) {
simProps.emplace_back(
"OIL_VISC",
Opm::UnitSystem::measure::viscosity,
adbToDoubleVector(rq[liquid_idx].mu));
adbToDoubleVector(sd.rq[liquid_idx].mu));
}
if (vapour_active) {
simProps.emplace_back(
"GAS_VISC",
Opm::UnitSystem::measure::viscosity,
adbToDoubleVector(rq[vapour_idx].mu));
adbToDoubleVector(sd.rq[vapour_idx].mu));
}
}
@ -568,29 +569,45 @@ namespace Opm
simProps.emplace_back(
"WATKR",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(rq[aqua_idx].kr));
adbToDoubleVector(sd.rq[aqua_idx].kr));
}
if (aqua_active && outKeywords["KRO"] > 0) {
simProps.emplace_back(
"OILKR",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(rq[liquid_idx].kr));
adbToDoubleVector(sd.rq[liquid_idx].kr));
}
if (aqua_active && outKeywords["KRG"] > 0) {
simProps.emplace_back(
"GASKR",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(rq[vapour_idx].kr));
adbToDoubleVector(sd.rq[vapour_idx].kr));
}
/**
* Vaporized and dissolved gas/oil ratio
*/
if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) {
//FIXME: This requires a separate structure instead of RQ. Perhaps solutionstate?
simProps.emplace_back(
"RVSAT",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(sd.rv));
}
if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) {
simProps.emplace_back(
"RSSAT",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(sd.rs));
}
/**
* Bubble point and dew point pressures
*/
if (vapour_active && liquid_active && outKeywords["PBPD"] > 0) {
Opm::OpmLog::warning("Bubble/dew point pressure output unsupported",
"Writing bubble points and dew points (PBPD) to file is unsupported, "
"as the simulator does not use these internally.");
}

2
opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
View File

@ -172,7 +172,7 @@ namespace Opm {
using Base::param_;
using Base::use_threshold_pressure_;
using Base::threshold_pressures_by_connection_;
using Base::rq_;
using Base::sd_;
using Base::phaseCondition_;
using Base::residual_;
using Base::terminal_output_;

68
opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
View File

@ -109,7 +109,7 @@ namespace Opm {
}
residual_.matbalscale.resize(fluid_.numPhases() + 1, 1.1169); // use the same as the water phase
// If deck has polymer, residual_ should contain polymer equation.
rq_.resize(fluid_.numPhases() + 1);
sd_.rq.resize(fluid_.numPhases() + 1);
residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
Base::material_name_.push_back("Polymer");
assert(poly_pos_ == fluid_.numPhases());
@ -245,7 +245,7 @@ namespace Opm {
const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], AutoDiffGrid::numCells(grid_));
const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
// Compute polymer accumulation term.
rq_[poly_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
sd_.rq[poly_pos_].accum[aix] = pv_mult * sd_.rq[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol)
+ pv_mult * rho_rock * (1. - phi) / phi * ads;
}
@ -282,10 +282,10 @@ namespace Opm {
// Compute b_p and the accumulation term b_p*s_p for each phase,
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
// These quantities are stored in rq_[phase].accum[1].
// These quantities are stored in sd_.rq[phase].accum[1].
// The corresponding accumulation terms from the start of
// the timestep (b^0_p*s^0_p etc.) were already computed
// on the initial call to assemble() and stored in rq_[phase].accum[0].
// on the initial call to assemble() and stored in sd_.rq[phase].accum[0].
computeAccum(state, 1);
// Set up the common parts of the mass balance equations
@ -308,12 +308,12 @@ namespace Opm {
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
const std::vector<PhasePresence>& cond = phaseCondition();
const ADB mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
const ADB rho = fluidDensity(canph_[phaseIdx], rq_[phaseIdx].b, state.rs, state.rv);
const ADB rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], mu, rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
residual_.material_balance_eq[ phaseIdx ] =
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])
+ ops_.div*rq_[phaseIdx].mflux;
pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
+ ops_.div*sd_.rq[phaseIdx].mflux;
}
// -------- Extra (optional) rs and rv contributions to the mass balance equations --------
@ -326,15 +326,15 @@ namespace Opm {
const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
const UpwindSelector<double> upwindOil(grid_, ops_,
rq_[po].dh.value());
sd_.rq[po].dh.value());
const ADB rs_face = upwindOil.select(state.rs);
const UpwindSelector<double> upwindGas(grid_, ops_,
rq_[pg].dh.value());
sd_.rq[pg].dh.value());
const ADB rv_face = upwindGas.select(state.rv);
residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux);
residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
// OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
@ -342,8 +342,8 @@ namespace Opm {
// Add polymer equation.
if (has_polymer_) {
residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
+ ops_.div*rq_[poly_pos_].mflux;
residual_.material_balance_eq[poly_pos_] = pvdt_ * (sd_.rq[poly_pos_].accum[1] - sd_.rq[poly_pos_].accum[0])
+ ops_.div*sd_.rq[poly_pos_].mflux;
}
@ -462,24 +462,24 @@ namespace Opm {
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr);
const ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
// Reduce mobility of water phase by relperm reduction and effective viscosity increase.
rq_[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
sd_.rq[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
// Compute polymer mobility.
const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value());
rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_poly_eff_visc;
rq_[poly_pos_].b = rq_[actph].b;
rq_[poly_pos_].dh = rq_[actph].dh;
UpwindSelector<double> upwind(grid_, ops_, rq_[poly_pos_].dh.value());
sd_.rq[poly_pos_].mob = tr_mult * mc * krw_eff * inv_poly_eff_visc;
sd_.rq[poly_pos_].b = sd_.rq[actph].b;
sd_.rq[poly_pos_].dh = sd_.rq[actph].dh;
UpwindSelector<double> upwind(grid_, ops_, sd_.rq[poly_pos_].dh.value());
// Compute polymer flux.
rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * (transi * rq_[poly_pos_].dh);
sd_.rq[poly_pos_].mflux = upwind.select(sd_.rq[poly_pos_].b * sd_.rq[poly_pos_].mob) * (transi * sd_.rq[poly_pos_].dh);
// Must recompute water flux since we have to use modified mobilities.
rq_[ actph ].mflux = upwind.select(rq_[actph].b * rq_[actph].mob) * (transi * rq_[actph].dh);
sd_.rq[ actph ].mflux = upwind.select(sd_.rq[actph].b * sd_.rq[actph].mob) * (transi * sd_.rq[actph].dh);
// applying the shear-thinning factors
if (has_plyshlog_) {
V shear_mult_faces_v = Eigen::Map<V>(shear_mult_faces_.data(), shear_mult_faces_.size());
ADB shear_mult_faces_adb = ADB::constant(shear_mult_faces_v);
rq_[poly_pos_].mflux = rq_[poly_pos_].mflux / shear_mult_faces_adb;
rq_[actph].mflux = rq_[actph].mflux / shear_mult_faces_adb;
sd_.rq[poly_pos_].mflux = sd_.rq[poly_pos_].mflux / shear_mult_faces_adb;
sd_.rq[actph].mflux = sd_.rq[actph].mflux / shear_mult_faces_adb;
}
}
}
@ -534,7 +534,7 @@ namespace Opm {
std::vector<ADB> mob_perfcells;
std::vector<ADB> b_perfcells;
wellModel().extractWellPerfProperties(state, rq_, mob_perfcells, b_perfcells);
wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
if (param_.solve_welleq_initially_ && initial_assembly) {
// solve the well equations as a pre-processing step
Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state);
@ -598,19 +598,19 @@ namespace Opm {
const ADB tr_mult = transMult(state.pressure);
const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond);
rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
sd_.rq[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rho = fluidDensity(canonicalPhaseIdx, rq_[phase].b, state.rs, state.rv);
const ADB rho = fluidDensity(canonicalPhaseIdx, sd_.rq[phase].b, state.rs, state.rv);
const ADB rhoavg = ops_.caver * rho;
rq_[ phase ].dh = ops_.ngrad * phasePressure[ canonicalPhaseIdx ] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
sd_.rq[ phase ].dh = ops_.ngrad * phasePressure[ canonicalPhaseIdx ] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
if (use_threshold_pressure_) {
applyThresholdPressures(rq_[ phase ].dh);
applyThresholdPressures(sd_.rq[ phase ].dh);
}
const ADB& b = rq_[ phase ].b;
const ADB& mob = rq_[ phase ].mob;
const ADB& dh = rq_[ phase ].dh;
const ADB& b = sd_.rq[ phase ].b;
const ADB& mob = sd_.rq[ phase ].mob;
const ADB& dh = sd_.rq[ phase ].dh;
UpwindSelector<double> upwind(grid_, ops_, dh.value());
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
@ -619,7 +619,7 @@ namespace Opm {
cmax,
kr[canonicalPhaseIdx]);
ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
rq_[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
sd_.rq[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;
const V& polymer_conc = state.concentration.value();
V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc);
@ -629,7 +629,7 @@ namespace Opm {
visc_mult.resize(nface);
std::copy(visc_mult_faces.data(), visc_mult_faces.data() + nface, visc_mult.begin());
rq_[ phase ].mflux = (transi * upwind.select(b * mob)) * dh;
sd_.rq[ phase ].mflux = (transi * upwind.select(b * mob)) * dh;
const auto& b_faces_adb = upwind.select(b);
@ -650,7 +650,7 @@ namespace Opm {
std::vector<double> phiavg(phiavg_adb.value().data(), phiavg_adb.value().data() + phiavg_adb.size());
water_vel.resize(nface);
std::copy(rq_[0].mflux.value().data(), rq_[0].mflux.value().data() + nface, water_vel.begin());
std::copy(sd_.rq[0].mflux.value().data(), sd_.rq[0].mflux.value().data() + nface, water_vel.begin());
for (size_t i = 0; i < nface; ++i) {
water_vel[i] = water_vel[i] / (b_faces[i] * phiavg[i] * internal_face_areas[i]);
@ -723,7 +723,7 @@ namespace Opm {
std::copy(visc_mult_wells_v.data(), visc_mult_wells_v.data() + visc_mult_wells_v.size(), visc_mult_wells.begin());
const int water_pos = fluid_.phaseUsage().phase_pos[Water];
ADB b_perfcells = subset(rq_[water_pos].b, well_cells);
ADB b_perfcells = subset(sd_.rq[water_pos].b, well_cells);
const ADB& p_perfcells = subset(state.pressure, well_cells);
const V& cdp = wellModel().wellPerforationPressureDiffs();

76
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@ -185,7 +185,7 @@ namespace {
, grav_ (gravityOperator(grid_, ops_, geo_))
, cmax_(V::Zero(grid.number_of_cells))
, phaseCondition_ (grid.number_of_cells)
, rq_ (fluid.numPhases() + 1)
, sd_ (fluid.numPhases() + 1)
, residual_ ( { std::vector<ADB>(fluid.numPhases() + 1, ADB::null()),
ADB::null(),
ADB::null(),
@ -519,17 +519,17 @@ namespace {
const ADB pv_mult = poroMult(press);
for (int phase = 0; phase < 2; ++phase) {
rq_[phase].b = fluidReciprocFVF(phase, pressure[phase], temp, cond, cells_);
sd_.rq[phase].b = fluidReciprocFVF(phase, pressure[phase], temp, cond, cells_);
}
rq_[0].accum[aix] = pv_mult * rq_[0].b * sat[0];
rq_[1].accum[aix] = pv_mult * rq_[1].b * sat[1];
sd_.rq[0].accum[aix] = pv_mult * sd_.rq[0].b * sat[0];
sd_.rq[1].accum[aix] = pv_mult * sd_.rq[1].b * sat[1];
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
const double rho_rock = polymer_props_ad_.rockDensity();
const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], grid_.number_of_cells, 1);
const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
rq_[2].accum[aix] = pv_mult * rq_[0].b * sat[0] * c * (1. - dead_pore_vol) + pv_mult * rho_rock * (1. - phi) / phi * ads;
sd_.rq[2].accum[aix] = pv_mult * sd_.rq[0].b * sat[0] * c * (1. - dead_pore_vol) + pv_mult * rho_rock * (1. - phi) / phi * ads;
}
@ -638,10 +638,10 @@ namespace {
// Compute b_p and the accumulation term b_p*s_p for each phase,
// except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o.
// These quantities are stored in rq_[phase].accum[1].
// These quantities are stored in sd_.rq[phase].accum[1].
// The corresponding accumulation terms from the start of
// the timestep (b^0_p*s^0_p etc.) were already computed
// in step() and stored in rq_[phase].accum[0].
// in step() and stored in sd_.rq[phase].accum[0].
computeAccum(state, 1);
// Set up the common parts of the mass balance equations
@ -651,19 +651,19 @@ namespace {
const std::vector<ADB> kr = computeRelPerm(state);
const auto& pu = fluid_.phaseUsage();
for (int phaseIdx=0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
rq_[phaseIdx].kr = kr[pu.phase_pos[phaseIdx]];
sd_.rq[phaseIdx].kr = kr[pu.phase_pos[phaseIdx]];
}
}
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, rq_[0].kr);
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, sd_.rq[0].kr);
const ADB mc = computeMc(state);
computeMassFlux(trans, mc, rq_[1].kr, krw_eff, state);
residual_.material_balance_eq[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0])
+ ops_.div*rq_[0].mflux;
residual_.material_balance_eq[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0])
+ ops_.div*rq_[1].mflux;
residual_.material_balance_eq[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0]) //+ cell / dt * (rq_[2].ads[1] - rq_[2].ads[0])
+ ops_.div*rq_[2].mflux;
computeMassFlux(trans, mc, sd_.rq[1].kr, krw_eff, state);
residual_.material_balance_eq[0] = pvdt*(sd_.rq[0].accum[1] - sd_.rq[0].accum[0])
+ ops_.div*sd_.rq[0].mflux;
residual_.material_balance_eq[1] = pvdt*(sd_.rq[1].accum[1] - sd_.rq[1].accum[0])
+ ops_.div*sd_.rq[1].mflux;
residual_.material_balance_eq[2] = pvdt*(sd_.rq[2].accum[1] - sd_.rq[2].accum[0]) //+ cell / dt * (sd_.rq[2].ads[1] - sd_.rq[2].ads[0])
+ ops_.div*sd_.rq[2].mflux;
// -------- Extra (optional) sg or rs equation, and rs contributions to the mass balance equations --------
@ -735,13 +735,13 @@ namespace {
// DUMP(nkgradp_well);
const Selector<double> cell_to_well_selector(nkgradp_well.value());
ADB well_rates_all = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern());
ADB perf_total_mob = subset(rq_[0].mob, well_cells) + subset(rq_[1].mob, well_cells);
ADB perf_total_mob = subset(sd_.rq[0].mob, well_cells) + subset(sd_.rq[1].mob, well_cells);
std::vector<ADB> well_contribs(np, ADB::null());
std::vector<ADB> well_perf_rates(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
const ADB& cell_b = rq_[phase].b;
const ADB& cell_b = sd_.rq[phase].b;
const ADB perf_b = subset(cell_b, well_cells);
const ADB& cell_mob = rq_[phase].mob;
const ADB& cell_mob = sd_.rq[phase].mob;
const V well_fraction = compi.col(phase);
// Using total mobilities for all phases for injection.
const ADB perf_mob_injector = (wops_.w2p * well_fraction.matrix()).array() * perf_total_mob;
@ -959,37 +959,37 @@ namespace {
{
const ADB tr_mult = transMult(state.pressure);
const std::vector<PhasePresence> cond = phaseCondition();
std::vector<ADB> press = computePressures(state);
const ADB& temp = state.temperature;
std::vector<ADB> press = computePressures(state);
const ADB& temp = state.temperature;
{
{
const ADB mu_w = fluidViscosity(0, press[0], temp, cond, cells_);
ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu_w.value());
rq_[0].mob = tr_mult * krw_eff * inv_wat_eff_vis;
rq_[2].mob = tr_mult * mc * krw_eff * inv_wat_eff_vis;
sd_.rq[0].mob = tr_mult * krw_eff * inv_wat_eff_vis;
sd_.rq[2].mob = tr_mult * mc * krw_eff * inv_wat_eff_vis;
const ADB mu_o = fluidViscosity(1, press[1], temp, cond, cells_);
rq_[1].mob = tr_mult * kro / mu_o;
sd_.rq[1].mob = tr_mult * kro / mu_o;
rq_[0].mu = mu_w;
rq_[1].mu = mu_o;
}
sd_.rq[0].mu = mu_w;
sd_.rq[1].mu = mu_o;
}
for (int phase = 0; phase < 2; ++phase) {
rq_[phase].rho = fluidDensity(phase, press[phase], temp, cond, cells_);
ADB& head = rq_[phase].head;
sd_.rq[phase].rho = fluidDensity(phase, press[phase], temp, cond, cells_);
ADB& head = sd_.rq[phase].head;
// compute gravity potensial using the face average as in eclipse and MRST
const ADB rhoavg = ops_.caver * rq_[phase].rho;
const ADB rhoavg = ops_.caver * sd_.rq[phase].rho;
const ADB dp = ops_.ngrad * press[phase] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
head = transi*dp;
UpwindSelector<double> upwind(grid_, ops_, head.value());
const ADB& b = rq_[phase].b;
const ADB& mob = rq_[phase].mob;
rq_[phase].mflux = upwind.select(b * mob) * head;
const ADB& b = sd_.rq[phase].b;
const ADB& mob = sd_.rq[phase].mob;
sd_.rq[phase].mflux = upwind.select(b * mob) * head;
}
rq_[2].b = rq_[0].b;
rq_[2].head = rq_[0].head;
UpwindSelector<double> upwind(grid_, ops_, rq_[2].head.value());
rq_[2].mflux = upwind.select(rq_[2].b * rq_[2].mob) * rq_[2].head;
sd_.rq[2].b = sd_.rq[0].b;
sd_.rq[2].head = sd_.rq[0].head;
UpwindSelector<double> upwind(grid_, ops_, sd_.rq[2].head.value());
sd_.rq[2].mflux = upwind.select(sd_.rq[2].b * sd_.rq[2].mob) * sd_.rq[2].head;
}

19
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp
View File

@ -75,6 +75,17 @@ namespace Opm {
std::vector<ADB> ads; // Adsorption term.
};
struct SimulatorData {
SimulatorData(int num_phases)
: rq(num_phases)
, rs(ADB::null())
, rv(ADB::null())
{
}
std::vector<ReservoirResidualQuant> rq;
ADB rs;
ADB rv;
};
/// Construct a solver. It will retain references to the
/// arguments of this functions, and they are expected to
@ -124,9 +135,9 @@ namespace Opm {
double relativeChange(const PolymerBlackoilState& previous,
const PolymerBlackoilState& current ) const;
/// Return reservoir residual quantitites (in particular for output functionality)
const std::vector<ReservoirResidualQuant>& getReservoirResidualQuantities() const {
return rq_;
/// Return reservoir simulation data (for output functionality)
const SimulatorData& getSimulatorData() const {
return sd_;
}
/// Compute fluid in place.
@ -173,7 +184,7 @@ namespace Opm {
const M grav_;
V cmax_;
std::vector<PhasePresence> phaseCondition_;
std::vector<ReservoirResidualQuant> rq_;
SimulatorData sd_;
// The mass_balance vector has one element for each active phase,
// each of which has size equal to the number of cells.
// The well_eq has size equal to the number of wells.