Extended the support for keywords for restart file output

This commit is contained in:
babrodtk 2016-09-02 09:38:25 +02:00
parent 35bed24465
commit 40b2b95d87
6 changed files with 221 additions and 85 deletions

View File

@ -117,6 +117,18 @@ 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)
};
typedef typename ModelTraits<Implementation>::ReservoirState ReservoirState;
typedef typename ModelTraits<Implementation>::WellState WellState;
typedef typename ModelTraits<Implementation>::ModelParameters ModelParameters;
@ -268,15 +280,9 @@ namespace Opm {
computeFluidInPlace(const ReservoirState& x,
const std::vector<int>& fipnum);
const ADB& getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex phase) const {
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
if (pu.phase_used[phase]) {
const int pos = pu.phase_pos[phase];
return rq_[pos].b;
}
else {
return ADB::null();
}
/// Return reservoir residual quantitites (in particular for output functionality)
const std::vector<ReservoirResidualQuant>& getReservoirResidualQuantities() const {
return rq_;
}
protected:
@ -288,14 +294,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 ---------

View File

@ -457,6 +457,9 @@ namespace detail {
: 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())
{
@ -830,13 +833,19 @@ namespace detail {
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) {
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);
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);
residual_.material_balance_eq[ phaseIdx ] =
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])

View File

@ -53,6 +53,9 @@ namespace Opm {
typedef BlackoilSequentialModelParameters ModelParameters;
typedef DefaultBlackoilSolutionState SolutionState;
typedef BlackoilTransportModel<Grid, WellModel> TransportModel;
typedef typename TransportModel::ReservoirResidualQuant ReservoirResidualQuant;
/// 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.
@ -260,14 +263,14 @@ namespace Opm {
}
const ADB& getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex phase) const {
return transport_model_->getReciprocalFormationVolumeFactor(phase);
/// Return reservoir residual quantitites (in particular for output functionality)
const std::vector<ReservoirResidualQuant>& getReservoirResidualQuantities() const {
return transport_model_->getReservoirResidualQuantities();
}
protected:
typedef BlackoilPressureModel<Grid, WellModel> PressureModel;
typedef BlackoilTransportModel<Grid, WellModel> TransportModel;
typedef NonlinearSolver<PressureModel> PressureSolver;
typedef NonlinearSolver<TransportModel> TransportSolver;

View File

@ -38,6 +38,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>
@ -412,35 +413,48 @@ namespace Opm
};
/**
* Converts an ADB into a standard vector by copy
*/
inline std::vector<double> adbToDoubleVector(const Opm::AutoDiffBlock<double>& input) {
const auto& b_v = input.value();
std::vector<double> b(b_v.data(), b_v.data() + b_v.size());
return b;
}
template<class Model>
std::vector<data::CellData> getCellData(const Model& model,
std::vector<data::CellData> getCellData(
const Opm::PhaseUsage& phaseUsage,
const Model& model,
const RestartConfig& restartConfig,
const int reportStepNum) {
std::vector<data::CellData> simProps;
std::vector<data::CellData> simProps;
std::map<const char*, int> outKeywords {
{"ALLPROPS", 0},
// Formation volume factors
{"BG", 0},
{"BO", 0},
{"BW", 0},
{"CONV", 0},
{"DEN", 0},
{"CONV", 0}, // < Cells with convergence problems
{"DEN", 0}, // < Densities
// Relative permeabilities
{"KRG", 0},
{"KRO", 0},
{"KRW", 0},
{"RVSAT", 0},
{"RSSAT", 0},
{"RVSAT", 0}, // < Vaporized gas/oil ratio
{"RSSAT", 0}, // < Dissolved gas/oil ratio
{"NORST", 0},
{"PBPD", 0},
{"VISC", 0}
{"NORST", 0}, // < Visualization restart file only
{"PBPD", 0}, // < Bubble point and dew point pressures
{"VISC", 0} // < Viscosities
};
//Get the value of each of the keys
@ -455,44 +469,143 @@ namespace Opm
outKeywords["BO"] = std::max(outKeywords["BO"], 1);
outKeywords["BW"] = std::max(outKeywords["BW"], 1);
outKeywords["DEN"] = std::max(outKeywords["DEN"], 1);
outKeywords["KRG"] = std::max(outKeywords["KRG"], 1);
outKeywords["KRO"] = std::max(outKeywords["KRO"], 1);
outKeywords["KRW"] = std::max(outKeywords["KRW"], 1);
outKeywords["DEN"] = std::max(outKeywords["DEN"], 1);
outKeywords["VISC"] = std::max(outKeywords["VISC"], 1);
}
const std::vector<typename Model::ReservoirResidualQuant>& rq = model.getReservoirResidualQuantities();
if (outKeywords["BW"] > 0) {
const auto& b_adb = model.getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex::Aqua);
const auto& b_v = b_adb.value();
const std::vector<double> b(b_v.data(), b_v.data() + b_v.size());
simProps.emplace_back("1OVERBW", Opm::UnitSystem::measure::volume, b);
//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) {
simProps.emplace_back(
"1OVERBW",
Opm::UnitSystem::measure::volume,
adbToDoubleVector(rq[aqua_idx].b));
}
if (outKeywords["BO"] > 0) {
const auto& b_adb = model.getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex::Liquid);
const auto& b_v = b_adb.value();
const std::vector<double> b(b_v.data(), b_v.data() + b_v.size());
simProps.emplace_back("1OVERBO", Opm::UnitSystem::measure::volume, b);
if (liquid_active && outKeywords["BO"] > 0) {
simProps.emplace_back(
"1OVERBO",
Opm::UnitSystem::measure::volume,
adbToDoubleVector(rq[liquid_idx].b));
}
if (outKeywords["BG"] > 0) {
const auto& b_adb = model.getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex::Vapour);
const auto& b_v = b_adb.value();
const std::vector<double> b(b_v.data(), b_v.data() + b_v.size());
simProps.emplace_back("1OVERBG", Opm::UnitSystem::measure::volume, b);
if (vapour_active && outKeywords["BG"] > 0) {
simProps.emplace_back(
"1OVERBG",
Opm::UnitSystem::measure::volume,
adbToDoubleVector(rq[vapour_idx].b));
}
/**
* Densities for water, oil gas
*/
if (outKeywords["DEN"] > 0) {
if (aqua_active) {
simProps.emplace_back(
"WAT_DEN",
Opm::UnitSystem::measure::density,
adbToDoubleVector(rq[aqua_idx].rho));
}
if (liquid_active) {
simProps.emplace_back(
"OIL_DEN",
Opm::UnitSystem::measure::density,
adbToDoubleVector(rq[liquid_idx].rho));
}
if (vapour_active) {
simProps.emplace_back(
"GAS_DEN",
Opm::UnitSystem::measure::density,
adbToDoubleVector(rq[vapour_idx].rho));
}
}
/**
* Viscosities for water, oil gas
*/
if (outKeywords["VISC"] > 0) {
if (aqua_active) {
simProps.emplace_back(
"WAT_VISC",
Opm::UnitSystem::measure::viscosity,
adbToDoubleVector(rq[aqua_idx].mu));
}
if (liquid_active) {
simProps.emplace_back(
"OIL_VISC",
Opm::UnitSystem::measure::viscosity,
adbToDoubleVector(rq[liquid_idx].mu));
}
if (vapour_active) {
simProps.emplace_back(
"GAS_VISC",
Opm::UnitSystem::measure::viscosity,
adbToDoubleVector(rq[vapour_idx].mu));
}
}
/**
* Relative permeabilities for water, oil, gas
*/
if (aqua_active && outKeywords["KRW"] > 0) {
simProps.emplace_back(
"WATKR",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(rq[aqua_idx].kr));
}
if (aqua_active && outKeywords["KRO"] > 0) {
simProps.emplace_back(
"OILKR",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(rq[liquid_idx].kr));
}
if (aqua_active && outKeywords["KRG"] > 0) {
simProps.emplace_back(
"GASKR",
Opm::UnitSystem::measure::permeability,
adbToDoubleVector(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?
}
if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) {
}
return simProps;
}
/**
* Template specialization to print raw cell data
* Template specialization to print raw cell data. That is, if the
* model argument is a vector of celldata, simply return that as-is.
*/
template<>
inline
std::vector<data::CellData> getCellData<std::vector<data::CellData> >(const std::vector<data::CellData>& model,
std::vector<data::CellData> getCellData<std::vector<data::CellData> >(
const Opm::PhaseUsage& phaseUsage,
const std::vector<data::CellData>& model,
const RestartConfig& restartConfig,
const int reportStepNum) {
return model;
@ -535,7 +648,7 @@ namespace Opm
const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
const RestartConfig& restartConfig = eclipseState_->getRestartConfig();
const int reportStepNum = timer.reportStepNum();
std::vector<data::CellData> cellData = detail::getCellData( physicalModel, restartConfig, reportStepNum );
std::vector<data::CellData> cellData = detail::getCellData( phaseUsage_, physicalModel, restartConfig, reportStepNum );
// serial output is only done on I/O rank
if( isIORank )

View File

@ -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())
@ -644,11 +647,17 @@ namespace {
// 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) {
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, rq_[0].kr);
const ADB mc = computeMc(state);
computeMassFlux(trans, mc, kr[1], krw_eff, 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])
@ -953,17 +962,23 @@ namespace {
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;
const ADB mu_o = fluidViscosity(1, press[1], temp, cond, cells_);
rq_[1].mob = tr_mult * kro / 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;
rq_[0].mu = mu_w;
rq_[1].mu = mu_o;
}
for (int phase = 0; phase < 2; ++phase) {
const ADB rho = fluidDensity(phase, press[phase], temp, cond, cells_);
ADB& head = rq_[ phase ].head;
rq_[phase].rho = fluidDensity(phase, press[phase], temp, cond, cells_);
ADB& head = 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 * 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());

View File

@ -62,6 +62,20 @@ 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 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.
};
/// 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.
@ -110,6 +124,11 @@ 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_;
}
/// Compute fluid in place.
/// \param[in] ReservoirState
/// \param[in] WellState
@ -121,27 +140,6 @@ namespace Opm {
private:
const ADB& getReciprocalFormationVolumeFactor(PhaseUsage::PhaseIndex phase) const {
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
if (pu.phase_used[phase]) {
const int pos = pu.phase_pos[phase];
return rq_[pos].b;
}
else {
return ADB::null();
}
}
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.
};
struct SolutionState {
SolutionState(const int np);
ADB pressure;