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

Do not make changes to eclthresholdpressure.hh and ecltracermodel.hh.

Browse Source 
Not required for first version, but should be added later.
This commit is contained in:
Atgeirr Flø Rasmussen
2022-08-09 09:19:26 +02:00
parent 322dcbfb36
commit 727bf8d01f
2 changed files with 27 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

60
ebos/eclthresholdpressure.hh
View File

@@ -61,15 +61,11 @@ class EclThresholdPressure : public EclGenericThresholdPressure<GetPropType<Type
GetPropType<TypeTag, Properties::GridView>,
GetPropType<TypeTag, Properties::ElementMapper>,
GetPropType<TypeTag, Properties::Scalar>>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
enum { dimWorld = GridView::dimensionworld };
enum { enableExperiments = getPropValue<TypeTag, Properties::EnableExperiments>() };
enum { numPhases = FluidSystem::numPhases };
@@ -99,36 +95,6 @@ public:
}
private:
template <class Face, class Stencil, class ElemCtx>
double calculateMaxDp(Face& face,
Stencil& stencil,
ElemCtx& elemCtx,
const unsigned& scvfIdx,
const unsigned& i,
const unsigned& j,
const unsigned& insideElemIdx,
const unsigned& outsideElemIdx)
{
typedef MathToolbox<Evaluation> Toolbox;
elemCtx.updateIntensiveQuantities(/*timeIdx=*/0);
elemCtx.updateExtensiveQuantities(/*timeIdx=*/0);
// determine the maximum difference of the pressure of any phase over the
// intersection
Scalar pth = 0.0;
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
const auto& up = elemCtx.intensiveQuantities(upIdx, /*timeIdx=*/0);
if (up.mobility(phaseIdx) > 0.0) {
Scalar phaseVal = Toolbox::value(extQuants.pressureDifference(phaseIdx));
pth = std::max(pth, std::abs(phaseVal));
}
}
return pth;
}
// compute the defaults of the threshold pressures using the initial condition
void computeDefaultThresholdPressures_()
{
@@ -141,14 +107,13 @@ private:
auto elemIt = gridView.template begin</*codim=*/ 0>();
const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
ElementContext elemCtx(simulator_);
simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
for (; elemIt != elemEndIt; ++elemIt) {
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
elemCtx.updateStencil(elem);
elemCtx.updateAll(elem);
const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
for (unsigned scvfIdx = 0; scvfIdx < stencil.numInteriorFaces(); ++ scvfIdx) {
@@ -162,19 +127,30 @@ private:
unsigned equilRegionInside = this->elemEquilRegion_[insideElemIdx];
unsigned equilRegionOutside = this->elemEquilRegion_[outsideElemIdx];
if (equilRegionInside == equilRegionOutside)
// the current face is not at the boundary between EQUIL regions!
continue;
const auto& problem = elemCtx.problem();
// don't include connections with negligible flow
const Evaluation& trans = problem.transmissibility(elemCtx, i, j);
const Evaluation& trans = simulator_.problem().transmissibility(elemCtx, i, j);
Scalar faceArea = face.area();
if (std::abs(faceArea*getValue(trans)) < 1e-18)
continue;
double pth = calculateMaxDp(face, stencil, elemCtx, scvfIdx,
i, j,
insideElemIdx, outsideElemIdx);
// don't include connections with negligible flow
// determine the maximum difference of the pressure of any phase over the
// intersection
Scalar pth = 0.0;
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
const auto& up = elemCtx.intensiveQuantities(upIdx, /*timeIdx=*/0);
if (up.mobility(phaseIdx) > 0.0) {
Scalar phaseVal = Toolbox::value(extQuants.pressureDifference(phaseIdx));
pth = std::max(pth, std::abs(phaseVal));
}
}
int offset1 = equilRegionInside*this->numEquilRegions_ + equilRegionOutside;
int offset2 = equilRegionOutside*this->numEquilRegions_ + equilRegionInside;

33
ebos/ecltracermodel.hh
View File

@@ -86,9 +86,6 @@ class EclTracerModel : public EclGenericTracerModel<GetPropType<TypeTag, Propert
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
using Eval = DenseAd::Evaluation<double, numEq>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
public:
EclTracerModel(Simulator& simulator)
: BaseType(simulator.vanguard().gridView(),
@@ -228,17 +225,6 @@ protected:
}
void getVolumeFlux(unsigned& upIdx,
Scalar& v,
const FvBaseElementContext<TypeTag>& elemCtx,
const int tracerPhaseIdx,
unsigned scvfIdx)
{
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx*/ 0);
upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
}
// evaluate the flux(es) over one face
void computeFlux_(TracerEvaluation & freeFlux,
bool & isUpFree,
@@ -250,17 +236,17 @@ protected:
{
const auto& stencil = elemCtx.stencil(timeIdx);
const auto& scvf = stencil.interiorFace(scvfIdx);
unsigned inIdx = scvf.interiorIndex();
unsigned upIdx;
Scalar v;
getVolumeFlux(upIdx,
v,
elemCtx,
tracerPhaseIdx,
scvfIdx);
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned inIdx = extQuants.interiorIndex();
unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar A = scvf.area();
Scalar v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx));
if (inIdx == upIdx) {
@@ -411,8 +397,7 @@ protected:
auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timIdx*/ 0.0);
elemCtx.updateAll(*elemIt);
int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timIdx=*/0);
Scalar fVolume;
computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/0);