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

ebos: determine the upstream degree of freedom consistently

Browse Source 
before, if the pressure gradient was zero the interior DOF was assumed
to be the upstream one. On the other side of the face it was done the
same which meant that the upstream cell was different depending on
which cell was looked at. This did not have any effect on the value of
the flux (because the pressure gradient was zero anyway), but when AD
was used this resulted in non-symmetric derivatives. In principle this
is okay because the point where the pressure difference between cells
is zero is a kink and thus the flux derivatives there are
undefined. In practice this made comparing with flow quite difficult,
so let's change it...
This commit is contained in:
Andreas Lauser 2016-06-14 11:19:32 +02:00
parent 4f539c0d4d
commit deff4aab74
1 changed files with 33 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
applications/ebos/eclfluxmodule.hh
View File

@ -184,7 +184,8 @@ protected:
unsigned upstreamIndex_(unsigned phaseIdx) const
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return (pressureDifferential_[phaseIdx] >= 0)?exteriorDofIdx_:interiorDofIdx_;
return upIdx_[phaseIdx];
}
/*!
@ -197,7 +198,8 @@ protected:
unsigned downstreamIndex_(unsigned phaseIdx) const
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
return (pressureDifferential_[phaseIdx] < 0)?exteriorDofIdx_:interiorDofIdx_;
return dnIdx_[phaseIdx];
}
/*!
@ -249,6 +251,29 @@ protected:
pressureDifferential_[phaseIdx] = pressureExterior - pressureInterior;
// 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
// is always the same: if the pressure is equal, the DOF with the lower
// global index is regarded to be the upstream one.
if (pressureDifferential_[phaseIdx] == 0) {
if (I > J) {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
}
else {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
}
}
else if (pressureDifferential_[phaseIdx] > 0) {
upIdx_[phaseIdx] = exteriorDofIdx_;
dnIdx_[phaseIdx] = interiorDofIdx_;
}
else {
upIdx_[phaseIdx] = interiorDofIdx_;
dnIdx_[phaseIdx] = exteriorDofIdx_;
}
// apply the threshold pressure for the intersection. note that the concept
// of threshold pressure is a quite big hack that only makes sense for ECL
// datasets. (and even there its physical justification is quite
@ -283,10 +308,6 @@ protected:
void calculateFluxes_(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
{ }
// the local indices of the interior and exterior degrees of freedom
unsigned interiorDofIdx_;
unsigned exteriorDofIdx_;
// transmissibility [m^3 s]
Scalar trans_;
@ -302,6 +323,12 @@ protected:
// the difference in effective pressure between the two degrees of
// freedom [Pa]
Evaluation pressureDifferential_[numPhases];
// the local indices of the interior and exterior degrees of freedom
unsigned short interiorDofIdx_;
unsigned short exteriorDofIdx_;
unsigned short upIdx_[numPhases];
unsigned short dnIdx_[numPhases];
};
} // namespace Ewoms