From 8250a815ccd4cc4af527ef183f354292c6e7336e Mon Sep 17 00:00:00 2001 From: hnil Date: Thu, 23 Mar 2023 10:12:02 +0100 Subject: [PATCH] - faster updateProperty - refactored for making local updating in inhereted classes --- ebos/eclproblem.hh | 263 +++++++++++++++++++++++++++++---------------- 1 file changed, 169 insertions(+), 94 deletions(-) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 04a70ba89..8ef6eddc3 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -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()) - 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 ++) { @@ -2081,24 +2067,50 @@ public: serializer(tracerModel_); serializer(*materialLawManager_); serializer(*eclWriter_); + } +private: + Implementation& asImp_() + { return *static_cast(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()) + updateMaxPolymerAdsorption_(); } -private: template void updateProperty_(const std::string& failureMsg, UpdateFunc func) { OPM_TIMEBLOCK(updateProperty); - ElementContext elemCtx(this->simulator()); + const auto& problem = this->simulator().problem(); + 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& 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::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_(fs, iq.pvtRegionIndex()); - else - this->lastRs_[compressedDofIdx] = std::numeric_limits::infinity(); - } + using FluidState = typename std::decay::type; - if (active[2]) { - const auto& fs = iq.fluidState(); - using FluidState = typename std::decay::type; - this->lastRv_[compressedDofIdx] = - BlackOil::template getRv_(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_(fs, iq.pvtRegionIndex()); + else + this->lastRs_[compressedDofIdx] = std::numeric_limits::infinity(); + } + + if (active[2]) { + const auto& fs = iq.fluidState(); + using FluidState = typename std::decay::type; + this->lastRv_[compressedDofIdx] + = BlackOil::template getRv_(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(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(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(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(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] thermalHalfTransIn;