fix the issues that emerged in the context of [at]totto82's review

This commit is contained in:
Andreas Lauser 2018-01-30 16:15:47 +01:00
parent 54c96aa1c2
commit 26228ec5f3
3 changed files with 40 additions and 33 deletions

View File

@ -77,12 +77,15 @@ class EclEquilInitializer
enum { waterCompIdx = FluidSystem::waterCompIdx };
enum { dimWorld = GridView::dimensionworld };
enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
// NB: setting the enableEnergy argument to true enables storage of enthalpy and
// internal energy!
typedef Opm::BlackOilFluidState<Scalar,
FluidSystem,
/*enableTemperature=*/true,
/*enableEnthalpy=*/enableEnergy> ScalarFluidState;
enableTemperature,
enableEnergy> ScalarFluidState;
public:
template <class EclMaterialLawManager>
@ -102,9 +105,6 @@ public:
vanguard.grid(),
simulator.problem().gravity()[dimWorld - 1]);
const std::vector<double>& tempiData =
eclState.get3DProperties().getDoubleGridProperty("TEMPI").getData();
// copy the result into the array of initial fluid states
initialFluidStates_.resize(numCartesianElems);
for (unsigned int elemIdx = 0; elemIdx < numElems; ++elemIdx) {
@ -135,18 +135,12 @@ public:
// set the temperature.
//
// TODO: setting temperature explicitly while computing static equilibirum
// for everything else is a bit inconsistent, i.e.,
// Opm::initStateEquil() should be extended to provide correct initial
// temperatures
assert(std::isfinite(tempiData[cartesianElemIdx]));
fluidState.setTemperature(tempiData[cartesianElemIdx]);
if (enableTemperature || enableEnergy)
fluidState.setTemperature(initialState.temperature()[elemIdx]);
// set the phase pressures.
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fluidState.setPressure(phaseIdx, initialState.press()[phaseIdx][elemIdx]);
}
}

View File

@ -268,6 +268,15 @@ SET_BOOL_PROP(EclBaseProblem, DisableWells, false);
// By default, we enable the debugging checks if we're compiled in debug mode
SET_BOOL_PROP(EclBaseProblem, EnableDebuggingChecks, true);
// store temperature (but do not conserve energy, as long as EnableEnergy is false)
SET_BOOL_PROP(EclBaseProblem, EnableTemperature, true);
// disable all extensions supported by black oil model. this should not really be
// necessary but it makes things a bit more explicit
SET_BOOL_PROP(EclBaseProblem, EnablePolymer, false);
SET_BOOL_PROP(EclBaseProblem, EnableSolvent, false);
SET_BOOL_PROP(EclBaseProblem, EnableEnergy, false);
} // namespace Properties
/*!
@ -296,6 +305,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
enum { numComponents = FluidSystem::numComponents };
enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
enum { enableTemperature = GET_PROP_VALUE(TypeTag, EnableTemperature) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
@ -325,8 +335,8 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
typedef Opm::BlackOilFluidState<Scalar,
FluidSystem,
/*enableTemperature=*/true,
/*enableEnthalpy=*/enableEnergy> InitialFluidState;
enableTemperature,
enableEnergy> InitialFluidState;
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;

View File

@ -649,18 +649,6 @@ phasePressures(const Grid& grid,
return press;
}
template <class Grid,
class Region,
class CellRange>
std::vector<double>
temperature(const Grid& /* G */,
const Region& /* reg */,
const CellRange& cells)
{
// use the standard temperature for everything for now
return std::vector<double>(cells.size(), 273.15 + 20.0);
}
/**
* Compute initial phase saturations by means of equilibration.
*
@ -945,7 +933,8 @@ public:
const Grid& grid,
const double grav = Opm::unit::gravity,
const bool applySwatInit = true)
: pp_(FluidSystem::numPhases,
: temperature_(grid.size(/*codim=*/0)),
pp_(FluidSystem::numPhases,
std::vector<double>(grid.size(/*codim=*/0))),
sat_(FluidSystem::numPhases,
std::vector<double>(grid.size(/*codim=*/0))),
@ -1069,8 +1058,11 @@ public:
}
}
// extract the initial temperature
updateInitialTemperature_(eclipseState);
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, materialLawManager, grid, grav);
calcPressSatRsRv(eclipseState, eqlmap, rec, materialLawManager, grid, grav);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
@ -1079,16 +1071,27 @@ public:
typedef std::vector<double> Vec;
typedef std::vector<Vec> PVec; // One per phase.
const Vec& temperature() const { return temperature_; }
const PVec& press() const { return pp_; }
const PVec& saturation() const { return sat_; }
const Vec& rs() const { return rs_; }
const Vec& rv() const { return rv_; }
private:
void updateInitialTemperature_(const Opm::EclipseState& eclState)
{
// Get the initial temperature data
const std::vector<double>& tempiData =
eclState.get3DProperties().getDoubleGridProperty("TEMPI").getData();
temperature_ = tempiData;
}
typedef EquilReg EqReg;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rsFunc_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rvFunc_;
std::vector<int> regionPvtIdx_;
Vec temperature_;
PVec pp_;
PVec sat_;
Vec rs_;
@ -1123,7 +1126,8 @@ private:
}
template <class RMap, class MaterialLawManager>
void calcPressSatRsRv(const RMap& reg,
void calcPressSatRsRv(const Opm::EclipseState& eclState,
const RMap& reg,
const std::vector< Opm::EquilRecord >& rec,
MaterialLawManager& materialLawManager,
const Grid& grid,
@ -1140,7 +1144,6 @@ private:
const EqReg eqreg(rec[r], rsFunc_[r], rvFunc_[r], regionPvtIdx_[r]);
PVec pressures = phasePressures<FluidSystem>(grid, eqreg, cells, grav);
const std::vector<double>& temp = temperature(grid, eqreg, cells);
const PVec sat = phaseSaturations<FluidSystem>(grid, eqreg, cells, materialLawManager, swatInit_, pressures);
const int np = FluidSystem::numPhases;
@ -1153,8 +1156,8 @@ private:
if (oil && gas) {
const int oilpos = FluidSystem::oilPhaseIdx;
const int gaspos = FluidSystem::gasPhaseIdx;
const Vec rsVals = computeRs(grid, cells, pressures[oilpos], temp, *(rsFunc_[r]), sat[gaspos]);
const Vec rvVals = computeRs(grid, cells, pressures[gaspos], temp, *(rvFunc_[r]), sat[oilpos]);
const Vec rsVals = computeRs(grid, cells, pressures[oilpos], temperature_, *(rsFunc_[r]), sat[gaspos]);
const Vec rvVals = computeRs(grid, cells, pressures[gaspos], temperature_, *(rvFunc_[r]), sat[oilpos]);
copyFromRegion(rsVals, cells, rs_);
copyFromRegion(rvVals, cells, rv_);
}