From 376de4c733bfd0dfba7b2a493ab03b87e9d12849 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 28 Feb 2025 09:27:13 +0100 Subject: [PATCH 1/5] cosmetics --- .../flow/GenericTracerModel_impl.hpp | 9 ++- opm/simulators/flow/TracerModel.hpp | 76 +++++++++---------- 2 files changed, 44 insertions(+), 41 deletions(-) diff --git a/opm/simulators/flow/GenericTracerModel_impl.hpp b/opm/simulators/flow/GenericTracerModel_impl.hpp index 5c7575d230..ac487eb5c6 100644 --- a/opm/simulators/flow/GenericTracerModel_impl.hpp +++ b/opm/simulators/flow/GenericTracerModel_impl.hpp @@ -116,8 +116,9 @@ template:: freeTracerConcentration(int tracerIdx, int globalDofIdx) const { - if (freeTracerConcentration_.empty()) + if (freeTracerConcentration_.empty()) { return 0.0; + } return freeTracerConcentration_[tracerIdx][globalDofIdx]; } @@ -126,8 +127,9 @@ template:: solTracerConcentration(int tracerIdx, int globalDofIdx) const { - if (solTracerConcentration_.empty()) + if (solTracerConcentration_.empty()) { return 0.0; + } return solTracerConcentration_[tracerIdx][globalDofIdx]; } @@ -225,8 +227,9 @@ doInit(bool rst, std::size_t numGridDof, { const auto& tracers = eclState_.tracer(); - if (tracers.size() == 0) + if (tracers.size() == 0) { return; // tracer treatment is supposed to be disabled + } // retrieve the number of tracers from the deck const std::size_t numTracers = tracers.size(); diff --git a/opm/simulators/flow/TracerModel.hpp b/opm/simulators/flow/TracerModel.hpp index de91537268..3d09e9104b 100644 --- a/opm/simulators/flow/TracerModel.hpp +++ b/opm/simulators/flow/TracerModel.hpp @@ -139,21 +139,24 @@ public: for (std::size_t tracerIdx = 0; tracerIdx < this->tracerPhaseIdx_.size(); ++tracerIdx) { if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::waterPhaseIdx) { if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){ - throw std::runtime_error("Water tracer specified for non-water fluid system:" + this->name(tracerIdx)); + throw std::runtime_error("Water tracer specified for non-water fluid system: " + + this->name(tracerIdx)); } wat_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); } else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::oilPhaseIdx) { if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){ - throw std::runtime_error("Oil tracer specified for non-oil fluid system:" + this->name(tracerIdx)); + throw std::runtime_error("Oil tracer specified for non-oil fluid system: " + + this->name(tracerIdx)); } oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); } else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::gasPhaseIdx) { if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){ - throw std::runtime_error("Gas tracer specified for non-gas fluid system:" + this->name(tracerIdx)); + throw std::runtime_error("Gas tracer specified for non-gas fluid system: " + + this->name(tracerIdx)); } gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]); @@ -241,10 +244,9 @@ protected: const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx); const auto& fs = intQuants.fluidState(); - const Scalar phaseVolume = - decay(fs.saturation(tracerPhaseIdx)) - *decay(fs.invB(tracerPhaseIdx)) - *decay(intQuants.porosity()); + const Scalar phaseVolume = decay(fs.saturation(tracerPhaseIdx)) * + decay(fs.invB(tracerPhaseIdx)) * + decay(intQuants.porosity()); return max(phaseVolume, 1e-10); } @@ -261,20 +263,18 @@ protected: // vaporized oil if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) { - phaseVolume = - decay(fs.saturation(FluidSystem::gasPhaseIdx)) - * decay(fs.invB(FluidSystem::gasPhaseIdx)) - * decay(fs.Rv()) - * decay(intQuants.porosity()); + phaseVolume = decay(fs.saturation(FluidSystem::gasPhaseIdx)) * + decay(fs.invB(FluidSystem::gasPhaseIdx)) * + decay(fs.Rv()) * + decay(intQuants.porosity()); } // dissolved gas else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) { - phaseVolume = - decay(fs.saturation(FluidSystem::oilPhaseIdx)) - * decay(fs.invB(FluidSystem::oilPhaseIdx)) - * decay(fs.Rs()) - * decay(intQuants.porosity()); + phaseVolume = decay(fs.saturation(FluidSystem::oilPhaseIdx)) * + decay(fs.invB(FluidSystem::oilPhaseIdx)) * + decay(fs.Rs()) * + decay(intQuants.porosity()); } else { phaseVolume = 0.0; @@ -302,9 +302,8 @@ protected: const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx); const auto& fs = intQuants.fluidState(); - const Scalar v = - decay(extQuants.volumeFlux(tracerPhaseIdx)) - * decay(fs.invB(tracerPhaseIdx)); + const Scalar v = decay(extQuants.volumeFlux(tracerPhaseIdx)) * + decay(fs.invB(tracerPhaseIdx)); const Scalar A = scvf.area(); if (inIdx == upIdx) { @@ -339,10 +338,9 @@ protected: const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx); const auto& fs = intQuants.fluidState(); - v = - decay(fs.invB(FluidSystem::gasPhaseIdx)) - * decay(extQuants.volumeFlux(FluidSystem::gasPhaseIdx)) - * decay(fs.Rv()); + v = decay(fs.invB(FluidSystem::gasPhaseIdx)) * + decay(extQuants.volumeFlux(FluidSystem::gasPhaseIdx)) * + decay(fs.Rv()); } // dissolved gas else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) { @@ -350,10 +348,9 @@ protected: const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx); const auto& fs = intQuants.fluidState(); - v = - decay(fs.invB(FluidSystem::oilPhaseIdx)) - * decay(extQuants.volumeFlux(FluidSystem::oilPhaseIdx)) - * decay(fs.Rs()); + v = decay(fs.invB(FluidSystem::oilPhaseIdx)) * + decay(extQuants.volumeFlux(FluidSystem::oilPhaseIdx)) * + decay(fs.Rs()); } else { upIdx = 0; @@ -525,7 +522,8 @@ protected: // Injection of free tracer only tr.residual_[tIdx][I][0] -= rate_f*wtracer[tIdx]; - // Store _injector_ tracer rate for reporting (can be done here since WTRACER is constant) + // Store _injector_ tracer rate for reporting + // (can be done here since WTRACER is constant) tracerRate[tIdx].rate += rate_f*wtracer[tIdx]; freeTracerRate[tIdx].rate += rate_f*wtracer[tIdx]; if (eclWell.isMultiSegment()) { @@ -536,13 +534,13 @@ protected: } else if (rate_f < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { - // Store _injector_ tracer rate for cross-flowing well connections (can be done here since WTRACER is constant) + // Store _injector_ tracer rate for cross-flowing well connections + // (can be done here since WTRACER is constant) tracerRate[tIdx].rate += rate_f*wtracer[tIdx]; freeTracerRate[tIdx].rate += rate_f*wtracer[tIdx]; // Production of free tracer tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0]; - } dVol_[0][tr.phaseIdx_][I] -= rate_f * dt; @@ -610,7 +608,7 @@ protected: // Note that we formulate the equations in terms of a concentration update // (compared to previous time step) and not absolute concentration. // This implies that current concentration (tr.concentration_[][]) contributes - // to the rhs both through storrage and flux terms. + // to the rhs both through storage and flux terms. // Compare also advanceTracerFields(...) below. DeferredLogger local_deferredLogger{}; @@ -870,17 +868,19 @@ protected: } } - // TODO: Some inconsistencies here that perhaps should be clarified. The "offical" rate as reported below is - // occasionally significant different from the sum over connections (as calculated above). Only observed - // for small values, neglible for the rate itself, but matters when used to calculate tracer concentrations. + // TODO: Some inconsistencies here that perhaps should be clarified. + // The "offical" rate as reported below is occasionally significant + // different from the sum over connections (as calculated above). Only observed + // for small values, neglible for the rate itself, but matters when used to + // calculate tracer concentrations. const Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_]; const Scalar rateWellTotal = official_well_rate_total; if (rateWellTotal > rateWellNeg) { // Cross flow - const Scalar bucketPrDay = 10.0/(1000.*3600.*24.); // ... keeps (some) trouble away - const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal/rateWellNeg : 0.0; + constexpr Scalar bucketPrDay = 10.0 / (1000. * 3600. * 24.); // ... keeps (some) trouble away + const Scalar factor = (rateWellTotal < -bucketPrDay) ? rateWellTotal / rateWellNeg : 0.0; for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { tracerRate[tIdx].rate *= factor; } @@ -940,7 +940,7 @@ protected: int numTracer() const { return idx_.size(); } - void addTracer(const int idx, const TV & concentration) + void addTracer(const int idx, const TV& concentration) { const int numGridDof = concentration.size(); idx_.emplace_back(idx); From ea67be3a3b255071edf3ea9da843f3dc27bb0c6c Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 28 Feb 2025 09:27:24 +0100 Subject: [PATCH 2/5] avoid repeated calculations --- opm/simulators/flow/TracerModel.hpp | 50 ++++++++++++++++------------- 1 file changed, 28 insertions(+), 22 deletions(-) diff --git a/opm/simulators/flow/TracerModel.hpp b/opm/simulators/flow/TracerModel.hpp index 3d09e9104b..b5b49bfcfa 100644 --- a/opm/simulators/flow/TracerModel.hpp +++ b/opm/simulators/flow/TracerModel.hpp @@ -519,25 +519,27 @@ protected: const Scalar rate_f = rate - rate_s; if (rate_f > 0) { 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] -= rate_f*wtracer[tIdx]; + tr.residual_[tIdx][I][0] -= delta; // Store _injector_ tracer rate for reporting // (can be done here since WTRACER is constant) - tracerRate[tIdx].rate += rate_f*wtracer[tIdx]; - freeTracerRate[tIdx].rate += rate_f*wtracer[tIdx]; + tracerRate[tIdx].rate += delta; + freeTracerRate[tIdx].rate += delta; if (eclWell.isMultiSegment()) { - (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += rate_f*wtracer[tIdx]; + (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta; } } dVol_[0][tr.phaseIdx_][I] -= rate_f * dt; } else if (rate_f < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { + const Scalar delta = rate_f * wtracer[tIdx]; // Store _injector_ tracer rate for cross-flowing well connections // (can be done here since WTRACER is constant) - tracerRate[tIdx].rate += rate_f*wtracer[tIdx]; - freeTracerRate[tIdx].rate += rate_f*wtracer[tIdx]; + tracerRate[tIdx].rate += delta; + freeTracerRate[tIdx].rate += delta; // Production of free tracer tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0]; @@ -583,23 +585,27 @@ protected: // 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) { - tr.residual_[tIdx][I][0] -= (dfVol / dt) * tr.concentration_[tIdx][I][0]; - tr.residual_[tIdx][I][1] += (dfVol / dt) * tr.concentration_[tIdx][I][0]; + const auto delta = (dfVol / dt) * tr.concentration_[tIdx][I][0]; + tr.residual_[tIdx][I][0] -= delta; + tr.residual_[tIdx][I][1] += delta; } else { - tr.residual_[tIdx][I][0] += (dsVol / dt) * tr.concentration_[tIdx][I][1]; - tr.residual_[tIdx][I][1] -= (dsVol / dt) * tr.concentration_[tIdx][I][1]; + const auto delta = (dsVol / dt) * tr.concentration_[tIdx][I][1]; + tr.residual_[tIdx][I][0] += delta; + tr.residual_[tIdx][I][1] -= delta; } } // Derivative matrix if (dsVol >= 0) { - (*tr.mat)[I][I][0][0] -= (dfVol / dt) * variable(1.0, 0).derivative(0); - (*tr.mat)[I][I][1][0] += (dfVol / dt) * variable(1.0, 0).derivative(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; } else { - (*tr.mat)[I][I][0][1] += (dsVol / dt) * variable(1.0, 0).derivative(0); - (*tr.mat)[I][I][1][1] -= (dsVol / dt) * variable(1.0, 0).derivative(0); + 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; } } @@ -840,22 +846,22 @@ protected: if (rate_f < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { // Store _producer_ free tracer rate for reporting - tracerRate[tIdx].rate += rate_f * tr.concentration_[tIdx][I][0]; - freeTracerRate[tIdx].rate += rate_f * tr.concentration_[tIdx][I][0]; + const Scalar delta = rate_f * tr.concentration_[tIdx][I][0]; + tracerRate[tIdx].rate += delta; + freeTracerRate[tIdx].rate += delta; if (eclWell.isMultiSegment()) { - (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += - rate_f * tr.concentration_[tIdx][I][0]; + (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta; } } } if (rate_s < 0) { for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) { // Store _producer_ solution tracer rate for reporting - tracerRate[tIdx].rate += rate_s * tr.concentration_[tIdx][I][1]; - solTracerRate[tIdx].rate += rate_s * tr.concentration_[tIdx][I][1]; + const Scalar delta = rate_s * tr.concentration_[tIdx][I][1]; + tracerRate[tIdx].rate += delta; + solTracerRate[tIdx].rate += delta; if (eclWell.isMultiSegment()) { - (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += - rate_s * tr.concentration_[tIdx][I][1]; + (*mswTracerRate)[tIdx].rate[eclWell.getConnections().get(i).segment()] += delta; } } } From dac7d41e70f21e00ab11c825b30e5a87187f92d1 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 28 Feb 2025 12:51:35 +0100 Subject: [PATCH 3/5] mark member functions const --- opm/simulators/flow/TracerModel.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/opm/simulators/flow/TracerModel.hpp b/opm/simulators/flow/TracerModel.hpp index b5b49bfcfa..de3c4d84b8 100644 --- a/opm/simulators/flow/TracerModel.hpp +++ b/opm/simulators/flow/TracerModel.hpp @@ -238,8 +238,8 @@ public: protected: // compute volume associated with free concentration Scalar computeFreeVolume_(const int tracerPhaseIdx, - unsigned globalDofIdx, - unsigned timeIdx) + const unsigned globalDofIdx, + const unsigned timeIdx) const { const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx); const auto& fs = intQuants.fluidState(); @@ -253,8 +253,8 @@ protected: // compute volume associated with solution concentration Scalar computeSolutionVolume_(const int tracerPhaseIdx, - unsigned globalDofIdx, - unsigned timeIdx) + const unsigned globalDofIdx, + const unsigned timeIdx) const { const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx); const auto& fs = intQuants.fluidState(); @@ -284,12 +284,12 @@ protected: } - void computeFreeFlux_(TracerEvaluation & freeFlux, - bool & isUp, + void computeFreeFlux_(TracerEvaluation& freeFlux, + bool& isUp, const int tracerPhaseIdx, const ElementContext& elemCtx, - unsigned scvfIdx, - unsigned timeIdx) + const unsigned scvfIdx, + const unsigned timeIdx) const { const auto& stencil = elemCtx.stencil(timeIdx); const auto& scvf = stencil.interiorFace(scvfIdx); @@ -320,8 +320,8 @@ protected: bool& isUp, const int tracerPhaseIdx, const ElementContext& elemCtx, - unsigned scvfIdx, - unsigned timeIdx) + const unsigned scvfIdx, + const unsigned timeIdx) const { const auto& stencil = elemCtx.stencil(timeIdx); const auto& scvf = stencil.interiorFace(scvfIdx); From c8d1e617f6da460b3a204589d4da543f28ce7673 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 5 Mar 2025 12:06:08 +0100 Subject: [PATCH 4/5] use ranged for --- opm/simulators/flow/GenericTracerModel_impl.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/opm/simulators/flow/GenericTracerModel_impl.hpp b/opm/simulators/flow/GenericTracerModel_impl.hpp index ac487eb5c6..a6e04b25f0 100644 --- a/opm/simulators/flow/GenericTracerModel_impl.hpp +++ b/opm/simulators/flow/GenericTracerModel_impl.hpp @@ -375,10 +375,8 @@ doInit(bool rst, std::size_t numGridDof, // all of its neighbors. (it also talks to itself since // degrees of freedom are sometimes quite egocentric.) for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) { - typename NeighborSet::iterator nIt = neighbors[dofIdx].begin(); - typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end(); - for (; nIt != nEndIt; ++nIt) { - tracerMatrix_->addindex(dofIdx, *nIt); + for (const auto& index : neighbors[dofIdx]) { + tracerMatrix_->addindex(dofIdx, index); } } tracerMatrix_->endindices(); From 0ac69fc4c5a0ab9fdc08e8377403b3bf8b44a92f Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 6 Mar 2025 09:40:22 +0100 Subject: [PATCH 5/5] avoid std::max call for water tracers --- opm/simulators/flow/TracerModel.hpp | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/opm/simulators/flow/TracerModel.hpp b/opm/simulators/flow/TracerModel.hpp index de3c4d84b8..cc9e5fe3bc 100644 --- a/opm/simulators/flow/TracerModel.hpp +++ b/opm/simulators/flow/TracerModel.hpp @@ -259,29 +259,27 @@ protected: const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx); const auto& fs = intQuants.fluidState(); - Scalar phaseVolume; + constexpr Scalar min_volume = 1e-10; // vaporized oil if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) { - phaseVolume = decay(fs.saturation(FluidSystem::gasPhaseIdx)) * - decay(fs.invB(FluidSystem::gasPhaseIdx)) * - decay(fs.Rv()) * - decay(intQuants.porosity()); + return std::max(decay(fs.saturation(FluidSystem::gasPhaseIdx)) * + decay(fs.invB(FluidSystem::gasPhaseIdx)) * + decay(fs.Rv()) * + decay(intQuants.porosity()), + min_volume); } // dissolved gas else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) { - phaseVolume = decay(fs.saturation(FluidSystem::oilPhaseIdx)) * - decay(fs.invB(FluidSystem::oilPhaseIdx)) * - decay(fs.Rs()) * - decay(intQuants.porosity()); - } - else { - phaseVolume = 0.0; + return std::max(decay(fs.saturation(FluidSystem::oilPhaseIdx)) * + decay(fs.invB(FluidSystem::oilPhaseIdx)) * + decay(fs.Rs()) * + decay(intQuants.porosity()), + min_volume); } - return max(phaseVolume, 1e-10); - + return min_volume; } void computeFreeFlux_(TracerEvaluation& freeFlux,