OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-11-02 20:30:55 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6058 from akva2/tracer_enumerate

Browse Source 
Tracer: enumerate components
This commit is contained in:
Bård Skaflestad authored and GitHub committed 2025-03-06 17:33:45 +01:00
2 parents 1498ef31cb 71a96821bc
commit 5baca97503
2 files changed
+70 -56
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

opm/simulators/flow/GenericTracerModel.hpp
+6
View File
@@ -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_;
opm/simulators/flow/TracerModel.hpp
+64 -56
View File
@@ -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<Scalar> 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<TracerEvaluation>(1.0, 0);
const TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(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<class TrRe>
@@ -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<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][Free][Free] -= rate_f * variable<TracerEvaluation>(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<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][Solution][Solution] -= rate_s * variable<TracerEvaluation>(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<TracerEvaluation>(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<TracerEvaluation>(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()) {