OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-12 03:05:38 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5330 from bska/run-udq-assign-after-action

Browse Source 
Rerun UDQ Assignments After Action Processing
This commit is contained in:
Markus Blatt 2024-05-21 17:19:54 +02:00 committed by GitHub
commit 15001cbf7e
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
1 changed files with 74 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
opm/simulators/flow/FlowProblem.hpp
View File

@ -632,76 +632,102 @@ public:
void endTimeStep()
{
OPM_TIMEBLOCK(endTimeStep);
#ifndef NDEBUG
if constexpr (getPropValue<TypeTag, Properties::EnableDebuggingChecks>()) {
// in debug mode, we don't care about performance, so we check if the model does
// the right thing (i.e., the mass change inside the whole reservoir must be
// equivalent to the fluxes over the grid's boundaries plus the source rates
// specified by the problem)
int rank = this->simulator().gridView().comm().rank();
if (rank == 0)
// in debug mode, we don't care about performance, so we check
// if the model does the right thing (i.e., the mass change
// inside the whole reservoir must be equivalent to the fluxes
// over the grid's boundaries plus the source rates specified by
// the problem).
const int rank = this->simulator().gridView().comm().rank();
if (rank == 0) {
std::cout << "checking conservativeness of solution\n";
}
this->model().checkConservativeness(/*tolerance=*/-1, /*verbose=*/true);
if (rank == 0)
if (rank == 0) {
std::cout << "solution is sufficiently conservative\n";
}
}
#endif // NDEBUG
auto& simulator = this->simulator();
wellModel_.endTimeStep();
aquiferModel_.endTimeStep();
tracerModel_.endTimeStep();
this->wellModel_.endTimeStep();
this->aquiferModel_.endTimeStep();
this->tracerModel_.endTimeStep();
// Compute flux for output
this->model().linearizer().updateFlowsInfo();
// deal with DRSDT and DRVDT
asImp_().updateCompositionChangeLimits_();
{
this->asImp_().updateCompositionChangeLimits_();
if (this->enableDriftCompensation_) {
OPM_TIMEBLOCK(driftCompansation);
if (enableDriftCompensation_) {
const auto& residual = this->model().linearizer().residual();
for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
drift_[globalDofIdx] = residual[globalDofIdx];
drift_[globalDofIdx] *= simulator.timeStepSize();
if constexpr (getPropValue<TypeTag, Properties::UseVolumetricResidual>())
drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
this->drift_[globalDofIdx] = residual[globalDofIdx] * simulator.timeStepSize();
if constexpr (getPropValue<TypeTag, Properties::UseVolumetricResidual>()) {
this->drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
}
}
}
bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>() &&
!this->simulator().episodeWillBeOver();
const bool isSubStep = !Parameters::get<TypeTag, Properties::EnableWriteAllSolutions>()
&& !this->simulator().episodeWillBeOver();
// For CpGrid with LGRs, ecl/vtk output is not supported yet.
const auto& grid = this->simulator().vanguard().gridView().grid();
using GridType = std::remove_cv_t< typename std::remove_reference<decltype(grid)>::type>;
bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
if ( !isCpGrid || (this->simulator().vanguard().gridView().grid().maxLevel()==0)) {
eclWriter_->evalSummaryState(isSubStep);
using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
if (!isCpGrid || (grid.maxLevel() == 0)) {
this->eclWriter_->evalSummaryState(isSubStep);
}
int episodeIdx = this->episodeIndex();
// Re-ordering in case of Alugrid
std::function<unsigned int(unsigned int)> gridToEquilGrid = [&simulator](unsigned int i) {
return simulator.vanguard().gridIdxToEquilGridIdx(i);
};
std::function<void(bool)> transUp =
[this,gridToEquilGrid](bool global) {
this->transmissibilities_.update(global,gridToEquilGrid);
};
{
OPM_TIMEBLOCK(applyActions);
actionHandler_.applyActions(episodeIdx,
simulator.time() + simulator.timeStepSize(),
transUp);
const int episodeIdx = this->episodeIndex();
// Re-ordering in case of Alugrid
this->actionHandler_
.applyActions(episodeIdx, simulator.time() + simulator.timeStepSize(),
[this](const bool global)
{
this->transmissibilities_
.update(global, [&vg = this->simulator().vanguard()]
(const unsigned int i)
{
return vg.gridIdxToEquilGridIdx(i);
});
});
// Rerun UDQ assignents following action processing to make sure
// that any UDQ ASSIGN operations triggered in action blocks
// take effect. This is mainly to work around a shortcoming of
// the ScheduleState copy constructor which clears pending UDQ
// assignments under the assumption that all such assignments
// have been processed. If an action block happens to trigger
// on the final time step of an episode and that action block
// runs a UDQ assignment, then that assignment would be dropped
// and the rest of the simulator will never see its effect
// without this hack.
this->actionHandler_
.evalUDQAssignments(episodeIdx, this->simulator().vanguard().udqState());
}
// deal with "clogging" for the MICP model
if constexpr (enableMICP){
if constexpr (enableMICP) {
auto& model = this->model();
const auto& residual = this->model().linearizer().residual();
for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) {
const auto& residual = model.linearizer().residual();
for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); ++globalDofIdx) {
auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx];
MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
}