diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake index ebe6339bf..89fc7a48b 100644 --- a/compareECLFiles.cmake +++ b/compareECLFiles.cmake @@ -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) diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh index 615cac4f2..1400f97d3 100644 --- a/ebos/eclproblem.hh +++ b/ebos/eclproblem.hh @@ -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 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 materialLawManager_;