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

Merge remote-tracking branch 'origin/master' into frankenstein

Browse Source
This commit is contained in:
Andreas Lauser 2016-09-12 23:18:02 +02:00
commit 62de30d9b2
35 changed files with 951 additions and 273 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/opm_init_check.cpp
View File

@ -208,7 +208,7 @@ public:
/*****************************************************************/
void initOPMTrans(TransGraph& opmTrans, DeckConstPtr deck, std::shared_ptr<const EclipseState> eclipseState) {
std::shared_ptr<GridManager> grid = std::make_shared<GridManager>( eclipseState->getInputGrid(),
std::shared_ptr<GridManager> grid = std::make_shared<GridManager>( *eclipseState->getInputGrid(),
eclipseState->get3DProperties().getDoubleGridProperty( "PORV" ).getData() );
const struct UnstructuredGrid * cGrid = grid->c_grid();
std::shared_ptr<BlackoilPropsAdInterface> props;

2
examples/sim_2p_comp_reorder.cpp
View File

@ -106,7 +106,7 @@ try
eclipseState.reset(new EclipseState(*deck, parseContext));
// Grid init
grid.reset(new GridManager(eclipseState->getInputGrid()));
grid.reset(new GridManager(*eclipseState->getInputGrid()));
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
state.reset( new BlackoilState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) ,2));

2
examples/sim_2p_incomp.cpp
View File

@ -117,7 +117,7 @@ try
deck = parser->parseFile(deck_filename , parseContext);
eclipseState.reset( new EclipseState(*deck, parseContext));
// Grid init
grid.reset(new GridManager(eclipseState->getInputGrid()));
grid.reset(new GridManager(*eclipseState->getInputGrid()));
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init

2
examples/sim_2p_incomp_ad.cpp
View File

@ -118,7 +118,7 @@ try
eclipseState.reset(new EclipseState(*deck , parseContext));
// Grid init
grid.reset(new GridManager(eclipseState->getInputGrid()));
grid.reset(new GridManager(*eclipseState->getInputGrid()));
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());

2
examples/sim_poly2p_comp_reorder.cpp
View File

@ -105,7 +105,7 @@ try
eclipseState.reset(new Opm::EclipseState(*deck , parseContext));
// Grid init
grid.reset(new GridManager(eclipseState->getInputGrid()));
grid.reset(new GridManager(*eclipseState->getInputGrid()));
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());

2
examples/sim_poly2p_incomp_reorder.cpp
View File

@ -106,7 +106,7 @@ try
eclipseState.reset(new Opm::EclipseState(*deck , parseContext));
// Grid init
grid.reset(new GridManager(eclipseState->getInputGrid()));
grid.reset(new GridManager(*eclipseState->getInputGrid()));
{
const UnstructuredGrid& ug_grid = *(grid->c_grid());
// Rock and fluid init

6
examples/sim_poly_fi2p_comp_ad.cpp
View File

@ -168,9 +168,9 @@ try
if (eclipseState->get3DProperties().hasDeckDoubleGridProperty("PORV")) {
const auto& porv = eclipseState->get3DProperties().getDoubleGridProperty("PORV").getData();
grid.reset(new GridManager(eclipseState->getInputGrid(), porv));
grid.reset(new GridManager(*eclipseState->getInputGrid(), porv));
} else {
grid.reset(new GridManager(eclipseState->getInputGrid()));
grid.reset(new GridManager(*eclipseState->getInputGrid()));
}
auto &cGrid = *grid->c_grid();
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
@ -240,7 +240,7 @@ try
<< std::flush;
Opm::BlackoilOutputWriter
outputWriter(cGrid, param, eclipseState, Opm::NNC(), pu,
outputWriter(cGrid, param, eclipseState, pu,
new_props->permeability() );
SimulatorReport fullReport;

42
opm/autodiff/BlackoilModelBase.hpp
View File

@ -78,6 +78,25 @@ namespace Opm {
typedef ADB::V V;
typedef ADB::M M;
struct ReservoirResidualQuant {
ReservoirResidualQuant();
std::vector<ADB> accum; // Accumulations
ADB mflux; // Mass flux (surface conditions)
ADB b; // Reciprocal FVF
ADB mu; // Viscosities
ADB rho; // Densities
ADB kr; // Permeabilities
ADB dh; // Pressure drop across int. interfaces
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;
@ -221,6 +240,19 @@ 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.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum);
protected:
// --------- Types and enums ---------
@ -230,14 +262,6 @@ namespace Opm {
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
struct ReservoirResidualQuant {
ReservoirResidualQuant();
std::vector<ADB> accum; // Accumulations
ADB mflux; // Mass flux (surface conditions)
ADB b; // Reciprocal FVF
ADB dh; // Pressure drop across int. interfaces
ADB mob; // Phase mobility (per cell)
};
// --------- Data members ---------
@ -260,7 +284,7 @@ namespace Opm {
bool use_threshold_pressure_;
V threshold_pressures_by_connection_;
std::vector<ReservoirResidualQuant> rq_;
SimulatorData sd_;
std::vector<PhasePresence> phaseCondition_;
// Well Model

161
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -123,7 +123,7 @@ typedef Eigen::Array<double,
, 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)))
@ -391,11 +391,23 @@ typedef Eigen::Array<double,
: accum(2, ADB::null())
, mflux( ADB::null())
, b ( ADB::null())
, mu ( ADB::null())
, rho ( ADB::null())
, kr ( ADB::null())
, dh ( ADB::null())
, mob ( ADB::null())
{
}
template <class Grid, class WellModel, class Implementation>
BlackoilModelBase<Grid, WellModel, Implementation>::
SimulatorData::SimulatorData(int num_phases)
: rq(num_phases)
, rs(ADB::null())
, rv(ADB::null())
{
}
@ -638,10 +650,10 @@ typedef Eigen::Array<double,
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]);
}
}
@ -652,11 +664,11 @@ typedef Eigen::Array<double,
// 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]);
}
}
@ -720,7 +732,7 @@ typedef Eigen::Array<double,
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);
@ -751,10 +763,10 @@ typedef Eigen::Array<double,
{
// 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
@ -764,17 +776,23 @@ typedef Eigen::Array<double,
V trans_all(transi.size() + trans_nnc.size());
trans_all << transi, trans_nnc;
const std::vector<ADB> kr = asImpl().computeRelPerm(state);
{
const std::vector<ADB> kr = asImpl().computeRelPerm(state);
for (int phaseIdx=0; phaseIdx < fluid_.numPhases(); ++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();
const ADB mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
const ADB rho = asImpl().fluidDensity(canph_[phaseIdx], rq_[phaseIdx].b, state.rs, state.rv);
asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], mu, 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 --------
@ -787,15 +805,15 @@ typedef Eigen::Array<double,
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 ]);
@ -822,7 +840,7 @@ typedef Eigen::Array<double,
{
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) )
@ -1366,21 +1384,21 @@ typedef Eigen::Array<double,
{
// 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);
}
@ -1602,7 +1620,7 @@ typedef Eigen::Array<double,
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();
@ -1616,7 +1634,7 @@ typedef Eigen::Array<double,
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;
@ -1737,7 +1755,7 @@ typedef Eigen::Array<double,
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;
@ -2100,6 +2118,83 @@ typedef Eigen::Array<double,
template <class Grid, class WellModel, class Implementation>
std::vector<V>
BlackoilModelBase<Grid, WellModel, Implementation>::
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum)
{
using namespace Opm::AutoDiffGrid;
const int nc = numCells(grid_);
std::vector<ADB> saturation(3, ADB::null());
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, x.numPhases());
const ADB pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
const ADB temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
saturation[Water] = active_[Water] ? ADB::constant(s.col(Water)) : ADB::null();
saturation[Oil] = active_[Oil] ? ADB::constant(s.col(Oil)) : ADB::constant(V::Zero(nc));
saturation[Gas] = active_[Gas] ? ADB::constant(s.col(Gas)) : ADB::constant(V::Zero(nc));
const ADB rs = ADB::constant(Eigen::Map<const V>(& x.gasoilratio()[0], nc, 1));
const ADB rv = ADB::constant(Eigen::Map<const V>(& x.rv()[0], nc, 1));
const auto canonical_phase_pressures = computePressures(pressure, saturation[Water], saturation[Oil], saturation[Gas]);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
const std::vector<PhasePresence> cond = phaseCondition();
const ADB pv_mult = poroMult(pressure);
const V& pv = geo_.poreVolume();
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
std::vector<V> fip(5, V::Zero(nc));
for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
const int pos = pu.phase_pos[ phase ];
const auto& b = asImpl().fluidReciprocFVF(phase, canonical_phase_pressures[phase], temperature, rs, rv, cond);
fip[phase] = ((pv_mult * b * saturation[pos] * pv).value());
}
}
if (active_[ Oil ] && active_[ Gas ]) {
// Account for gas dissolved in oil and vaporized oil
const int po = pu.phase_pos[Oil];
const int pg = pu.phase_pos[Gas];
fip[3] = rs.value() * fip[po];
fip[4] = rv.value() * fip[pg];
}
const int dims = *std::max_element(fipnum.begin(), fipnum.end());
std::vector<V> values(dims, V::Zero(7));
for (int i = 0; i < 5; ++i) {
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) {
values[fipnum[c]-1][i] += fip[i][c];
}
}
}
// compute PAV and PORV for every regions.
const V hydrocarbon = saturation[Oil].value() + saturation[Gas].value();
V hcpv = V::Zero(nc);
V pres = V::Zero(nc);
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) {
hcpv[fipnum[c]-1] += pv[c] * hydrocarbon[c];
pres[fipnum[c]-1] += pv[c] * pressure.value()[c];
values[fipnum[c]-1][5] += pv[c];
values[fipnum[c]-1][6] += pv[c] * pressure.value()[c] * hydrocarbon[c];
}
}
for (int reg = 0; reg < dims; ++reg) {
if (hcpv[reg] != 0) {
values[reg][6] /= hcpv[reg];
} else {
values[reg][6] = pres[reg] / values[reg][5];
}
}
return values;
}
} // namespace Opm
#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED

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

@ -137,7 +137,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_;
@ -168,11 +168,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
@ -236,9 +236,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];

21
opm/autodiff/BlackoilSequentialModel.hpp
View File

@ -53,6 +53,10 @@ namespace Opm {
typedef BlackoilSequentialModelParameters ModelParameters;
typedef DefaultBlackoilSolutionState SolutionState;
typedef BlackoilPressureModel<Grid, WellModel> PressureModel;
typedef BlackoilTransportModel<Grid, WellModel> TransportModel;
typedef typename TransportModel::SimulatorData SimulatorData;
/// Construct the model. It will retain references to the
/// arguments of this functions, and they are expected to
/// remain in scope for the lifetime of the solver.
@ -247,13 +251,26 @@ namespace Opm {
}
/// Compute fluid in place.
/// \param[in] ReservoirState
/// \param[in] WellState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum) const
{
return transport_solver_.computeFluidInPlace(x, fipnum);
}
/// Return reservoir simulation data (for output functionality)
const SimulatorData& getSimulatorData() const {
return transport_model_->getSimulatorData();
}
protected:
typedef BlackoilPressureModel<Grid, WellModel> PressureModel;
typedef BlackoilTransportModel<Grid, WellModel> TransportModel;
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;

50
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);
@ -160,8 +160,6 @@ namespace Opm {
/// Solve the Jacobian system Jx = r where J is the Jacobian and
/// r is the residual.
V solveJacobianSystem() const
@ -205,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_;
@ -275,14 +273,18 @@ namespace Opm {
void assembleMassBalanceEq(const SolutionState& state)
{
// 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
@ -309,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);
}
}
@ -329,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.
@ -344,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);
}
@ -363,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)
@ -374,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 --------
@ -386,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_) {
@ -410,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();
@ -439,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) {
@ -540,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();
@ -554,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;

17
opm/autodiff/FlowMain.hpp
View File

@ -145,6 +145,7 @@ namespace Opm
asImpl().extractMessages();
asImpl().runDiagnostics();
asImpl().setupState();
asImpl().writeInit();
asImpl().distributeData();
asImpl().setupOutputWriter();
asImpl().setupLinearSolver();
@ -649,7 +650,19 @@ namespace Opm
}
void writeInit()
{
bool output = param_.getDefault("output", true);
bool output_ecl = param_.getDefault("output_ecl", true);
const Grid& grid = grid_init_->grid();
if( output && output_ecl && output_cout_)
{
const EclipseGrid& inputGrid = *eclipse_state_->getInputGrid();
EclipseWriter writer(eclipse_state_, UgGridHelpers::createEclipseGrid( grid , inputGrid ));
writer.writeInitAndEgrid(geoprops_->simProps(grid),
geoprops_->nonCartesianConnections());
}
}
// Setup output writer.
@ -663,7 +676,6 @@ namespace Opm
output_writer_.reset(new BlackoilOutputWriter(grid_init_->grid(),
param_,
eclipse_state_,
geoprops_->nonCartesianConnections(),
Opm::phaseUsageFromDeck(deck_),
fluidprops_->permeability()));
}
@ -747,7 +759,6 @@ namespace Opm
fullReport.reportParam(tot_os);
}
} else {
output_writer_->writeInit(geoprops_->simProps(grid_init_->grid()) , geoprops_->nonCartesianConnections( ));
if (output_cout_) {
std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
}

3
opm/autodiff/GeoProps.hpp
View File

@ -382,6 +382,7 @@ namespace Opm
const int* cartdims = Opm::UgGridHelpers::cartDims(grid);
EclipseGridConstPtr eclgrid = eclState->getInputGrid();
const auto& porv = eclState->get3DProperties().getDoubleGridProperty("PORV").getData();
const auto& actnum = eclState->get3DProperties().getIntGridProperty("ACTNUM").getData();
for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
const int nx = cartdims[0];
const int ny = cartdims[1];
@ -395,7 +396,7 @@ namespace Opm
// that has pore volume less than the MINPV threshold
int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
while ( cartesianCellIdxAbove >= 0 &&
porv[cartesianCellIdxAbove] > 0 &&
actnum[cartesianCellIdxAbove] > 0 &&
porv[cartesianCellIdxAbove] < eclgrid->getMinpvValue() ) {
// Volume weighted arithmetic average of NTG

4
opm/autodiff/GridInit.hpp
View File

@ -54,7 +54,7 @@ namespace Opm
public:
/// Initialize from a deck and/or an eclipse state and (logical cartesian) specified pore volumes.
GridInit(EclipseStateConstPtr eclipse_state, const std::vector<double>& porv)
: grid_manager_(eclipse_state->getInputGrid(), porv)
: grid_manager_(*eclipse_state->getInputGrid(), porv)
{
}
/// Access the created grid.
@ -76,7 +76,7 @@ namespace Opm
/// Initialize from a deck and/or an eclipse state and (logical cartesian) specified pore volumes.
GridInit(EclipseStateConstPtr eclipse_state, const std::vector<double>& porv)
{
grid_.processEclipseFormat(eclipse_state->getInputGrid(), false, false, false, porv);
grid_.processEclipseFormat(*eclipse_state->getInputGrid(), false, false, false, porv);
}
/// Access the created grid. Note that mutable access may be required for load balancing.
Dune::CpGrid& grid()

2
opm/autodiff/MissingFeatures.cpp
View File

@ -77,7 +77,7 @@ namespace MissingFeatures {
"GRUPNET", "IMKRVD", "IMPES", "IMPTVD", "MAPUNITS",
"MAXVALUE", "MESSAGES", "MINVALUE", "MONITOR", "MSGFILE",
"MULT_XYZ", "NETBALAN", "NEXTSTEP", "NOCASC", "NOECHO",
"NOGGF", "NOINSPEC", "NOMONITO", "NONNC", "NORSSPEC", "NOSIM",
"NOGGF", "NOINSPEC", "NOMONITO", "NONNC", "NORSSPEC",
"NSTACK", "NUMRES", "NUPCOL", "OILVISCT", "OLDTRAN", "OPTIONS",
"PARALLEL", "PBVD", "PCG", "PERMXY", "PERMYZ",
"PERMZX", "PIMULTAB", "PLYADSS", "PLYDHFLF",

8
opm/autodiff/NonlinearSolver.hpp
View File

@ -125,6 +125,14 @@ namespace Opm {
/// Number of well iterations used in all calls to step().
int wellIterationsLastStep() const;
/// Compute fluid in place.
/// \param[in] ReservoirState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum) const;
/// Reference to physical model.
const PhysicalModel& model() const;

9
opm/autodiff/NonlinearSolver_impl.hpp
View File

@ -99,6 +99,15 @@ namespace Opm
return wellIterationsLast_;
}
template <class PhysicalModel>
std::vector<V>
NonlinearSolver<PhysicalModel>::computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum) const
{
return model_->computeFluidInPlace(x, fipnum);
}
template <class PhysicalModel>
int
NonlinearSolver<PhysicalModel>::

18
opm/autodiff/ParallelDebugOutput.hpp
View File

@ -45,10 +45,12 @@ namespace Opm
public:
virtual ~ParallelDebugOutputInterface() {}
// gather solution to rank 0 for EclipseWriter
//! \brief gather solution to rank 0 for EclipseWriter
//! \param localWellState The well state
//! \param wellStateStepNumber The step number of the well state.
virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
const WellState& localWellState,
const int reportStep ) = 0;
const int wellStateStepNumber ) = 0;
virtual const SimulationDataContainer& globalReservoirState() const = 0 ;
virtual const WellState& globalWellState() const = 0 ;
@ -77,7 +79,7 @@ namespace Opm
// gather solution to rank 0 for EclipseWriter
virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
const WellState& localWellState,
const int /* reportStep */)
const int /* wellStateStepNumber */)
{
globalState_ = &localReservoirState;
wellState_ = &localWellState;
@ -382,9 +384,11 @@ namespace Opm
void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
{
// write all cell data registered in local state
for (auto& pair : globalState_.cellData()) {
// we loop over the data of the local state as
// its order governs the order the data got received.
for (auto& pair : localState_.cellData()) {
const std::string& key = pair.first;
auto& data = pair.second;
auto& data = globalState_.getCellData(key);
const size_t stride = globalState_.numCellDataComponents( key );
for( size_t i=0; i<stride; ++i )
@ -519,7 +523,7 @@ namespace Opm
// gather solution to rank 0 for EclipseWriter
bool collectToIORank( const SimulationDataContainer& localReservoirState,
const WellState& localWellState,
const int reportStep )
const int wellStateStepNumber )
{
if( isIORank() )
{
@ -530,7 +534,7 @@ namespace Opm
const DynamicListEconLimited dynamic_list_econ_limited;
// Create wells and well state.
WellsManager wells_manager(eclipseState_,
reportStep,
wellStateStepNumber,
Opm::UgGridHelpers::numCells( globalGrid ),
Opm::UgGridHelpers::globalCell( globalGrid ),
Opm::UgGridHelpers::cartDims( globalGrid ),

10
opm/autodiff/SimulatorBase.hpp
View File

@ -160,6 +160,16 @@ namespace Opm
const BlackoilState& x,
WellState& xw);
void
FIPUnitConvert(const UnitSystem& units,
std::vector<V>& fip);
V
FIPTotals(const std::vector<V>& fip, const ReservoirState& state);
void
outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg);
void computeWellPotentials(const Wells* wells,
const WellState& xw,
std::vector<double>& well_potentials);

137
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -111,10 +111,6 @@ namespace Opm
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
}
output_writer_.writeInit( geo_.simProps(grid_) , geo_.nonCartesianConnections( ) );
std::string restorefilename = param_.getDefault("restorefile", std::string("") );
if( ! restorefilename.empty() )
{
@ -135,6 +131,18 @@ namespace Opm
std::vector<double> well_potentials;
DynamicListEconLimited dynamic_list_econ_limited;
bool ooip_computed = false;
std::vector<int> fipnum_global = eclipse_state_->get3DProperties().getIntGridProperty("FIPNUM").getData();
//Get compressed cell fipnum.
std::vector<int> fipnum(AutoDiffGrid::numCells(grid_));
if (fipnum_global.empty()) {
std::fill(fipnum.begin(), fipnum.end(), 0);
} else {
for (size_t c = 0; c < fipnum.size(); ++c) {
fipnum[c] = fipnum_global[AutoDiffGrid::globalCell(grid_)[c]];
}
}
std::vector<V> OOIP;
// Main simulation loop.
while (!timer.done()) {
// Report timestep.
@ -168,7 +176,9 @@ namespace Opm
// write the inital state at the report stage
if (timer.initialStep()) {
output_writer_.writeTimeStep( timer, state, well_state );
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state );
}
// Max oil saturation (for VPPARS), hysteresis update.
@ -183,7 +193,14 @@ namespace Opm
const WellModel well_model(wells);
auto solver = asImpl().createSolver(well_model);
std::unique_ptr<Solver> solver = asImpl().createSolver(well_model);
// Compute orignal FIP;
if (!ooip_computed) {
OOIP = solver->computeFluidInPlace(state, fipnum);
FIPUnitConvert(eclipse_state_->getUnits(), OOIP);
ooip_computed = true;
}
if( terminal_output_ )
{
@ -248,10 +265,21 @@ namespace Opm
// Report timing.
const double st = solver_timer.secsSinceStart();
// Compute current FIP.
std::vector<V> COIP;
COIP = solver->computeFluidInPlace(state, fipnum);
FIPUnitConvert(eclipse_state_->getUnits(), COIP);
V OOIP_totals = FIPTotals(OOIP, state);
V COIP_totals = FIPTotals(COIP, state);
outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
for (size_t reg = 0; reg < OOIP.size(); ++reg) {
outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
}
// accumulate total time
stime += st;
if ( terminal_output_ )
{
std::string msg;
@ -270,7 +298,8 @@ namespace Opm
++timer;
// write simulation state at the report stage
output_writer_.writeTimeStep( timer, state, well_state );
const auto& physicalModel = solver->model();
output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
prev_well_state = well_state;
// The well potentials are only computed if they are needed
@ -626,8 +655,100 @@ namespace Opm
}
template <class Implementation>
void
SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units,
std::vector<V>& fip)
{
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
for (size_t i = 0; i < fip.size(); ++i) {
fip[i][0] = unit::convert::to(fip[i][0], unit::stb);
fip[i][1] = unit::convert::to(fip[i][1], unit::stb);
fip[i][2] = unit::convert::to(fip[i][2], 1000*unit::cubic(unit::feet));
fip[i][3] = unit::convert::to(fip[i][3], 1000*unit::cubic(unit::feet));
fip[i][4] = unit::convert::to(fip[i][4], unit::stb);
fip[i][5] = unit::convert::to(fip[i][5], unit::stb);
fip[i][6] = unit::convert::to(fip[i][6], unit::psia);
}
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
for (size_t i = 0; i < fip.size(); ++i) {
fip[i][6] = unit::convert::to(fip[i][6], unit::barsa);
}
}
}
template <class Implementation>
V
SimulatorBase<Implementation>::FIPTotals(const std::vector<V>& fip, const ReservoirState& state)
{
V totals(V::Zero(7));
for (int i = 0; i < 5; ++i) {
for (size_t reg = 0; reg < fip.size(); ++reg) {
totals[i] += fip[reg][i];
}
}
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = state.numPhases();
const PhaseUsage& pu = props_.phaseUsage();
const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
const V so = pu.phase_used[BlackoilPhases::Liquid] ? V(s.col(BlackoilPhases::Liquid)) : V::Zero(nc);
const V sg = pu.phase_used[BlackoilPhases::Vapour] ? V(s.col(BlackoilPhases::Vapour)) : V::Zero(nc);
const V hydrocarbon = so + sg;
const V p = Eigen::Map<const V>(& state.pressure()[0], nc);
totals[5] = geo_.poreVolume().sum();
totals[6] = unit::convert::to((p * geo_.poreVolume() * hydrocarbon).sum() / ((geo_.poreVolume() * hydrocarbon).sum()), unit::barsa);
return totals;
}
template <class Implementation>
void
SimulatorBase<Implementation>::outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg)
{
std::ostringstream ss;
if (!reg) {
ss << " ===================================================\n"
<< " : Field Totals :\n";
} else {
ss << " ===================================================\n"
<< " : FIPNUM report region "
<< std::setw(2) << reg << " :\n";
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
<< " : Porv volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
}
if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
<< std::fixed << std::setprecision(0)
<< " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
if (!reg) {
ss << " : Pressure is weighted by hydrocarbon pore voulme :\n"
<< " : Pore volumes are taken at reference conditions :\n";
}
ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
}
ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
<< ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
<< std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
<< ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
<< ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
<< std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
<< ":========================:==========================================:================:==========================================:\n";
OpmLog::note(ss.str());
}
template <class Implementation>
void

11
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
View File

@ -72,9 +72,6 @@ namespace Opm
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
}
// init output writer
output_writer_.writeInit( geo_.simProps(grid_) , geo_.nonCartesianConnections( ) );
std::string restorefilename = param_.getDefault("restorefile", std::string("") );
if( ! restorefilename.empty() )
{
@ -128,7 +125,9 @@ namespace Opm
// write the inital state at the report stage
if (timer.initialStep()) {
output_writer_.writeTimeStep( timer, state, well_state );
// No per cell data is written for initial step, but will be
// for subsequent steps, when we have started simulating
output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state );
}
// Max oil saturation (for VPPARS), hysteresis update.
@ -182,8 +181,10 @@ namespace Opm
// Increment timer, remember well state.
++timer;
// write simulation state at the report stage
output_writer_.writeTimeStep( timer, state, well_state );
const auto& physicalModel = solver->model();
output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
prev_well_state = well_state;
}

90
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
View File

@ -42,6 +42,11 @@
#include <boost/filesystem.hpp>
//For OutputWriterHelper
#include <map>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#ifdef HAVE_OPM_GRID
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/version.hh>
@ -52,6 +57,7 @@ namespace Opm
{
void outputStateVtk(const UnstructuredGrid& grid,
const SimulationDataContainer& state,
const int step,
@ -245,14 +251,6 @@ namespace Opm
}
#endif
void
BlackoilOutputWriter::
writeInit(const std::vector<data::CellData>& simProps, const NNC& nnc)
{
if( eclWriter_ ) {
eclWriter_->writeInitAndEgrid(simProps, nnc);
}
}
@ -264,17 +262,20 @@ namespace Opm
std::unique_ptr< SimulatorTimerInterface > timer_;
const SimulationDataContainer state_;
const WellState wellState_;
std::vector<data::CellData> simProps_;
const bool substep_;
explicit WriterCall( BlackoilOutputWriter& writer,
const SimulatorTimerInterface& timer,
const SimulationDataContainer& state,
const WellState& wellState,
const std::vector<data::CellData>& simProps,
bool substep )
: writer_( writer ),
timer_( timer.clone() ),
state_( state ),
wellState_( wellState ),
simProps_( simProps ),
substep_( substep )
{
}
@ -283,18 +284,38 @@ namespace Opm
void run ()
{
// write data
writer_.writeTimeStepSerial( *timer_, state_, wellState_, substep_ );
writer_.writeTimeStepSerial( *timer_, state_, wellState_, simProps_, substep_ );
}
};
}
void
BlackoilOutputWriter::
writeTimeStep(const SimulatorTimerInterface& timer,
writeTimeStepWithoutCellProperties(
const SimulatorTimerInterface& timer,
const SimulationDataContainer& localState,
const WellState& localWellState,
bool substep)
{
std::vector<data::CellData> noCellProperties;
writeTimeStepWithCellProperties(timer, localState, localWellState, noCellProperties, substep);
}
void
BlackoilOutputWriter::
writeTimeStepWithCellProperties(
const SimulatorTimerInterface& timer,
const SimulationDataContainer& localState,
const WellState& localWellState,
const std::vector<data::CellData>& cellData,
bool substep)
{
// VTK output (is parallel if grid is parallel)
if( vtkWriter_ ) {
@ -304,8 +325,15 @@ namespace Opm
bool isIORank = output_ ;
if( parallelOutput_ && parallelOutput_->isParallel() )
{
// If this is not the initial write and no substep, then the well
// state used in the computation is actually the one of the last
// step. We need that well state for the gathering. Otherwise
// It an exception with a message like "global state does not
// contain well ..." might be thrown.
int wellStateStepNumber = ( ! substep && timer.reportStepNum() > 0) ?
(timer.reportStepNum() - 1) : timer.reportStepNum();
// collect all solutions to I/O rank
isIORank = parallelOutput_->collectToIORank( localState, localWellState, timer.reportStepNum() );
isIORank = parallelOutput_->collectToIORank( localState, localWellState, wellStateStepNumber );
}
const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
@ -316,20 +344,23 @@ namespace Opm
{
if( asyncOutput_ ) {
// dispatch the write call to the extra thread
asyncOutput_->dispatch( detail::WriterCall( *this, timer, state, wellState, substep ) );
asyncOutput_->dispatch( detail::WriterCall( *this, timer, state, wellState, cellData, substep ) );
}
else {
// just write the data to disk
writeTimeStepSerial( timer, state, wellState, substep );
writeTimeStepSerial( timer, state, wellState, cellData, substep );
}
}
}
void
BlackoilOutputWriter::
writeTimeStepSerial(const SimulatorTimerInterface& timer,
const SimulationDataContainer& state,
const WellState& wellState,
const std::vector<data::CellData>& simProps,
bool substep)
{
// Matlab output
@ -344,33 +375,7 @@ namespace Opm
if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
} else {
std::vector<data::CellData> simProps;
/*
The simProps vector can be passed to the writeTimestep routine
to add more properties to the restart file. Examples of the
elements for the simProps vector can be the relative
permeabilites KRO, KRG and KRW and the fluxes.
Which properties are requested are configured with the RPTRST
keyword, which is internalized in the RestartConfig class in
EclipseState.
*/
/*
Assuming we already have correctly initialized
std::vector<double> instances kro,krw and krg with the oil,
water and gas relative permeabilities. Then we can write those
to the restart file with:
std::vector<data::CellData> simProps;
simProps.emplace_back( {"KRO" , UnitSystem::measure::identity , kro} );
simProps.emplace_back( {"KRG" , UnitSystem::measure::identity , krg} );
simProps.emplace_back( {"KRW" , UnitSystem::measure::identity , krw} );
*/
eclWriter_->writeTimeStep(timer.currentStepNum(),
eclWriter_->writeTimeStep(timer.reportStepNum(),
substep,
timer.simulationTimeElapsed(),
simToSolution( state, phaseUsage_ ),
@ -441,7 +446,10 @@ namespace Opm
restorefile >> state;
restorefile >> wellState;
writeTimeStep( timer, state, wellState );
// No per cell data is written for restore steps, but will be
// for subsequent steps, when we have started simulating
writeTimeStepWithoutCellProperties( timer, state, wellState );
// some output
std::cout << "Restored step " << timer.reportStepNum() << " at day "
<< unit::convert::to(timer.simulationTimeElapsed(),unit::day) << std::endl;

271
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@ -25,13 +25,13 @@
#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>
#include <opm/core/wells/DynamicListEconLimited.hpp>
#include <opm/output/Cells.hpp>
#include <opm/output/OutputWriter.hpp>
#include <opm/output/eclipse/EclipseWriter.hpp>
#include <opm/autodiff/GridHelpers.hpp>
@ -39,6 +39,7 @@
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/ThreadHandle.hpp>
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
@ -212,23 +213,58 @@ namespace Opm
BlackoilOutputWriter(const Grid& grid,
const parameter::ParameterGroup& param,
Opm::EclipseStateConstPtr eclipseState,
const NNC&,
const Opm::PhaseUsage &phaseUsage,
const double* permeability );
/** \copydoc Opm::OutputWriter::writeInit */
void writeInit(const std::vector<data::CellData>& simProps, const NNC& nnc);
/** \copydoc Opm::OutputWriter::writeTimeStep */
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight. This function will extract the
* requested output cell properties specified by the RPTRST keyword
* and write these to file.
*/
template<class Model>
void writeTimeStep(const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const Opm::WellState& wellState,
const Model& physicalModel,
bool substep = false);
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight. This function will write all
* CellData in simProps to the file as well.
*/
void writeTimeStepWithCellProperties(
const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const Opm::WellState& wellState,
const std::vector<data::CellData>& simProps,
bool substep = false);
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight. This function will not write
* any cell properties (e.g., those requested by RPTRST keyword)
*/
void writeTimeStepWithoutCellProperties(
const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const Opm::WellState& wellState,
bool substep = false);
/** \copydoc Opm::OutputWriter::writeTimeStep */
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight. This is the function which does
* the actual write to file.
*/
void writeTimeStepSerial(const SimulatorTimerInterface& timer,
const SimulationDataContainer& reservoirState,
const Opm::WellState& wellState,
const std::vector<data::CellData>& simProps,
bool substep);
/** \brief return output directory */
@ -285,7 +321,6 @@ namespace Opm
BlackoilOutputWriter(const Grid& grid,
const parameter::ParameterGroup& param,
Opm::EclipseStateConstPtr eclipseState,
const NNC& nnc,
const Opm::PhaseUsage &phaseUsage,
const double* permeability )
: output_( param.getDefault("output", true) ),
@ -301,9 +336,7 @@ namespace Opm
new BlackoilMatlabWriter< Grid >( grid, outputDir_ ) : 0 ),
eclWriter_( output_ && parallelOutput_->isIORank() &&
param.getDefault("output_ecl", true) ?
new EclipseWriter(eclipseState,
parallelOutput_->numCells(),
parallelOutput_->globalCell())
new EclipseWriter(eclipseState,UgGridHelpers::createEclipseGrid( grid , *eclipseState->getInputGrid()))
: 0 ),
eclipseState_(eclipseState),
asyncOutput_()
@ -377,5 +410,227 @@ namespace Opm
}
namespace detail {
/**
* Converts an ADB into a standard vector by copy
*/
inline std::vector<double> adbToDoubleVector(const Opm::AutoDiffBlock<double>& adb) {
const auto& adb_v = adb.value();
std::vector<double> vec(adb_v.data(), adb_v.data() + adb_v.size());
return vec;
}
template<class Model>
std::vector<data::CellData> getCellData(
const Opm::PhaseUsage& phaseUsage,
const Model& model,
const RestartConfig& restartConfig,
const int reportStepNum) {
std::vector<data::CellData> simProps;
//Get the value of each of the keys
std::map<std::string, int> outKeywords = restartConfig.getRestartKeywords(reportStepNum);
for (auto& keyValue : outKeywords) {
keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
}
const typename Model::SimulatorData& sd = model.getSimulatorData();
//Get shorthands for water, oil, gas
const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
const int aqua_idx = phaseUsage.phase_pos[Opm::PhaseUsage::Aqua];
const int liquid_idx = phaseUsage.phase_pos[Opm::PhaseUsage::Liquid];
const int vapour_idx = phaseUsage.phase_pos[Opm::PhaseUsage::Vapour];
/**
* Formation volume factors for water, oil, gas
*/
if (aqua_active && outKeywords["BW"] > 0) {
outKeywords["BW"] = 0;
simProps.emplace_back(data::CellData{
"1OVERBW",
Opm::UnitSystem::measure::water_inverse_formation_volume_factor,
std::move(adbToDoubleVector(sd.rq[aqua_idx].b))});
}
if (liquid_active && outKeywords["BO"] > 0) {
outKeywords["BO"] = 0;
simProps.emplace_back(data::CellData{
"1OVERBO",
Opm::UnitSystem::measure::oil_inverse_formation_volume_factor,
std::move(adbToDoubleVector(sd.rq[liquid_idx].b))});
}
if (vapour_active && outKeywords["BG"] > 0) {
outKeywords["BG"] = 0;
simProps.emplace_back(data::CellData{
"1OVERBG",
Opm::UnitSystem::measure::gas_inverse_formation_volume_factor,
std::move(adbToDoubleVector(sd.rq[vapour_idx].b))});
}
/**
* Densities for water, oil gas
*/
if (outKeywords["DEN"] > 0) {
outKeywords["DEN"] = 0;
if (aqua_active) {
simProps.emplace_back(data::CellData{
"WAT_DEN",
Opm::UnitSystem::measure::density,
std::move(adbToDoubleVector(sd.rq[aqua_idx].rho))});
}
if (liquid_active) {
simProps.emplace_back(data::CellData{
"OIL_DEN",
Opm::UnitSystem::measure::density,
std::move(adbToDoubleVector(sd.rq[liquid_idx].rho))});
}
if (vapour_active) {
simProps.emplace_back(data::CellData{
"GAS_DEN",
Opm::UnitSystem::measure::density,
std::move(adbToDoubleVector(sd.rq[vapour_idx].rho))});
}
}
/**
* Viscosities for water, oil gas
*/
if (outKeywords["VISC"] > 0) {
outKeywords["VISC"] = 0;
if (aqua_active) {
simProps.emplace_back(data::CellData{
"WAT_VISC",
Opm::UnitSystem::measure::viscosity,
std::move(adbToDoubleVector(sd.rq[aqua_idx].mu))});
}
if (liquid_active) {
simProps.emplace_back(data::CellData{
"OIL_VISC",
Opm::UnitSystem::measure::viscosity,
std::move(adbToDoubleVector(sd.rq[liquid_idx].mu))});
}
if (vapour_active) {
simProps.emplace_back(data::CellData{
"GAS_VISC",
Opm::UnitSystem::measure::viscosity,
std::move(adbToDoubleVector(sd.rq[vapour_idx].mu))});
}
}
/**
* Relative permeabilities for water, oil, gas
*/
if (aqua_active && outKeywords["KRW"] > 0) {
if (sd.rq[aqua_idx].kr.size() > 0) {
outKeywords["KRW"] = 0;
simProps.emplace_back(data::CellData{
"WATKR",
Opm::UnitSystem::measure::permeability,
std::move(adbToDoubleVector(sd.rq[aqua_idx].kr))});
}
else {
Opm::OpmLog::warning("Empty:WATKR",
"Not emitting empty Water Rel-Perm");
}
}
if (liquid_active && outKeywords["KRO"] > 0) {
if (sd.rq[liquid_idx].kr.size() > 0) {
outKeywords["KRO"] = 0;
simProps.emplace_back(data::CellData{
"OILKR",
Opm::UnitSystem::measure::permeability,
std::move(adbToDoubleVector(sd.rq[liquid_idx].kr))});
}
else {
Opm::OpmLog::warning("Empty:OILKR",
"Not emitting empty Oil Rel-Perm");
}
}
if (vapour_active && outKeywords["KRG"] > 0) {
if (sd.rq[vapour_idx].kr.size() > 0) {
outKeywords["KRG"] = 0;
simProps.emplace_back(data::CellData{
"GASKR",
Opm::UnitSystem::measure::permeability,
std::move(adbToDoubleVector(sd.rq[vapour_idx].kr))});
}
else {
Opm::OpmLog::warning("Empty:GASKR",
"Not emitting empty Gas Rel-Perm");
}
}
/**
* Vaporized and dissolved gas/oil ratio
*/
if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) {
outKeywords["RSSAT"] = 0;
simProps.emplace_back(data::CellData{
"RSSAT",
Opm::UnitSystem::measure::gas_oil_ratio,
std::move(adbToDoubleVector(sd.rs))});
}
if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) {
outKeywords["RVSAT"] = 0;
simProps.emplace_back(data::CellData{
"RVSAT",
Opm::UnitSystem::measure::oil_gas_ratio,
std::move(adbToDoubleVector(sd.rv))});
}
/**
* Bubble point and dew point pressures
*/
if (vapour_active && liquid_active && outKeywords["PBPD"] > 0) {
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.");
}
//Warn for any unhandled keyword
for (auto& keyValue : outKeywords) {
if (keyValue.second > 0) {
std::string logstring = "Keyword '";
logstring.append(keyValue.first);
logstring.append("' is unhandled for output to file.");
Opm::OpmLog::warning("Unhandled output keyword", logstring);
}
}
return simProps;
}
}
template<class Model>
inline void
BlackoilOutputWriter::
writeTimeStep(const SimulatorTimerInterface& timer,
const SimulationDataContainer& localState,
const WellState& localWellState,
const Model& physicalModel,
bool substep)
{
const RestartConfig& restartConfig = eclipseState_->getRestartConfig();
const int reportStepNum = timer.reportStepNum();
std::vector<data::CellData> cellData = detail::getCellData( phaseUsage_, physicalModel, restartConfig, reportStepNum );
writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep);
}
}
#endif

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();

158
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(),
@ -287,6 +287,9 @@ namespace {
: accum(2, ADB::null())
, mflux( ADB::null())
, b ( ADB::null())
, mu ( ADB::null())
, rho ( ADB::null())
, kr ( ADB::null())
, head ( ADB::null())
, mob ( ADB::null())
, ads (2, ADB::null())
@ -516,17 +519,84 @@ 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];
const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
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;
}
std::vector<V>
FullyImplicitCompressiblePolymerSolver::computeFluidInPlace(const PolymerBlackoilState& x,
const std::vector<int>& fipnum)
{
const int np = x.numPhases();
const int nc = grid_.number_of_cells;
SolutionState state(np);
state.pressure = ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
state.temperature = ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
for (int phase = 0; phase < np; ++phase) {
state.saturation[phase] = ADB::constant(s.col(phase));
}
const ADB& press = state.pressure;
const ADB& temp = state.temperature;
const std::vector<ADB>& sat = state.saturation;
const std::vector<PhasePresence> cond = phaseCondition();
std::vector<ADB> pressure = computePressures(state);
const ADB pv_mult = poroMult(press);
const V& pv = geo_.poreVolume();
std::vector<V> fip(5, V::Zero(nc));
for (int phase = 0; phase < 2; ++phase) {
const ADB& b = fluidReciprocFVF(phase, pressure[phase], temp, cond, cells_);
fip[phase] = (pv_mult * b * sat[phase] * pv).value();
}
const int dims = *std::max_element(fipnum.begin(), fipnum.end());
std::vector<V> values(dims, V::Zero(7));
V hcpv = V::Zero(nc);
V pres = V::Zero(nc);
for (int i = 0; i < 5; ++i) {
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) {
values[fipnum[c]-1][i] += fip[i][c];
}
}
}
// compute PAV and PORV or every regions.
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) {
hcpv[fipnum[c]-1] += pv[c] * s.col(Oil)[c];
pres[fipnum[c]-1] += pv[c] * state.pressure.value()[c];
values[fipnum[c]-1][5] += pv[c];
values[fipnum[c]-1][6] += pv[c] * state.pressure.value()[c] * s.col(Oil)[c];
}
}
for (int reg = 0; reg < dims; ++reg) {
if (hcpv[reg] != 0) {
values[reg][6] /= hcpv[reg];
} else {
values[reg][6] = pres[reg] / values[reg][5];
}
}
return values;
}
@ -568,26 +638,32 @@ 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
// for each active phase.
const V trans = subset(geo_.transmissibility(), ops_.internal_faces);
const std::vector<ADB> kr = computeRelPerm(state);
{
const std::vector<ADB> kr = computeRelPerm(state);
const auto& pu = fluid_.phaseUsage();
for (int phaseIdx=0; phaseIdx < fluid_.numPhases(); ++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, kr[0]);
const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, sd_.rq[0].kr);
const ADB mc = computeMc(state);
computeMassFlux(trans, mc, kr[1], 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 --------
@ -659,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;
@ -883,31 +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());
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_);
sd_.rq[1].mob = tr_mult * kro / mu_o;
sd_.rq[0].mu = mu_w;
sd_.rq[1].mu = mu_o;
}
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;
const ADB mu_o = fluidViscosity(1, press[1], temp, cond, cells_);
rq_[1].mob = tr_mult * kro / mu_o;
for (int phase = 0; phase < 2; ++phase) {
const ADB 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 * 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;
}

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

@ -54,6 +54,39 @@ namespace Opm {
class FullyImplicitCompressiblePolymerSolver
{
public:
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
struct ReservoirResidualQuant {
ReservoirResidualQuant();
std::vector<ADB> accum; // Accumulations
ADB mflux; // Mass flux (surface conditions)
ADB b; // Reciprocal FVF
ADB mu; // Viscosities
ADB rho; // Densities
ADB kr; // Permeabilities
ADB head; // Pressure drop across int. interfaces
ADB mob; // Phase mobility (per cell)
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
/// remain in scope for the lifetime of the solver.
@ -102,24 +135,21 @@ namespace Opm {
double relativeChange(const PolymerBlackoilState& previous,
const PolymerBlackoilState& current ) const;
private:
typedef AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::Array<double,
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor> DataBlock;
/// Return reservoir simulation data (for output functionality)
const SimulatorData& getSimulatorData() const {
return sd_;
}
struct ReservoirResidualQuant {
ReservoirResidualQuant();
std::vector<ADB> accum; // Accumulations
ADB mflux; // Mass flux (surface conditions)
ADB b; // Reciprocal FVF
ADB head; // Pressure drop across int. interfaces
ADB mob; // Phase mobility (per cell)
std::vector<ADB> ads; // Adsorption term.
};
/// Compute fluid in place.
/// \param[in] ReservoirState
/// \param[in] WellState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>
computeFluidInPlace(const PolymerBlackoilState& x,
const std::vector<int>& fipnum);
private:
struct SolutionState {
SolutionState(const int np);
@ -154,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.

4
tests/test_boprops_ad.cpp
View File

@ -71,7 +71,7 @@ struct TestFixture : public Setup
{
TestFixture()
: Setup()
, grid (eclState->getInputGrid())
, grid (*eclState->getInputGrid())
, boprops_ad(deck, eclState, *grid.c_grid(), param.getDefault("init_rock", false))
{
}
@ -89,7 +89,7 @@ struct TestFixtureAd : public Setup
{
TestFixtureAd()
: Setup()
, grid (eclState->getInputGrid())
, grid (*eclState->getInputGrid())
, props(deck, eclState, *grid.c_grid(),
param.getDefault("init_rock", false))
{

2
tests/test_rateconverter.cpp
View File

@ -71,7 +71,7 @@ struct TestFixture : public Setup
{
TestFixture()
: Setup()
, grid (eclState->getInputGrid())
, grid (*eclState->getInputGrid())
, ad_props(deck, eclState, *grid.c_grid(), param.getDefault("init_rock", false))
{
}

16
tests/test_transmissibilitymultipliers.cpp
View File

@ -165,7 +165,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
Opm::DeckConstPtr origDeck = parser->parseString(origDeckString, parseContext);
Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(*origDeck , parseContext));
auto origGridManager = std::make_shared<Opm::GridManager>(origEclipseState->getInputGrid());
auto origGridManager = std::make_shared<Opm::GridManager>(*origEclipseState->getInputGrid());
auto origProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(origDeck, origEclipseState, *(origGridManager->c_grid()));
Opm::DerivedGeology origGeology(*(origGridManager->c_grid()), *origProps, origEclipseState, false);
@ -176,7 +176,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
Opm::DeckConstPtr multDeck = parser->parseString(multDeckString, parseContext);
Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(*multDeck, parseContext));
auto multGridManager = std::make_shared<Opm::GridManager>(multEclipseState->getInputGrid());
auto multGridManager = std::make_shared<Opm::GridManager>(*multEclipseState->getInputGrid());
auto multProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multDeck, multEclipseState, *(multGridManager->c_grid()));
Opm::DerivedGeology multGeology(*(multGridManager->c_grid()), *multProps, multEclipseState, false);
@ -188,7 +188,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString, parseContext);
Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(*multMinusDeck , parseContext));
auto multMinusGridManager = std::make_shared<Opm::GridManager>(multMinusEclipseState->getInputGrid());
auto multMinusGridManager = std::make_shared<Opm::GridManager>(*multMinusEclipseState->getInputGrid());
auto multMinusProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multMinusDeck, multMinusEclipseState, *(multMinusGridManager->c_grid()));
Opm::DerivedGeology multMinusGeology(*(multMinusGridManager->c_grid()), *multMinusProps, multMinusEclipseState, false);
@ -199,7 +199,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseContext);
Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(*ntgDeck, parseContext));
auto ntgGridManager = std::make_shared<Opm::GridManager>(ntgEclipseState->getInputGrid());
auto ntgGridManager = std::make_shared<Opm::GridManager>(*ntgEclipseState->getInputGrid());
auto ntgProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(ntgDeck, ntgEclipseState, *(ntgGridManager->c_grid()));
Opm::DerivedGeology ntgGeology(*(ntgGridManager->c_grid()), *ntgProps, ntgEclipseState, false);
@ -281,7 +281,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(*origDeck , parseContext));
auto origGrid = std::make_shared<Dune::CpGrid>();
origGrid->processEclipseFormat(origEclipseState->getInputGrid(), 0.0, false);
origGrid->processEclipseFormat(*origEclipseState->getInputGrid(), 0.0, false);
auto origProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(origDeck,
origEclipseState,
@ -296,7 +296,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(*multDeck, parseContext));
auto multGrid = std::make_shared<Dune::CpGrid>();
multGrid->processEclipseFormat(multEclipseState->getInputGrid(), 0.0, false);
multGrid->processEclipseFormat(*multEclipseState->getInputGrid(), 0.0, false);
auto multProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multDeck, multEclipseState, *multGrid);
@ -310,7 +310,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(*multMinusDeck, parseContext));
auto multMinusGrid = std::make_shared<Dune::CpGrid>();
multMinusGrid->processEclipseFormat(multMinusEclipseState->getInputGrid(), 0.0, false);
multMinusGrid->processEclipseFormat(*multMinusEclipseState->getInputGrid(), 0.0, false);
auto multMinusProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multMinusDeck, multMinusEclipseState, *multMinusGrid);
@ -324,7 +324,7 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(*ntgDeck, parseContext));
auto ntgGrid = std::make_shared<Dune::CpGrid>();
ntgGrid->processEclipseFormat(ntgEclipseState->getInputGrid(), 0.0, false);
ntgGrid->processEclipseFormat(*ntgEclipseState->getInputGrid(), 0.0, false);
auto ntgProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(ntgDeck, ntgEclipseState, *ntgGrid);