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

Merge pull request #4544 from hnil/faster_beginstep

Browse Source 
faster begintimestep
This commit is contained in:
Bård Skaflestad 2023-03-24 13:52:25 +01:00 committed by GitHub
commit 3ea6f714db
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 170 additions and 95 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

265
ebos/eclproblem.hh
View File

@ -1052,22 +1052,7 @@ public:
// update maximum water saturation and minimum pressure
// used when ROCKCOMP is activated
const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
const bool invalidateFromMinPressure = updateMinPressure_();
// update hysteresis and max oil saturation used in vappars
const bool invalidateFromHyst = updateHysteresis_();
const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
// the derivatives may have change
bool invalidateIntensiveQuantities = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
if (invalidateIntensiveQuantities){
OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
}
if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
updateMaxPolymerAdsorption_();
asImp_().updateExplicitQuantitites_();
wellModel_.beginTimeStep();
if (enableAquifers_)
@ -1126,9 +1111,10 @@ public:
tracerModel_.endTimeStep();
// deal with DRSDT and DRVDT
updateCompositionChangeLimits_();
asImp_().updateCompositionChangeLimits_();
{
OPM_TIMEBLOCK(driftCompansation);
asImp_().updateCompositionChangeLimits_();
if (enableDriftCompensation_) {
const auto& residual = this->model().linearizer().residual();
for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
@ -1468,7 +1454,7 @@ public:
*/
std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
{ return materialLawManager_; }
template <class FluidState>
void updateRelperms(
std::array<Evaluation,numPhases> &mobility,
@ -1749,6 +1735,7 @@ public:
*/
void initialSolutionApplied()
{
this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/0);
// initialize the wells. Note that this needs to be done after initializing the
// intrinsic permeabilities and the after applying the initial solution because
// the well model uses these...
@ -2081,24 +2068,49 @@ public:
serializer(tracerModel_);
serializer(*materialLawManager_);
serializer(*eclWriter_);
}
private:
Implementation& asImp_()
{ return *static_cast<Implementation *>(this); }
void updateExplicitQuantitites_()
{
OPM_TIMEBLOCK(updateExplicitQuantities);
const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_();
const bool invalidateFromMinPressure = updateMinPressure_();
// update hysteresis and max oil saturation used in vappars
const bool invalidateFromHyst = updateHysteresis_();
const bool invalidateFromMaxOilSat = updateMaxOilSaturation_();
// the derivatives may have change
bool invalidateIntensiveQuantities
= invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat;
if (invalidateIntensiveQuantities) {
OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities);
this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
}
if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>())
updateMaxPolymerAdsorption_();
}
private:
template<class UpdateFunc>
void updateProperty_(const std::string& failureMsg,
UpdateFunc func)
{
OPM_TIMEBLOCK(updateProperty);
ElementContext elemCtx(this->simulator());
const auto& model = this->simulator().model();
const auto& primaryVars = model.solution(/*timeIdx*/0);
const auto& vanguard = this->simulator().vanguard();
size_t numGridDof = primaryVars.size();
OPM_BEGIN_PARALLEL_TRY_CATCH();
for (const auto& elem : elements(vanguard.gridView())) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
unsigned compressedDofIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& iq = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
func(compressedDofIdx, iq);
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++dofIdx) {
const auto& iq = *model.cachedIntensiveQuantities(dofIdx, /*timeIdx=*/ 0);
func(dofIdx, iq);
}
OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm());
}
@ -2119,60 +2131,71 @@ private:
this->updateProperty_("EclProblem::updateCompositionChangeLimits_()) failed:",
[this,episodeIdx,active](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
auto& simulator = this->simulator();
auto& vanguard = simulator.vanguard();
if (active[0]) {
// This implements the convective DRSDT as described in
// Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow simulator"
// Submitted to TCCS 11, 2021
const Scalar g = this->gravity_[dim - 1];
const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
const Scalar distZ = vanguard.cellThickness(compressedDofIdx);
const auto& fs = iq.fluidState();
const Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx));
const Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
const Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx));
const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(),t,p);
const Scalar saturatedInvB = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p);
const Scalar rsZero = 0.0;
const Scalar pureDensity = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(),t,p,rsZero) * FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
const Scalar saturatedDensity = saturatedInvB * (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex()) + rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
const Scalar deltaDensity = saturatedDensity - pureDensity;
const Scalar rs = getValue(fs.Rs());
const Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(),t,p,rs);
const Scalar poro = getValue(iq.porosity());
// Note that for so = 0 this gives no limits (inf) for the dissolution rate
// Also we restrict the effect of convective mixing to positive density differences
// i.e. we only allow for fingers moving downward
this->convectiveDrs_[compressedDofIdx] = permz * rssat * max(0.0, deltaDensity) * g / ( so * visc * distZ * poro);
}
this->updateCompositionChangeLimits_(compressedDofIdx,
iq,
episodeIdx,
active);
}
);
}
if (active[1]) {
const auto& fs = iq.fluidState();
void updateCompositionChangeLimits_(unsigned compressedDofIdx, const IntensiveQuantities& iq,int episodeIdx, const std::array<bool,3>& active)
{
auto& simulator = this->simulator();
auto& vanguard = simulator.vanguard();
if (active[0]) {
// This implements the convective DRSDT as described in
// Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow
// simulator" Submitted to TCCS 11, 2021
const Scalar g = this->gravity_[dim - 1];
const DimMatrix& perm = intrinsicPermeability(compressedDofIdx);
const Scalar permz = perm[dim - 1][dim - 1]; // The Z permeability
const Scalar distZ = vanguard.cellThickness(compressedDofIdx);
const auto& fs = iq.fluidState();
const Scalar t = getValue(fs.temperature(FluidSystem::oilPhaseIdx));
const Scalar p = getValue(fs.pressure(FluidSystem::oilPhaseIdx));
const Scalar so = getValue(fs.saturation(FluidSystem::oilPhaseIdx));
const Scalar rssat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), t, p);
const Scalar saturatedInvB
= FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p);
const Scalar rsZero = 0.0;
const Scalar pureDensity
= FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), t, p, rsZero)
* FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex());
const Scalar saturatedDensity = saturatedInvB
* (FluidSystem::oilPvt().oilReferenceDensity(fs.pvtRegionIndex())
+ rssat * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, fs.pvtRegionIndex()));
const Scalar deltaDensity = saturatedDensity - pureDensity;
const Scalar rs = getValue(fs.Rs());
const Scalar visc = FluidSystem::oilPvt().viscosity(fs.pvtRegionIndex(), t, p, rs);
const Scalar poro = getValue(iq.porosity());
// Note that for so = 0 this gives no limits (inf) for the dissolution rate
// Also we restrict the effect of convective mixing to positive density differences
// i.e. we only allow for fingers moving downward
this->convectiveDrs_[compressedDofIdx]
= permz * rssat * max(0.0, deltaDensity) * g / (so * visc * distZ * poro);
}
using FluidState = typename std::decay<decltype(fs)>::type;
if (active[1]) {
const auto& fs = iq.fluidState();
int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
const auto& oilVaporizationControl = vanguard.schedule()[episodeIdx].oilvap();
if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
this->lastRs_[compressedDofIdx] =
BlackOil::template getRs_<FluidSystem,
FluidState,
Scalar>(fs, iq.pvtRegionIndex());
else
this->lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
}
using FluidState = typename std::decay<decltype(fs)>::type;
if (active[2]) {
const auto& fs = iq.fluidState();
using FluidState = typename std::decay<decltype(fs)>::type;
this->lastRv_[compressedDofIdx] =
BlackOil::template getRv_<FluidSystem,
FluidState,
Scalar>(fs, iq.pvtRegionIndex());
}
});
int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx);
const auto& oilVaporizationControl = vanguard.schedule()[episodeIdx].oilvap();
if (oilVaporizationControl.getOption(pvtRegionIdx) || fs.saturation(gasPhaseIdx) > freeGasMinSaturation_)
this->lastRs_[compressedDofIdx]
= BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
else
this->lastRs_[compressedDofIdx] = std::numeric_limits<Scalar>::infinity();
}
if (active[2]) {
const auto& fs = iq.fluidState();
using FluidState = typename std::decay<decltype(fs)>::type;
this->lastRv_[compressedDofIdx]
= BlackOil::template getRv_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
}
}
bool updateMaxOilSaturation_()
@ -2183,19 +2206,30 @@ private:
// we use VAPPARS
if (this->vapparsActive(episodeIdx)) {
this->updateProperty_("EclProblem::updateMaxOilSaturation_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
const auto& fs = iq.fluidState();
const Scalar So = decay<Scalar>(fs.saturation(oilPhaseIdx));
auto& mos = this->maxOilSaturation_;
mos[compressedDofIdx] = std::max(mos[compressedDofIdx], So);
});
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
this->updateMaxOilSaturation_(compressedDofIdx,iq);
});
return true;
}
return false;
}
bool updateMaxOilSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
OPM_TIMEBLOCK_LOCAL(updateMaxOilSaturation);
const auto& fs = iq.fluidState();
const Scalar So = decay<Scalar>(fs.saturation(oilPhaseIdx));
auto& mos = this->maxOilSaturation_;
if(mos[compressedDofIdx] < So){
mos[compressedDofIdx] = So;
return true;
}else{
return false;
}
}
bool updateMaxWaterSaturation_()
{
OPM_TIMEBLOCK(updateMaxWaterSaturation);
@ -2207,14 +2241,26 @@ private:
this->updateProperty_("EclProblem::updateMaxWaterSaturation_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
const auto& fs = iq.fluidState();
const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
auto& mow = this->maxWaterSaturation_;
mow[compressedDofIdx] = std::max(mow[compressedDofIdx], Sw);
this->updateMaxWaterSaturation_(compressedDofIdx,iq);
});
return true;
}
bool updateMaxWaterSaturation_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
OPM_TIMEBLOCK_LOCAL(updateMaxWaterSaturation);
const auto& fs = iq.fluidState();
const Scalar Sw = decay<Scalar>(fs.saturation(waterPhaseIdx));
auto& mow = this->maxWaterSaturation_;
if(mow[compressedDofIdx]< Sw){
mow[compressedDofIdx] = Sw;
return true;
}else{
return false;
}
}
bool updateMinPressure_()
{
OPM_TIMEBLOCK(updateMinPressure);
@ -2225,14 +2271,24 @@ private:
this->updateProperty_("EclProblem::updateMinPressure_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
const auto& fs = iq.fluidState();
const Scalar mo = getValue(fs.pressure(oilPhaseIdx));
auto& mos = this->minOilPressure_;
mos[compressedDofIdx] = std::min(mos[compressedDofIdx], mo);
this->updateMinPressure_(compressedDofIdx,iq);
});
return true;
}
bool updateMinPressure_(unsigned compressedDofIdx, const IntensiveQuantities& iq){
OPM_TIMEBLOCK_LOCAL(updateMinPressure);
const auto& fs = iq.fluidState();
const Scalar mo = getValue(fs.pressure(oilPhaseIdx));
auto& mos = this->minOilPressure_;
if(mos[compressedDofIdx]> mo){
mos[compressedDofIdx] = mo;
return true;
}else{
return false;
}
}
void readMaterialParameters_()
{
OPM_TIMEBLOCK(readMaterialParameters);
@ -2718,18 +2774,37 @@ private:
return true;
}
bool updateHysteresis_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
OPM_TIMEBLOCK_LOCAL(updateHysteresis_);
materialLawManager_->updateHysteresis(iq.fluidState(), compressedDofIdx);
//TODO change materials to give a bool
return true;
}
void updateMaxPolymerAdsorption_()
{
// we need to update the max polymer adsoption data for all elements
this->updateProperty_("EclProblem::updateMaxPolymerAdsorption_() failed:",
[this](unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
const Scalar pa = scalarValue(iq.polymerAdsorption());
auto& mpa = this->maxPolymerAdsorption_;
mpa[compressedDofIdx] = std::max(mpa[compressedDofIdx], pa);
this->updateMaxPolymerAdsorption_(compressedDofIdx,iq);
});
}
bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq)
{
const Scalar pa = scalarValue(iq.polymerAdsorption());
auto& mpa = this->maxPolymerAdsorption_;
if(mpa[compressedDofIdx]<pa){
mpa[compressedDofIdx] = pa;
return true;
}else{
return false;
}
}
struct PffDofData_
{
ConditionalStorage<enableEnergy, Scalar> thermalHalfTransIn;