OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-09-05 04:40:19 -05:00
Code Issues Packages Projects Releases Wiki Activity

[WIP] Refactor calculateBoundaryGradients_()

Browse Source 
Not addressing solvent yet.
This commit is contained in:
Atgeirr Flø Rasmussen
2022-09-28 08:47:04 +02:00
parent b65f5e9808
commit 68e1479caf
1 changed files with 59 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files
ebos/eclfluxmodule.hh
+59 -32
View File
@@ -418,34 +418,61 @@ protected:
unsigned timeIdx, unsigned timeIdx,
const FluidState& exFluidState) const FluidState& exFluidState)
{ {
const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
Scalar faceArea = scvf.area();
Scalar zEx = scvf.integrationPos()[dimWorld - 1];
const auto& problem = elemCtx.problem(); const auto& problem = elemCtx.problem();
calculateBoundaryGradients_(problem,
elemCtx.globalSpaceIndex(0, timeIdx),
elemCtx.intensiveQuantities(0, timeIdx),
scvfIdx,
timeIdx,
faceArea,
zEx,
exFluidState,
upIdx_,
dnIdx_,
volumeFlux_,
pressureDifference_);
}
public:
/*!
* \brief Update the required gradients for boundary faces
*/
template <class Problem, class FluidState, class EvaluationContainer>
static void calculateBoundaryGradients_(const Problem& problem,
const unsigned globalSpaceIdx,
const IntensiveQuantities& intQuantsIn,
const unsigned bfIdx,
const unsigned timeIdx,
const double faceArea,
const double zEx,
const FluidState& exFluidState,
short (&upIdx)[numPhases],
short (&dnIdx)[numPhases],
EvaluationContainer& volumeFlux,
EvaluationContainer& pressureDifference)
{
bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions(); bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
if (!enableBoundaryMassFlux) if (!enableBoundaryMassFlux)
return; return;
const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.boundaryFace(scvfIdx);
unsigned interiorDofIdx = scvf.interiorIndex(); Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx);
Scalar trans = problem.transmissibilityBoundary(elemCtx, scvfIdx);
Scalar faceArea = scvf.area();
// estimate the gravity correction: for performance reasons we use a simplified // estimate the gravity correction: for performance reasons we use a simplified
// approach for this flux module that assumes that gravity is constant and always // approach for this flux module that assumes that gravity is constant and always
// acts into the downwards direction. (i.e., no centrifuge experiments, sorry.) // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
Scalar g = elemCtx.problem().gravity()[dimWorld - 1]; Scalar g = problem.gravity()[dimWorld - 1];
const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
// this is quite hacky because the dune grid interface does not provide a // this is quite hacky because the dune grid interface does not provide a
// cellCenterDepth() method (so we ask the problem to provide it). The "good" // cellCenterDepth() method (so we ask the problem to provide it). The "good"
// solution would be to take the Z coordinate of the element centroids, but since // solution would be to take the Z coordinate of the element centroids, but since
// ECL seems to like to be inconsistent on that front, it needs to be done like // ECL seems to like to be inconsistent on that front, it needs to be done like
// here... // here...
Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx); Scalar zIn = problem.dofCenterDepth(globalSpaceIdx);
Scalar zEx = scvf.integrationPos()[dimWorld - 1];
// the distances from the DOF's depths. (i.e., the additional depth of the // the distances from the DOF's depths. (i.e., the additional depth of the
// exterior DOF) // exterior DOF)
@@ -465,62 +492,62 @@ protected:
Evaluation pressureExterior = exFluidState.pressure(phaseIdx); Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
pressureExterior += rhoAvg*(distZ*g); pressureExterior += rhoAvg*(distZ*g);
pressureDifference_[phaseIdx] = pressureExterior - pressureInterior; pressureDifference[phaseIdx] = pressureExterior - pressureInterior;
// decide the upstream index for the phase. for this we make sure that the // decide the upstream index for the phase. for this we make sure that the
// degree of freedom which is regarded upstream if both pressures are equal // degree of freedom which is regarded upstream if both pressures are equal
// is always the same: if the pressure is equal, the DOF with the lower // is always the same: if the pressure is equal, the DOF with the lower
// global index is regarded to be the upstream one. // global index is regarded to be the upstream one.
if (pressureDifference_[phaseIdx] > 0.0) { const unsigned interiorDofIdx = 0; // Valid only for cell-centered FV.
upIdx_[phaseIdx] = -1; if (pressureDifference[phaseIdx] > 0.0) {
dnIdx_[phaseIdx] = interiorDofIdx; upIdx[phaseIdx] = -1;
dnIdx[phaseIdx] = interiorDofIdx;
} }
else { else {
upIdx_[phaseIdx] = interiorDofIdx; upIdx[phaseIdx] = interiorDofIdx;
dnIdx_[phaseIdx] = -1; dnIdx[phaseIdx] = -1;
} }
Evaluation transModified = trans; Evaluation transModified = trans;
unsigned upstreamIdx = upstreamIndex_(phaseIdx); if (upIdx[phaseIdx] == interiorDofIdx) {
if (upstreamIdx == interiorDofIdx) {
// this is slightly hacky because in the automatic differentiation case, it // this is slightly hacky because in the automatic differentiation case, it
// only works for the element centered finite volume method. for ebos this // only works for the element centered finite volume method. for ebos this
// does not matter, though. // does not matter, though.
const auto& up = elemCtx.intensiveQuantities(upstreamIdx, timeIdx); const auto& up = intQuantsIn;
// deal with water induced rock compaction // deal with water induced rock compaction
const double transMult = Toolbox::value(up.rockCompTransMultiplier()); const double transMult = Toolbox::value(up.rockCompTransMultiplier());
transModified *= transMult; transModified *= transMult;
volumeFlux_[phaseIdx] = volumeFlux[phaseIdx] =
pressureDifference_[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea); pressureDifference[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea);
if (enableSolvent && phaseIdx == gasPhaseIdx) // TODO: Figure out if this did have any effect. It should?
asImp_().setSolventVolumeFlux( pressureDifference_[phaseIdx]*up.solventMobility()*(-transModified/faceArea)); // if (enableSolvent && phaseIdx == gasPhaseIdx)
// asImp_().setSolventVolumeFlux( pressureDifference[phaseIdx]*up.solventMobility()*(-transModified/faceArea));
} }
else { else {
// compute the phase mobility using the material law parameters of the // compute the phase mobility using the material law parameters of the
// interior element. TODO: this could probably be done more efficiently // interior element. TODO: this could probably be done more efficiently
const auto& matParams = const auto& matParams = problem.materialLawParams(globalSpaceIdx);
elemCtx.problem().materialLawParams(elemCtx,
interiorDofIdx,
/*timeIdx=*/0);
std::array<typename FluidState::Scalar,numPhases> kr; std::array<typename FluidState::Scalar,numPhases> kr;
MaterialLaw::relativePermeabilities(kr, matParams, exFluidState); MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);
const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx); const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
volumeFlux_[phaseIdx] = volumeFlux[phaseIdx] =
pressureDifference_[phaseIdx]*mob*(-transModified/faceArea); pressureDifference[phaseIdx]*mob*(-transModified/faceArea);
// Solvent inflow is not yet supported // Solvent inflow is not yet supported
if (enableSolvent && phaseIdx == gasPhaseIdx) // if (enableSolvent && phaseIdx == gasPhaseIdx)
asImp_().setSolventVolumeFlux(0.0); // asImp_().setSolventVolumeFlux(0.0);
} }
} }
} }
protected:
/*! /*!
* \brief Update the volumetric fluxes for all fluid phases on the interior faces of the context * \brief Update the volumetric fluxes for all fluid phases on the interior faces of the context
*/ */