From 4e5d85188c81748c5a458c951ba082933f0efeee Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 30 May 2023 16:15:19 +0200 Subject: [PATCH 1/5] Add and use shiftStorageCache(). --- .../common/fvbasediscretization.hh | 30 +++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/opm/models/discretization/common/fvbasediscretization.hh b/opm/models/discretization/common/fvbasediscretization.hh index 65f4d7566..d7fa9211c 100644 --- a/opm/models/discretization/common/fvbasediscretization.hh +++ b/opm/models/discretization/common/fvbasediscretization.hh @@ -878,6 +878,31 @@ public: return storageCache_[timeIdx][globalIdx]; } + /*! + * \brief Move the storage cache for a given time index to the back. + * + * This method should only be called by the time discretization. + * + * \param numSlots The number of time step slots for which the + * hints should be shifted. + */ + void shiftStorageCache(unsigned numSlots = 1) + { + if (!enableStorageCache()) { + return; + } + + if (simulator_.problem().recycleFirstIterationStorage()) { + return; + } + + assert(numSlots > 0); + + for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++ timeIdx) { + storageCache_[timeIdx + numSlots] = storageCache_[timeIdx]; + } + } + /*! * \brief Set an entry of the cache for the storage term. * @@ -1472,9 +1497,10 @@ public: // make the current solution the previous one. solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0); - // shift the intensive quantities cache by one position in the - // history + // shift the intensive quantities and storage caches by one + // position in the history (if active and required) asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1); + asImp_().shiftStorageCache(/*numSlots=*/1); } /*! From 5e96e4362ec18151709585fb0e30adfa404300ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 30 May 2023 16:01:58 +0200 Subject: [PATCH 2/5] Adapt TpfaLinearizer to possible shifting of storage cache. --- .../discretization/common/tpfalinearizer.hh | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 644227e28..02ceb287f 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -608,22 +608,18 @@ private: setResAndJacobi(res, bMat, adres); // Either use cached storage term, or compute it on the fly. if (model_().enableStorageCache()) { - // The cached storage for timeIdx 0 (current time) is not - // used, but after storage cache is shifted at the end of the - // timestep, it will become cached storage for timeIdx 1. - model_().updateCachedStorage(globI, /*timeIdx=*/0, res); - if (model_().newtonMethod().numIterations() == 0) { - // Need to update the storage cache. - if (problem_().recycleFirstIterationStorage()) { + if (problem_().recycleFirstIterationStorage()) { + if (model_().newtonMethod().numIterations() == 0) { // Assumes nothing have changed in the system which // affects masses calculated from primary variables. model_().updateCachedStorage(globI, /*timeIdx=*/1, res); - } else { - Dune::FieldVector tmp; - IntensiveQuantities intQuantOld = model_().intensiveQuantities(globI, 1); - LocalResidual::computeStorage(tmp, intQuantOld); - model_().updateCachedStorage(globI, /*timeIdx=*/1, tmp); } + } else { + // The cached storage for timeIdx 0 (current time) is + // not used at this step, but after storage cache is + // shifted at the end of the timestep, it will become + // cached storage for timeIdx 1. + model_().updateCachedStorage(globI, /*timeIdx=*/0, res); } res -= model_().cachedStorage(globI, 1); } else { From a35a5a00ac6ec2ac63ff1d7ea66d115707617b48 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 31 May 2023 15:06:28 +0200 Subject: [PATCH 3/5] Bugfix storage cache in FvBaseLocalResidual. --- .../common/fvbaselocalresidual.hh | 58 ++++++++----------- 1 file changed, 24 insertions(+), 34 deletions(-) diff --git a/opm/models/discretization/common/fvbaselocalresidual.hh b/opm/models/discretization/common/fvbaselocalresidual.hh index e8c4d0716..e4a4b9de9 100644 --- a/opm/models/discretization/common/fvbaselocalresidual.hh +++ b/opm/models/discretization/common/fvbaselocalresidual.hh @@ -535,40 +535,30 @@ protected: if (elemCtx.enableStorageCache()) { const auto& model = elemCtx.model(); unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - if (model.newtonMethod().numIterations() == 0 && - !elemCtx.haveStashedIntensiveQuantities()) - { - if (!elemCtx.problem().recycleFirstIterationStorage()) { - // we re-calculate the storage term for the solution of the - // previous time step from scratch instead of using the one of - // the first iteration of the current time step. - tmp2 = 0.0; - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/1); - asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/1); - } - else { - // if the storage term is cached and we're in the first iteration - // of the time step, use the storage term of the first iteration - // as the one as the solution of the last time step (this assumes - // that the initial guess for the solution at the end of the time - // step is the same as the solution at the beginning of the time - // step. This is usually true, but some fancy preprocessing - // scheme might invalidate that assumption.) - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]); - } - - Valgrind::CheckDefined(tmp2); - - model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); - } - else { - // if the mass storage at the beginning of the time step is not cached, - // if the storage term is cached and we're not looking at the first - // iteration of the time step, we take the cached data. - tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1); - Valgrind::CheckDefined(tmp2); - } + for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) { + tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]); + } + if(elemCtx.problem().recycleFirstIterationStorage()) { + if (model.newtonMethod().numIterations() == 0) + { + // if the storage term is cached and we're in the first iteration + // of the time step, use the storage term of the first iteration + // as the one as the solution of the last time step (this assumes + // that the initial guess for the solution at the end of the time + // step is the same as the solution at the beginning of the time + // step. This is usually true, but some fancy preprocessing + // scheme might invalidate that assumption.) + Valgrind::CheckDefined(tmp2); + model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); + } + } else { + model.updateCachedStorage(globalDofIdx, /*timeIdx=*/0, tmp2); + } + // if the storage term at the beginning of the time step is cached + // from the last time step or we're not looking at the first + // iteration of the time step, we take the cached data. + tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1); + Valgrind::CheckDefined(tmp2); } else { // if the mass storage at the beginning of the time step is not cached, From f65ea9410c0463d6f19e43f5754336d7bb6a134b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 2 Jun 2023 16:25:25 +0200 Subject: [PATCH 4/5] Fix for finite difference (not AD) usage. --- .../common/fvbaselocalresidual.hh | 39 ++++++++++--------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/opm/models/discretization/common/fvbaselocalresidual.hh b/opm/models/discretization/common/fvbaselocalresidual.hh index e4a4b9de9..eb99e91a0 100644 --- a/opm/models/discretization/common/fvbaselocalresidual.hh +++ b/opm/models/discretization/common/fvbaselocalresidual.hh @@ -535,25 +535,26 @@ protected: if (elemCtx.enableStorageCache()) { const auto& model = elemCtx.model(); unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0); - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) { - tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]); - } - if(elemCtx.problem().recycleFirstIterationStorage()) { - if (model.newtonMethod().numIterations() == 0) - { - // if the storage term is cached and we're in the first iteration - // of the time step, use the storage term of the first iteration - // as the one as the solution of the last time step (this assumes - // that the initial guess for the solution at the end of the time - // step is the same as the solution at the beginning of the time - // step. This is usually true, but some fancy preprocessing - // scheme might invalidate that assumption.) - Valgrind::CheckDefined(tmp2); - model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); - } - } else { - model.updateCachedStorage(globalDofIdx, /*timeIdx=*/0, tmp2); - } + for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) { + tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]); + } + if (!elemCtx.haveStashedIntensiveQuantities()) { + if (elemCtx.problem().recycleFirstIterationStorage()) { + if (model.newtonMethod().numIterations() == 0) { + // if the storage term is cached and we're in the first iteration + // of the time step, use the storage term of the first iteration + // as the one as the solution of the last time step (this assumes + // that the initial guess for the solution at the end of the time + // step is the same as the solution at the beginning of the time + // step. This is usually true, but some fancy preprocessing + // scheme might invalidate that assumption.) + Valgrind::CheckDefined(tmp2); + model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); + } + } else { + model.updateCachedStorage(globalDofIdx, /*timeIdx=*/0, tmp2); + } + } // if the storage term at the beginning of the time step is cached // from the last time step or we're not looking at the first // iteration of the time step, we take the cached data. From 84baba140505a8baee35d8278ffbb15f632f5748 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 5 Jun 2023 09:07:15 +0200 Subject: [PATCH 5/5] Set storage cache for timeIdx 1 on first step. --- opm/models/discretization/common/fvbaselocalresidual.hh | 5 +++++ opm/models/discretization/common/tpfalinearizer.hh | 5 +++++ 2 files changed, 10 insertions(+) diff --git a/opm/models/discretization/common/fvbaselocalresidual.hh b/opm/models/discretization/common/fvbaselocalresidual.hh index eb99e91a0..c9a133d5c 100644 --- a/opm/models/discretization/common/fvbaselocalresidual.hh +++ b/opm/models/discretization/common/fvbaselocalresidual.hh @@ -553,6 +553,11 @@ protected: } } else { model.updateCachedStorage(globalDofIdx, /*timeIdx=*/0, tmp2); + // If we are at the very first timestep, there has + // never been a storage cache shift, and we must also set the timeIdx 1 cache here. + if (elemCtx.simulator().timeStepIndex() == 0 && model.newtonMethod().numIterations() == 0) { + model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2); + } } } // if the storage term at the beginning of the time step is cached diff --git a/opm/models/discretization/common/tpfalinearizer.hh b/opm/models/discretization/common/tpfalinearizer.hh index 02ceb287f..23d750243 100644 --- a/opm/models/discretization/common/tpfalinearizer.hh +++ b/opm/models/discretization/common/tpfalinearizer.hh @@ -620,6 +620,11 @@ private: // shifted at the end of the timestep, it will become // cached storage for timeIdx 1. model_().updateCachedStorage(globI, /*timeIdx=*/0, res); + // If we are at the very first timestep, there has + // never been a storage cache shift, and we must also set the timeIdx 1 cache here. + if (simulator_().timeStepIndex() == 0 && model_().newtonMethod().numIterations() == 0) { + model_().updateCachedStorage(globI, /*timeIdx=*/1, res); + } } res -= model_().cachedStorage(globI, 1); } else {