mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge pull request #3968 from totto82/drift_fix
Adapt drift compensation code to the way Flow checks convergence.
This commit is contained in:
commit
3068424382
@ -1225,7 +1225,7 @@ add_test_compare_restarted_simulation(CASENAME spe9
|
||||
ABS_TOL ${abs_tol_restart}
|
||||
REL_TOL ${rel_tol_restart}
|
||||
RESTART_STEP 15
|
||||
TEST_ARGS --sched-restart=false)
|
||||
TEST_ARGS --sched-restart=false --tolerance-mb=1e-7)
|
||||
|
||||
add_test_compare_restarted_simulation(CASENAME ctaquifer_2d_oilwater
|
||||
FILENAME 2D_OW_CTAQUIFER
|
||||
@ -1248,8 +1248,8 @@ add_test_compare_restarted_simulation(CASENAME fetkovich_2d
|
||||
add_test_compare_restarted_simulation(CASENAME numerical_aquifer_3d_1aqu
|
||||
FILENAME 3D_1AQU_3CELLS
|
||||
SIMULATOR flow
|
||||
ABS_TOL ${abs_tol_restart}
|
||||
REL_TOL ${rel_tol_restart}
|
||||
ABS_TOL 0.4
|
||||
REL_TOL 4.0e-3
|
||||
RESTART_STEP 3
|
||||
DIR aquifer-num
|
||||
TEST_ARGS --sched-restart=true --enable-tuning=true)
|
||||
|
@ -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;
|
||||
@ -1836,29 +1840,32 @@ public:
|
||||
|
||||
// if requested, compensate systematic mass loss for cells which were "well
|
||||
// behaved" in the last time step
|
||||
if (enableDriftCompensation_) {
|
||||
// Note that we don't allow for drift compensation if there is
|
||||
// no source terms i.e. no wells and no aquifers.
|
||||
const auto& schedule = this->simulator().vanguard().schedule();
|
||||
int episodeIdx = this->simulator().episodeIndex();
|
||||
const auto& aquifer = this->simulator().vanguard().eclState().aquifer();
|
||||
bool compensateDrift = schedule.numWells(episodeIdx) > 0 && aquifer.active();
|
||||
if (enableDriftCompensation_ && compensateDrift) {
|
||||
const auto& simulator = this->simulator();
|
||||
const auto& model = this->model();
|
||||
|
||||
// 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();
|
||||
|
||||
// we use a lower tolerance for the compensation too
|
||||
// assure the added drift from the last step does not
|
||||
// cause convergence issues on the current step
|
||||
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 +2950,38 @@ private:
|
||||
return dtNext;
|
||||
}
|
||||
|
||||
void computeAndSetEqWeights_() {
|
||||
std::vector<Scalar> sumInvB(numPhases, 0.0);
|
||||
const auto& gridView = this->gridView();
|
||||
ElementContext elemCtx(this->simulator());
|
||||
for(const auto& elem: elements(gridView, Dune::Partitions::interior)) {
|
||||
elemCtx.updatePrimaryStencil(elem);
|
||||
int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
|
||||
const auto& dofFluidState = initialFluidStates_[elemIdx];
|
||||
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
|
||||
if (!FluidSystem::phaseIsActive(phaseIdx))
|
||||
continue;
|
||||
|
||||
sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx);
|
||||
}
|
||||
}
|
||||
|
||||
size_t numDof = this->model().numGridDof();
|
||||
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_;
|
||||
|
Loading…
Reference in New Issue
Block a user