Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Tracer: some janitoring #6043

Merged
merged 2 commits into from
Mar 3, 2025
Merged

Tracer: some janitoring #6043

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fec172c
mark variables const where appropriate
akva2 Feb 27, 2025
224e705
remove unused variable
akva2 Feb 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 49 additions & 52 deletions opm/simulators/flow/TracerModel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -236,7 +236,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
const auto& fs = intQuants.fluidState();

Scalar phaseVolume =
const Scalar phaseVolume =
decay<Scalar>(fs.saturation(tracerPhaseIdx))
*decay<Scalar>(fs.invB(tracerPhaseIdx))
*decay<Scalar>(intQuants.porosity());
Expand Down Expand Up @@ -290,18 +290,18 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
const auto& scvf = stencil.interiorFace(scvfIdx);

const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned inIdx = extQuants.interiorIndex();
const unsigned inIdx = extQuants.interiorIndex();

unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
const unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);

const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
const auto& fs = intQuants.fluidState();

Scalar v =
const Scalar v =
decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx))
* decay<Scalar>(fs.invB(tracerPhaseIdx));

Scalar A = scvf.area();
const Scalar A = scvf.area();
if (inIdx == upIdx) {
freeFlux = A*v*variable<TracerEvaluation>(1.0, 0);
isUp = true;
Expand All @@ -323,7 +323,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
const auto& scvf = stencil.interiorFace(scvfIdx);

const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned inIdx = extQuants.interiorIndex();
const unsigned inIdx = extQuants.interiorIndex();

Scalar v;
unsigned upIdx;
Expand Down Expand Up @@ -355,7 +355,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
v = 0.0;
}

Scalar A = scvf.area();
const Scalar A = scvf.area();
if (inIdx == upIdx) {
solFlux = A*v*variable<TracerEvaluation>(1.0, 0);
isUp = true;
Expand All @@ -380,7 +380,6 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
}

// Storage terms at previous time step (timeIdx = 1)
std::vector<Scalar> storageOfTimeIndex1(tr.numTracer());
std::vector<Scalar> fStorageOfTimeIndex1(tr.numTracer());
std::vector<Scalar> sStorageOfTimeIndex1(tr.numTracer());
if (elemCtx.enableStorageCache()) {
Expand All @@ -390,28 +389,28 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
}
}
else {
Scalar fVolume1 = computeFreeVolume_(tr.phaseIdx_, I1, 1);
Scalar sVolume1 = computeSolutionVolume_(tr.phaseIdx_, I1, 1);
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];
}
}

TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
const TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
const TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
dsVol_[tr.phaseIdx_][I] += sVol.value() * scvVolume - sVol1_[tr.phaseIdx_][I];
dfVol_[tr.phaseIdx_][I] += fVol.value() * scvVolume - fVol1_[tr.phaseIdx_][I];
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
// Free part
Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][0];
Scalar fLocalStorage = (fStorageOfTimeIndex0 - fStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
tr.residual_[tIdx][I][0] += fLocalStorage; //residual + flux
const Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][0];
const Scalar fLocalStorage = (fStorageOfTimeIndex0 - fStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
tr.residual_[tIdx][I][0] += fLocalStorage; // residual + flux

// Solution part
Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][1];
Scalar sLocalStorage = (sStorageOfTimeIndex0 - sStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
tr.residual_[tIdx][I][1] += sLocalStorage; //residual + flux
const Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][1];
const Scalar sLocalStorage = (sStorageOfTimeIndex0 - sStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
tr.residual_[tIdx][I][1] += sLocalStorage; // residual + flux
}

// Derivative matrix
Expand Down Expand Up @@ -439,12 +438,12 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
computeSolFlux_(sFlux, isUpS, tr.phaseIdx_, elemCtx, scvfIdx, 0);
dsVol_[tr.phaseIdx_][I] += sFlux.value() * dt;
dfVol_[tr.phaseIdx_][I] += fFlux.value() * dt;
int fGlobalUpIdx = isUpF ? I : J;
int sGlobalUpIdx = isUpS ? I : J;
for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
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][0] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][0]; // residual + flux
tr.residual_[tIdx][I][1] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][1]; // residual + flux
}

// Derivative matrix
Expand Down Expand Up @@ -487,12 +486,12 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
wtracer[tIdx] = this->currentConcentration_(eclWell, this->name(tr.idx_[tIdx]));
}

Scalar dt = simulator_.timeStepSize();
std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value();
const Scalar dt = simulator_.timeStepSize();
const std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value();
const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
const auto I = ws.perf_data.cell_index[i];
Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
const Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
Scalar rate_s;
if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
Expand All @@ -504,7 +503,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
rate_s = 0.0;
}

Scalar rate_f = rate - rate_s;
const Scalar rate_f = rate - rate_s;
if (rate_f > 0) {
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
// Injection of free tracer only
Expand Down Expand Up @@ -619,10 +618,11 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
}

ElementContext elemCtx(simulator_);
const Scalar dt = elemCtx.simulator().timeStepSize();
for (const auto& elem : elements(simulator_.gridView())) {
elemCtx.updateStencil(elem);

std::size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/0);
const std::size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/0);

if (elem.partitionType() != Dune::InteriorEntity) {
// Dirichlet boundary conditions needed for the parallel matrix
Expand All @@ -637,27 +637,26 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
elemCtx.updateAllIntensiveQuantities();
elemCtx.updateAllExtensiveQuantities();

Scalar extrusionFactor =
const Scalar extrusionFactor =
elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
Valgrind::CheckDefined(extrusionFactor);
assert(isfinite(extrusionFactor));
assert(extrusionFactor > 0.0);
Scalar scvVolume =
const Scalar scvVolume =
elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume()
* extrusionFactor;
Scalar dt = elemCtx.simulator().timeStepSize();

std::size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/1);
const std::size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timeIdx=*/1);

for (auto& tr : tbatch) {
this->assembleTracerEquationVolume(tr, elemCtx, scvVolume, dt, I, I1);
}

std::size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
const std::size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
unsigned j = face.exteriorIndex();
unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
const unsigned j = face.exteriorIndex();
const unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
for (auto& tr : tbatch) {
this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt);
}
Expand All @@ -667,7 +666,6 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
for (auto& tr : tbatch) {
this->assembleTracerEquationSource(tr, dt, I);
}

}

// Communicate overlap using grid Communication
Expand All @@ -694,16 +692,16 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
for (const auto& elem : elements(simulator_.gridView())) {
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
Scalar extrusionFactor = elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0);
const Scalar extrusionFactor = elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
const Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
const int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0);
for (auto& tr : tbatch) {
if (tr.numTracer() == 0) {
continue;
}

Scalar fVol1 = computeFreeVolume_(tr.phaseIdx_, globalDofIdx, 0);
Scalar sVol1 = computeSolutionVolume_(tr.phaseIdx_, globalDofIdx, 0);
const Scalar fVol1 = computeFreeVolume_(tr.phaseIdx_, globalDofIdx, 0);
const Scalar sVol1 = computeSolutionVolume_(tr.phaseIdx_, globalDofIdx, 0);
fVol1_[tr.phaseIdx_][globalDofIdx] = fVol1 * scvVolume;
sVol1_[tr.phaseIdx_][globalDofIdx] = sVol1 * scvVolume;
dsVol_[tr.phaseIdx_][globalDofIdx] = 0.0;
Expand Down Expand Up @@ -732,7 +730,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
dx[tIdx] = 0.0;
}

bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
const bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
if (!converged) {
OpmLog::warning("### Tracer model: Linear solver did not converge. ###");
}
Expand All @@ -745,7 +743,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
// present are set to zero
const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, 0);
const auto& fs = intQuants.fluidState();
Scalar Sf = decay<Scalar>(fs.saturation(tr.phaseIdx_));
const Scalar Sf = decay<Scalar>(fs.saturation(tr.phaseIdx_));
Scalar Ss = 0.0;
if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
Ss = decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx));
Expand All @@ -754,7 +752,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
Ss = decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx));
}

const Scalar tol_gas_sat = 1e-6;
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;
}
Expand Down Expand Up @@ -786,11 +784,11 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G

Scalar rateWellPos = 0.0;
Scalar rateWellNeg = 0.0;
std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value();
const std::size_t well_index = simulator_.problem().wellModel().wellState().index(eclWell.name()).value();
const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
const auto I = ws.perf_data.cell_index[i];
Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
const Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);

Scalar rate_s;
if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
Expand All @@ -803,7 +801,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
rate_s = 0.0;
}

Scalar rate_f = rate - rate_s;
const Scalar rate_f = rate - rate_s;
if (rate_f < 0) {
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
// Store _producer_ free tracer rate for reporting
Expand Down Expand Up @@ -843,14 +841,13 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
}
}

Scalar rateWellTotal = rateWellNeg + rateWellPos;

// 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.
Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
const Scalar official_well_rate_total =
simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];

rateWellTotal = official_well_rate_total;
const Scalar rateWellTotal = official_well_rate_total;

if (rateWellTotal > rateWellNeg) { // Cross flow
const Scalar bucketPrDay = 10.0/(1000.*3600.*24.); // ... keeps (some) trouble away
Expand Down Expand Up @@ -916,7 +913,7 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G

void addTracer(const int idx, const TV & concentration)
{
int numGridDof = concentration.size();
const int numGridDof = concentration.size();
idx_.emplace_back(idx);
concentrationInitial_.emplace_back(concentration);
concentration_.emplace_back(concentration);
Expand Down