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

Restrict drift compensation by CNV tolerance

Browse Source
This commit is contained in:
Tor Harald Sandve
2022-05-05 13:43:16 +02:00
parent 1ca731a388
commit 1f54e90f33
1 changed files with 39 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
ebos/eclproblem.hh
View File

@@ -884,6 +884,7 @@ public:
readEclRestartSolution_();
else
readInitialCondition_();
tracerModel_.prepareTracerBatches();
updatePffDofData_();
@@ -897,6 +898,9 @@ public:
readBoundaryConditions_();
// compute and set eq weights based on initial b values
computeAndSetEqWeights_();
if (enableDriftCompensation_) {
drift_.resize(this->model().numGridDof());
drift_ = 0.0;
@@ -1842,23 +1846,19 @@ public:
// we need a higher maxCompensation than the Newton tolerance because the
// current time step might be shorter than the last one
Scalar maxCompensation = 10.0*model.newtonMethod().tolerance();
Scalar maxCompensation = model.newtonMethod().tolerance()/10;
Scalar poro = this->porosity(globalDofIdx, timeIdx);
Scalar dt = simulator.timeStepSize();
EqVector dofDriftRate = drift_[globalDofIdx];
dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx);
// compute the weighted total drift rate
Scalar totalDriftRate = 0.0;
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
totalDriftRate +=
std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
// make sure that we do not exceed the maximum rate of drift compensation
if (totalDriftRate > maxCompensation)
dofDriftRate *= maxCompensation/totalDriftRate;
// restrict drift compensation to the CNV tolerance
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro;
if (cnv > maxCompensation) {
dofDriftRate[eqIdx] *= maxCompensation/cnv;
}
}
for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
rate[eqIdx] -= dofDriftRate[eqIdx];
@@ -2943,6 +2943,33 @@ private:
return dtNext;
}
void computeAndSetEqWeights_() {
std::vector<Scalar> sumInvB(numPhases, 0.0);
size_t numDof = this->model().numGridDof();
for (size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
const auto& dofFluidState = initialFluidStates_[dofIdx];
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
}
}
const auto& comm = this->simulator().vanguard().grid().comm();
comm.sum(sumInvB.data(),sumInvB.size());
Scalar numTotalDof = comm.sum(numDof);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
Scalar avgB = numTotalDof / sumInvB[phaseIdx];
unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
this->model().setEqWeight(activeSolventCompIdx, avgB);
}
}
typename Vanguard::TransmissibilityType transmissibilities_;
std::shared_ptr<EclMaterialLawManager> materialLawManager_;