Skip to content

Commit

Permalink
Fix mass exchange term.
Browse files Browse the repository at this point in the history 
Additionally, store separate well terms for free and solution tracers
  • Loading branch information
@svenn-t
svenn-t committed May 21, 2024
1 parent aae1ccb commit 587df5e
Show file tree
Hide file tree
Showing 4 changed files with 107 additions and 33 deletions.
10 changes: 10 additions & 0 deletions opm/simulators/flow/GenericTracerModel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,8 @@ class GenericTracerModel {
const std::string& name(int tracerIdx) const;
std::string fname(int tracerIdx) const;
std::string sname(int tracerIdx) const;
std::string wellfname(int tracerIdx) const;
std::string wellsname(int tracerIdx) const;


/*!
Expand All @@ -83,12 +85,18 @@ class GenericTracerModel {
*/
const std::map<std::pair<std::string, std::string>, double>&
getWellTracerRates() const {return wellTracerRate_;}
const std::map<std::pair<std::string, std::string>, double>&
getWellFreeTracerRates() const {return wellFreeTracerRate_;}
const std::map<std::pair<std::string, std::string>, double>&
getWellSolTracerRates() const {return wellSolTracerRate_;}

template<class Serializer>
void serializeOp(Serializer& serializer)
{
serializer(tracerConcentration_);
serializer(wellTracerRate_);
serializer(wellFreeTracerRate_);
serializer(wellSolTracerRate_);
}

protected:
Expand Down Expand Up @@ -126,6 +134,8 @@ class GenericTracerModel {

// <wellName, tracerIdx> -> wellRate
std::map<std::pair<std::string, std::string>, double> wellTracerRate_;
std::map<std::pair<std::string, std::string>, double> wellFreeTracerRate_;
std::map<std::pair<std::string, std::string>, double> wellSolTracerRate_;
/// \brief Function returning the cell centers
std::function<std::array<double,dimWorld>(int)> centroids_;
};
Expand Down
16 changes: 15 additions & 1 deletion opm/simulators/flow/GenericTracerModel_impl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,20 @@ sname(int tracerIdx) const
return this->eclState_.tracer()[tracerIdx].sname();
}

template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
wellfname(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].wellfname();
}

template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
wellsname(int tracerIdx) const
{
return this->eclState_.tracer()[tracerIdx].wellsname();
}

template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
double GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
currentConcentration_(const Well& eclWell, const std::string& name) const
Expand Down Expand Up @@ -376,7 +390,7 @@ template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar
bool GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b)
{
Scalar tolerance = 1e-2;
Scalar tolerance = 1e-6;
int maxIter = 100;

int verbosity = 0;
Expand Down
110 changes: 78 additions & 32 deletions opm/simulators/flow/TracerModel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,8 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
// resize free and solution volume storages
fVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->freeTracerConcentration_[tracerIdx].size());
sVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size());
dsVol_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size());
dfVol_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size());
}

// will be valid after we move out of tracerMatrix_
Expand Down Expand Up @@ -391,6 +393,8 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G

TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
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];
Expand All @@ -415,7 +419,8 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned I,
unsigned J)
unsigned J,
const Scalar dt)
{
if (tr.numTracer() == 0)
return;
Expand All @@ -426,6 +431,8 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
bool isUpS;
computeFreeFlux_(fFlux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
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) {
Expand Down Expand Up @@ -455,25 +462,27 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
const auto& eclWell = well.wellEcl();
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
this->wellTracerRate_[std::make_pair(eclWell.name(), this->name(tr.idx_[tIdx]))] = 0.0;
this->wellFreeTracerRate_[std::make_pair(eclWell.name(), this->wellfname(tr.idx_[tIdx]))] = 0.0;
this->wellSolTracerRate_[std::make_pair(eclWell.name(), this->wellsname(tr.idx_[tIdx]))] = 0.0;
}

std::vector<double> wtracer(tr.numTracer());
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
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 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_);

Scalar rate_s;
if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
}
else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
}
else {
rate_s = 0.0;
Expand All @@ -488,44 +497,58 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G

// Store _injector_ tracer rate for reporting
this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate_f*wtracer[tIdx];
this->wellFreeTracerRate_.at(std::make_pair(eclWell.name(),this->wellfname(tr.idx_[tIdx]))) += rate_f*wtracer[tIdx];
}
dfVol_[tr.phaseIdx_][I] -= rate_f * dt;
}
else if (rate_f < 0) {
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
// Free and solution tracer production
tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0];
tr.residual_[tIdx][I][1] -= rate_s * tr.concentration_[tIdx][I][1];

// Store _producer_ tracer rate for reporting
this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) +=
rate_f * tr.concentration_[tIdx][I][0] + rate_s * tr.concentration_[tIdx][I][1];
rate_f * tr.concentration_[tIdx][I][0];
this->wellFreeTracerRate_.at(std::make_pair(eclWell.name(),this->wellfname(tr.idx_[tIdx]))) +=
rate_f * tr.concentration_[tIdx][I][0];
}

dfVol_[tr.phaseIdx_][I] -= rate_f * dt;

// Derivative matrix for producer
(*tr.mat)[I][I][0][0] -= rate_f * variable<TracerEvaluation>(1.0, 0).derivative(0);
}
if (rate_s < 0) {
for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
tr.residual_[tIdx][I][1] -= rate_s * tr.concentration_[tIdx][I][1];

this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) +=
rate_s * tr.concentration_[tIdx][I][1];
this->wellSolTracerRate_.at(std::make_pair(eclWell.name(),this->wellsname(tr.idx_[tIdx]))) +=
rate_s * tr.concentration_[tIdx][I][1];
}
dsVol_[tr.phaseIdx_][I] -= rate_s * dt;

(*tr.mat)[I][I][1][1] -= rate_s * variable<TracerEvaluation>(1.0, 0).derivative(0);
}
}
}

template<class TrRe>
void assembleTracerEquationSource(TrRe& tr,
const Scalar scvVolume,
const Scalar dt,
unsigned I)
{
if (tr.numTracer() == 0)
return;

// Diff dissolved phase volume
Scalar sVol0 = computeSolutionVolume_(tr.phaseIdx_, I, 0) * scvVolume;
Scalar dsVol = sVol0 - sVol1_[tr.phaseIdx_][I];

const Scalar& dsVol = dsVol_[tr.phaseIdx_][I];
const Scalar& dfVol = dfVol_[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) {
tr.residual_[tIdx][I][0] += (dsVol / dt) * tr.concentration_[tIdx][I][0];
tr.residual_[tIdx][I][1] -= (dsVol / dt) * tr.concentration_[tIdx][I][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];
}
else {
tr.residual_[tIdx][I][0] += (dsVol / dt) * tr.concentration_[tIdx][I][1];
Expand All @@ -535,12 +558,12 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G

// Derivative matrix
if (dsVol >= 0) {
(*tr.mat)[I][I][0][0] += (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][1][0] -= (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][0][0] -= (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][1][0] += (dfVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
}
else {
(*tr.mat)[I][I][1][1] += (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][0][1] -= (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][0][1] += (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
(*tr.mat)[I][I][1][1] -= (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
}
}

Expand All @@ -561,6 +584,14 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
}
}

// Well terms
const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
for (const auto& wellPtr : wellPtrs) {
for (auto& tr : tbatch) {
this->assembleTracerEquationWell(tr, *wellPtr);
}
}

ElementContext elemCtx(simulator_);
for (const auto& elem : elements(simulator_.gridView())) {
elemCtx.updateStencil(elem);
Expand Down Expand Up @@ -603,25 +634,17 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
unsigned j = face.exteriorIndex();
unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
for (auto& tr : tbatch) {
this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J);
this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J, dt);
}
}

// Source terms (mass transfer between free and solution tracer)
for (auto& tr : tbatch) {
this->assembleTracerEquationSource(tr, scvVolume, dt, I);
this->assembleTracerEquationSource(tr, dt, I);
}

}

// Well terms
const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
for (const auto& wellPtr : wellPtrs) {
for (auto& tr : tbatch) {
this->assembleTracerEquationWell(tr, *wellPtr);
}
}

// Communicate overlap using grid Communication
for (auto& tr : tbatch) {
if (tr.numTracer() == 0)
Expand Down Expand Up @@ -656,6 +679,8 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
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;
dfVol_[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];
Expand Down Expand Up @@ -685,11 +710,30 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
}

for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
// New concetration
tr.concentration_[tIdx] -= dx[tIdx];

// Partition concentration into free and solution tracers for output
for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) {
// New concetration. Concentrations that are negative or where free/solution phase is not
// 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_));
Scalar Ss = 0.0;
if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas())
Ss = decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx));
else if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil())
Ss = decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx));

if ((tr.concentration_[tIdx][globalDofIdx][0] - dx[tIdx][globalDofIdx][0] < 0.0)
|| Sf < 1e-6)
tr.concentration_[tIdx][globalDofIdx][0] = 0.0;
else
tr.concentration_[tIdx][globalDofIdx][0] -= dx[tIdx][globalDofIdx][0];
if (tr.concentration_[tIdx][globalDofIdx][1] - dx[tIdx][globalDofIdx][1] < 0.0
|| Ss < 1e-6)
tr.concentration_[tIdx][globalDofIdx][1] = 0.0;
else
tr.concentration_[tIdx][globalDofIdx][1] -= dx[tIdx][globalDofIdx][1];

// 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];
}
Expand Down Expand Up @@ -805,6 +849,8 @@ class TracerModel : public GenericTracerModel<GetPropType<TypeTag, Properties::G
TracerBatch<TracerVector>& gas_;
std::array<std::vector<Scalar>, 3> fVol1_;
std::array<std::vector<Scalar>, 3> sVol1_;
std::array<std::vector<Scalar>, 3> dsVol_;
std::array<std::vector<Scalar>, 3> dfVol_;
};

} // namespace Opm
Expand Down
4 changes: 4 additions & 0 deletions opm/simulators/wells/BlackoilWellModel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,6 +266,10 @@ class WellContributions;

const auto& tracerRates = simulator_.problem().tracerModel().getWellTracerRates();
this->assignWellTracerRates(wsrpt, tracerRates);
const auto& freeTracerRates = simulator_.problem().tracerModel().getWellFreeTracerRates();
this->assignWellTracerRates(wsrpt, freeTracerRates);
const auto& solTracerRates = simulator_.problem().tracerModel().getWellSolTracerRates();
this->assignWellTracerRates(wsrpt, solTracerRates);


BlackoilWellModelGuideRates(*this).assignWellGuideRates(wsrpt, this->reportStepIndex());
Expand Down

0 comments on commit 587df5e

Please sign in to comment.