From 71a96821bcba127edf46332f8c17036cefb8985b Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 28 Feb 2025 13:17:55 +0100 Subject: [PATCH] enumerate tracer components --- opm/simulators/flow/GenericTracerModel.hpp | 6 ++ opm/simulators/flow/TracerModel.hpp | 120 +++++++++++---------- 2 files changed, 70 insertions(+), 56 deletions(-) diff --git a/opm/simulators/flow/GenericTracerModel.hpp b/opm/simulators/flow/GenericTracerModel.hpp index 0e59a7b306..ac954be18b 100644 --- a/opm/simulators/flow/GenericTracerModel.hpp +++ b/opm/simulators/flow/GenericTracerModel.hpp @@ -140,6 +140,12 @@ protected: Scalar currentConcentration_(const Well& eclWell, const std::string& name) const; + //! \brief Tracer type index + enum TracerTypeIdx { + Free = 0, + Solution = 1, + }; + const GridView& gridView_; const EclipseState& eclState_; const CartesianIndexMapper& cartMapper_; diff --git a/opm/simulators/flow/TracerModel.hpp b/opm/simulators/flow/TracerModel.hpp index cc9e5fe3bc..6a8aa6192b 100644 --- a/opm/simulators/flow/TracerModel.hpp +++ b/opm/simulators/flow/TracerModel.hpp @@ -236,6 +236,10 @@ public: } protected: + using TracerTypeIdx = typename BaseType::TracerTypeIdx; + using BaseType::Free; + using BaseType::Solution; + // compute volume associated with free concentration Scalar computeFreeVolume_(const int tracerPhaseIdx, const unsigned globalDofIdx, @@ -384,38 +388,38 @@ protected: std::vector sStorageOfTimeIndex1(tr.numTracer()); if (elemCtx.enableStorageCache()) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - fStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][0]; - sStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][1]; + fStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][Free]; + sStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][Solution]; } } else { const Scalar fVolume1 = computeFreeVolume_(tr.phaseIdx_, I1, 1); const Scalar sVolume1 = computeSolutionVolume_(tr.phaseIdx_, I1, 1); for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - fStorageOfTimeIndex1[tIdx] = fVolume1 * tr.concentration_[tIdx][I1][0]; - sStorageOfTimeIndex1[tIdx] = sVolume1 * tr.concentration_[tIdx][I1][1]; + fStorageOfTimeIndex1[tIdx] = fVolume1 * tr.concentration_[tIdx][I1][Free]; + sStorageOfTimeIndex1[tIdx] = sVolume1 * tr.concentration_[tIdx][I1][Solution]; } } const TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable(1.0, 0); const TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable(1.0, 0); - dVol_[1][tr.phaseIdx_][I] += sVol.value() * scvVolume - vol1_[1][tr.phaseIdx_][I]; - dVol_[0][tr.phaseIdx_][I] += fVol.value() * scvVolume - vol1_[0][tr.phaseIdx_][I]; + dVol_[Solution][tr.phaseIdx_][I] += sVol.value() * scvVolume - vol1_[1][tr.phaseIdx_][I]; + dVol_[Free][tr.phaseIdx_][I] += fVol.value() * scvVolume - vol1_[0][tr.phaseIdx_][I]; for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { // Free part - const Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][0]; + const Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][Free]; const Scalar fLocalStorage = (fStorageOfTimeIndex0 - fStorageOfTimeIndex1[tIdx]) * scvVolume/dt; - tr.residual_[tIdx][I][0] += fLocalStorage; // residual + flux + tr.residual_[tIdx][I][Free] += fLocalStorage; // residual + flux // Solution part - const Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][1]; + const Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][Solution]; const Scalar sLocalStorage = (sStorageOfTimeIndex0 - sStorageOfTimeIndex1[tIdx]) * scvVolume/dt; - tr.residual_[tIdx][I][1] += sLocalStorage; // residual + flux + tr.residual_[tIdx][I][Solution] += sLocalStorage; // residual + flux } // Derivative matrix - (*tr.mat)[I][I][0][0] += fVol.derivative(0) * scvVolume/dt; - (*tr.mat)[I][I][1][1] += sVol.derivative(0) * scvVolume/dt; + (*tr.mat)[I][I][Free][Free] += fVol.derivative(0) * scvVolume/dt; + (*tr.mat)[I][I][Solution][Solution] += sVol.derivative(0) * scvVolume/dt; } template @@ -436,24 +440,24 @@ protected: bool isUpS; computeFreeFlux_(fFlux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0); computeSolFlux_(sFlux, isUpS, tr.phaseIdx_, elemCtx, scvfIdx, 0); - dVol_[1][tr.phaseIdx_][I] += sFlux.value() * dt; - dVol_[0][tr.phaseIdx_][I] += fFlux.value() * dt; + dVol_[Solution][tr.phaseIdx_][I] += sFlux.value() * dt; + dVol_[Free][tr.phaseIdx_][I] += fFlux.value() * dt; const int fGlobalUpIdx = isUpF ? I : J; const int sGlobalUpIdx = isUpS ? I : J; for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { // Free and solution fluxes - tr.residual_[tIdx][I][0] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][0]; // residual + flux - tr.residual_[tIdx][I][1] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][1]; // residual + flux + tr.residual_[tIdx][I][Free] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][Free]; // residual + flux + tr.residual_[tIdx][I][Solution] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][Solution]; // residual + flux } // Derivative matrix if (isUpF){ - (*tr.mat)[J][I][0][0] = -fFlux.derivative(0); - (*tr.mat)[I][I][0][0] += fFlux.derivative(0); + (*tr.mat)[J][I][Free][Free] = -fFlux.derivative(0); + (*tr.mat)[I][I][Free][Free] += fFlux.derivative(0); } if (isUpS) { - (*tr.mat)[J][I][1][1] = -sFlux.derivative(0); - (*tr.mat)[I][I][1][1] += sFlux.derivative(0); + (*tr.mat)[J][I][Solution][Solution] = -sFlux.derivative(0); + (*tr.mat)[I][I][Solution][Solution] += sFlux.derivative(0); } } @@ -519,7 +523,7 @@ protected: for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { const Scalar delta = rate_f * wtracer[tIdx]; // Injection of free tracer only - tr.residual_[tIdx][I][0] -= delta; + tr.residual_[tIdx][I][Free] -= delta; // Store _injector_ tracer rate for reporting // (can be done here since WTRACER is constant) @@ -529,7 +533,7 @@ protected: (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta; } } - dVol_[0][tr.phaseIdx_][I] -= rate_f * dt; + dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt; } else if (rate_f < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { @@ -540,22 +544,22 @@ protected: freeTracerRate[tIdx].rate += delta; // Production of free tracer - tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0]; + tr.residual_[tIdx][I][Free] -= rate_f * tr.concentration_[tIdx][I][Free]; } - dVol_[0][tr.phaseIdx_][I] -= rate_f * dt; + dVol_[Free][tr.phaseIdx_][I] -= rate_f * dt; // Derivative matrix for free tracer producer - (*tr.mat)[I][I][0][0] -= rate_f * variable(1.0, 0).derivative(0); + (*tr.mat)[I][I][Free][Free] -= rate_f * variable(1.0, 0).derivative(0); } if (rate_s < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { // Production of solution tracer - tr.residual_[tIdx][I][1] -= rate_s * tr.concentration_[tIdx][I][1]; + tr.residual_[tIdx][I][Solution] -= rate_s * tr.concentration_[tIdx][I][Solution]; } - dVol_[1][tr.phaseIdx_][I] -= rate_s * dt; + dVol_[Solution][tr.phaseIdx_][I] -= rate_s * dt; // Derivative matrix for solution tracer producer - (*tr.mat)[I][I][1][1] -= rate_s * variable(1.0, 0).derivative(0); + (*tr.mat)[I][I][Solution][Solution] -= rate_s * variable(1.0, 0).derivative(0); } } } @@ -577,33 +581,33 @@ protected: return; } - const Scalar& dsVol = dVol_[1][tr.phaseIdx_][I]; - const Scalar& dfVol = dVol_[0][tr.phaseIdx_][I]; + const Scalar& dsVol = dVol_[Solution][tr.phaseIdx_][I]; + const Scalar& dfVol = dVol_[Free][tr.phaseIdx_][I]; // Source term determined by sign of dsVol: if dsVol > 0 then ms -> mf, else mf -> ms for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { if (dsVol >= 0) { - const auto delta = (dfVol / dt) * tr.concentration_[tIdx][I][0]; - tr.residual_[tIdx][I][0] -= delta; - tr.residual_[tIdx][I][1] += delta; + const auto delta = (dfVol / dt) * tr.concentration_[tIdx][I][Free]; + tr.residual_[tIdx][I][Free] -= delta; + tr.residual_[tIdx][I][Solution] += delta; } else { - const auto delta = (dsVol / dt) * tr.concentration_[tIdx][I][1]; - tr.residual_[tIdx][I][0] += delta; - tr.residual_[tIdx][I][1] -= delta; + const auto delta = (dsVol / dt) * tr.concentration_[tIdx][I][Solution]; + tr.residual_[tIdx][I][Free] += delta; + tr.residual_[tIdx][I][Solution] -= delta; } } // Derivative matrix if (dsVol >= 0) { const auto delta = (dfVol / dt) * variable(1.0, 0).derivative(0); - (*tr.mat)[I][I][0][0] -= delta; - (*tr.mat)[I][I][1][0] += delta; + (*tr.mat)[I][I][Free][Free] -= delta; + (*tr.mat)[I][I][Solution][Free] += delta; } else { const auto delta = (dsVol / dt) * variable(1.0, 0).derivative(0); - (*tr.mat)[I][I][0][1] += delta; - (*tr.mat)[I][I][1][1] -= delta; + (*tr.mat)[I][I][Free][Solution] += delta; + (*tr.mat)[I][I][Solution][Solution] -= delta; } } @@ -737,13 +741,13 @@ protected: const Scalar fVol1 = computeFreeVolume_(tr.phaseIdx_, globalDofIdx, 0); const Scalar sVol1 = computeSolutionVolume_(tr.phaseIdx_, globalDofIdx, 0); - vol1_[0][tr.phaseIdx_][globalDofIdx] = fVol1 * scvVolume; - vol1_[1][tr.phaseIdx_][globalDofIdx] = sVol1 * scvVolume; - dVol_[0][tr.phaseIdx_][globalDofIdx] = 0.0; - dVol_[1][tr.phaseIdx_][globalDofIdx] = 0.0; + vol1_[Free][tr.phaseIdx_][globalDofIdx] = fVol1 * scvVolume; + vol1_[Solution][tr.phaseIdx_][globalDofIdx] = sVol1 * scvVolume; + dVol_[Free][tr.phaseIdx_][globalDofIdx] = 0.0; + dVol_[Solution][tr.phaseIdx_][globalDofIdx] = 0.0; for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - tr.storageOfTimeIndex1_[tIdx][globalDofIdx][0] = fVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][0]; - tr.storageOfTimeIndex1_[tIdx][globalDofIdx][1] = sVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][1]; + tr.storageOfTimeIndex1_[tIdx][globalDofIdx][Free] = fVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][Free]; + tr.storageOfTimeIndex1_[tIdx][globalDofIdx][Solution] = sVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][Solution]; } } } @@ -788,22 +792,26 @@ protected: } constexpr Scalar tol_gas_sat = 1e-6; - if (tr.concentration_[tIdx][globalDofIdx][0] - dx[tIdx][globalDofIdx][0] < 0.0|| Sf < tol_gas_sat) { - tr.concentration_[tIdx][globalDofIdx][0] = 0.0; + if (tr.concentration_[tIdx][globalDofIdx][Free] - dx[tIdx][globalDofIdx][Free] < 0.0 || + Sf < tol_gas_sat) + { + tr.concentration_[tIdx][globalDofIdx][Free] = 0.0; } else { - tr.concentration_[tIdx][globalDofIdx][0] -= dx[tIdx][globalDofIdx][0]; + tr.concentration_[tIdx][globalDofIdx][Free] -= dx[tIdx][globalDofIdx][Free]; } - if (tr.concentration_[tIdx][globalDofIdx][1] - dx[tIdx][globalDofIdx][1] < 0.0 || Ss < tol_gas_sat) { - tr.concentration_[tIdx][globalDofIdx][1] = 0.0; + if (tr.concentration_[tIdx][globalDofIdx][Solution] - dx[tIdx][globalDofIdx][Solution] < 0.0 || + Ss < tol_gas_sat) + { + tr.concentration_[tIdx][globalDofIdx][Solution] = 0.0; } else { - tr.concentration_[tIdx][globalDofIdx][1] -= dx[tIdx][globalDofIdx][1]; + tr.concentration_[tIdx][globalDofIdx][Solution] -= dx[tIdx][globalDofIdx][Solution]; } // Partition concentration into free and solution tracers for output - this->freeTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][0]; - this->solTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][1]; + this->freeTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][Free]; + this->solTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][Solution]; } } @@ -844,7 +852,7 @@ protected: if (rate_f < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { // Store _producer_ free tracer rate for reporting - const Scalar delta = rate_f * tr.concentration_[tIdx][I][0]; + const Scalar delta = rate_f * tr.concentration_[tIdx][I][Free]; tracerRate[tIdx].rate += delta; freeTracerRate[tIdx].rate += delta; if (eclWell.isMultiSegment()) { @@ -855,7 +863,7 @@ protected: if (rate_s < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { // Store _producer_ solution tracer rate for reporting - const Scalar delta = rate_s * tr.concentration_[tIdx][I][1]; + const Scalar delta = rate_s * tr.concentration_[tIdx][I][Solution]; tracerRate[tIdx].rate += delta; solTracerRate[tIdx].rate += delta; if (eclWell.isMultiSegment()) {