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

Modification to reduce PR to minimal changes.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen
2022-08-03 09:52:03 +02:00
parent 9de8d12276
commit 5fba14373b
9 changed files with 36 additions and 183 deletions
Download Patch File Download Diff File Expand all files Collapse all files

91
ebos/eclthresholdpressure.hh
View File

@@ -39,7 +39,6 @@
#include <algorithm>
#include <vector>
#include <opm/models/discretization/common/smallelementcontext.hh>
namespace Opm {
/*!
@@ -100,10 +99,16 @@ 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){
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);
@@ -114,7 +119,7 @@ private:
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));
@@ -122,72 +127,8 @@ private:
}
return pth;
}
template<class Face,class Stencil>
double calculateMaxDp(Face& face, Stencil& stencil,
SmallElementContext<TypeTag>& elemCtx,const unsigned& scvfIdx,
const unsigned& i,const unsigned& j,
const unsigned& insideElemIdx,const unsigned& outsideElemIdx){
typedef MathToolbox<Evaluation> Toolbox;
// 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);
Scalar Vin = elemCtx.dofVolume(i, /*timeIdx=*/0);
Scalar Vex = elemCtx.dofVolume(j, /*timeIdx=*/0);
Scalar thpres = 0.0;//NB ??problem.thresholdPressure(globalIndexIn, globalIndexEx);
// estimate the gravity correction: for performance reasons we use a simplified
// approach for this flux module that assumes that gravity is constant and always
// acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
const auto& problem = elemCtx.problem();
Scalar g = problem.gravity()[dimWorld - 1];
const auto& intQuantsIn = elemCtx.intensiveQuantities(i, /*timeIdx*/0);
const auto& intQuantsEx = elemCtx.intensiveQuantities(j, /*timeIdx*/0);
// 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"
// 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
// here...
Scalar zIn = problem.dofCenterDepth(elemCtx, i, /*timeIdx*/0);
Scalar zEx = problem.dofCenterDepth(elemCtx, j, /*timeIdx*/0);
// the distances from the DOF's depths. (i.e., the additional depth of the
// exterior DOF)
Scalar distZ = zIn - zEx;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
short dnIdx;
//
short upIdx;
Evaluation pressureDifference;
ExtensiveQuantities::calculatePhasePressureDiff_(upIdx,
dnIdx,
pressureDifference,
intQuantsIn,
intQuantsEx,
/*timeIdx*/0,//input
phaseIdx,//input
i,//input
j,//intput
Vin,
Vex,
insideElemIdx,
outsideElemIdx,
distZ*g,
thpres);
const IntensiveQuantities& up = (upIdx == i) ? intQuantsIn : intQuantsEx;
if (up.mobility(phaseIdx) > 0.0) {
Scalar phaseVal = Toolbox::value(pressureDifference);
pth = std::max(pth, std::abs(phaseVal));
}
}
return pth;
}
// compute the defaults of the threshold pressures using the initial condition
void computeDefaultThresholdPressures_()
{
@@ -208,13 +149,9 @@ private:
continue;
elemCtx.updateStencil(elem);
//
const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
for (unsigned scvfIdx = 0; scvfIdx < stencil.numInteriorFaces(); ++ scvfIdx) {
const auto& face = stencil.interiorFace(scvfIdx);
unsigned i = face.interiorIndex();
@@ -225,7 +162,6 @@ 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;
@@ -235,7 +171,6 @@ private:
Scalar faceArea = face.area();
if (std::abs(faceArea*getValue(trans)) < 1e-18)
continue;
double pth = calculateMaxDp(face, stencil, elemCtx, scvfIdx,
i, j,
insideElemIdx, outsideElemIdx);

31
ebos/ecltracermodel.hh
View File

@@ -227,39 +227,18 @@ protected:
freeVolume = phaseVolume * variable<LhsEval>(1.0, 0);
}
//template<class TypeTag>
void getVolumeFlux(unsigned& upIdx,
Scalar& v,
const FvBaseElementContext<TypeTag>& elemCtx,
const int tracerPhaseIdx,
unsigned scvfIdx
){
unsigned scvfIdx)
{
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx*/ 0);
upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
}
//template <class TypeTag>
void getVolumeFlux(unsigned& upIdx,
Scalar& v,
const SmallElementContext<TypeTag>& elemCtx,
const int tracerPhaseIdx,
unsigned scvfIdx
){
short upIdxV[numPhases];
short dnIdxV[numPhases];
Eval volumFlux[numPhases];
Eval pressureDifferences[numPhases];
ExtensiveQuantities::volumeAndPhasePressureDifferences(upIdxV,
dnIdxV,
volumFlux,
pressureDifferences,
elemCtx,
scvfIdx,
/*timeIdx*/ 0);
v = decay<Scalar>(volumFlux[tracerPhaseIdx]);
upIdx = upIdxV[tracerPhaseIdx] ;
}
// evaluate the flux(es) over one face
void computeFlux_(TracerEvaluation & freeFlux,
bool & isUpFree,
@@ -278,11 +257,10 @@ protected:
v,
elemCtx,
tracerPhaseIdx,
scvfIdx);
scvfIdx);
const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
const auto& fs = intQuants.fluidState();
Scalar A = scvf.area();
Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx));
if (inIdx == upIdx) {
@@ -433,7 +411,6 @@ protected:
auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
for (; elemIt != elemEndIt; ++ elemIt) {
//elemCtx.updateAll(*elemIt);
elemCtx.updatePrimaryStencil(*elemIt);
elemCtx.updatePrimaryIntensiveQuantities(/*timIdx*/ 0.0);
int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timIdx=*/0);

9
flow/flow_blackoil_tpfa.cpp
View File

@@ -20,10 +20,9 @@
#include <opm/material/common/ResetLocale.hpp>
#include <opm/simulators/flow/SimulatorFullyImplicitBlackoilEbos.hpp>
#include <opm/simulators/flow/Main.hpp>
// modifications from standard
#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
#include <opm/models/blackoil/blackoilintensivequantitiessimple.hh>
#include <opm/models/discretization/common/smallelementcontext.hh>
#include <opm/models/discretization/common/tpfalinearizer.hh>
namespace Opm {
@@ -45,12 +44,6 @@ namespace Opm {
template<class TypeTag>
struct LocalResidual<TypeTag, TTag::EclFlowProblemTPFA> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
// template<class TypeTag>
// struct ElementContext<TypeTag, TTag::EclFlowProblemTPFA> { using type = SmallElementContext<TypeTag>; };
// template<class TypeTag>
// struct IntensiveQuantities<TypeTag, TTag::EclFlowProblemTPFA> { using type = BlackOilIntensiveQuantitiesSimple<TypeTag>; };
}
}

23
opm/simulators/aquifers/AquiferNumerical.hpp
View File

@@ -252,31 +252,14 @@ private:
return sum_pressure_watervolume / sum_watervolume;
}
template<class ElemCtx>
const double getWaterFlux(ElemCtx& elem_ctx,unsigned face_idx) const{
template <class ElemCtx>
const double getWaterFlux(ElemCtx& elem_ctx, unsigned face_idx) const
{
const auto& exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
const double water_flux = Toolbox::value(exQuants.volumeFlux(this->phaseIdx_()));
return water_flux;
}
const double getWaterFlux(SmallElementContext<TypeTag>& elem_ctx,unsigned face_idx) const{
short upIdx[numPhases];
short dnIdx[numPhases];
Eval volumFlux[numPhases];
Eval pressureDifferences[numPhases];
ExtensiveQuantities::volumeAndPhasePressureDifferences(upIdx,
dnIdx,
volumFlux,
pressureDifferences,
elem_ctx,
face_idx,
/*timeIdx*/ 0);
return Toolbox::value(volumFlux[this->phaseIdx_()]);
}
double calculateAquiferFluxRate() const
{
double aquifer_flux = 0.0;

53
opm/simulators/flow/BlackoilModelEbos.hpp
View File

@@ -52,8 +52,6 @@
#include <opm/simulators/linalg/ISTLSolverEbos.hpp>
#include <opm/models/blackoil/blackoilintensivequantitiessimple.hh>
#include <dune/istl/owneroverlapcopy.hh>
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7)
#include <dune/common/parallel/communication.hh>
@@ -172,8 +170,6 @@ namespace Opm {
using Indices = GetPropType<TypeTag, Properties::Indices>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
typedef double Scalar;
static const int numEq = Indices::numEq;
@@ -572,56 +568,23 @@ namespace Opm {
ebosSolver.solve(x);
}
template<class IntensiveQuantity>
void updateIntensiveQuantity(const Problem& /*problem*/,
SolutionVector& solution,
const BVector& dx,
const IntensiveQuantity*)
{
ebosSimulator_.model().newtonMethod().update_(/*nextSolution=*/solution,
/*curSolution=*/solution,
/*update=*/dx,
/*resid=*/dx); // the update routines of the black
// oil model do not care about the
// residual
ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
}
void updateIntensiveQuantity(const Problem& problem,
SolutionVector& solution,
const BVector& dx,
const BlackOilIntensiveQuantitiesSimple<TypeTag>* /*intensive*/)
{
auto& model = ebosSimulator_.model();
auto& ebosNewtonMethod = model.newtonMethod();
ebosNewtonMethod.update_(/*nextSolution*/solution,
/*curSolution=*/solution,
/*update=*/dx,
/*resid=*/dx); // the update routines of the black
// oil model do not care about the
// residual
model.invalidateAndUpdateIntensiveQuantitiesSimple(problem,
solution,
/*timeIdx*/0);
}
/// Apply an update to the primary variables.
void updateSolution(const BVector& dx)
{
auto& problem = ebosSimulator_.problem();
auto& ebosNewtonMethod = ebosSimulator_.model().newtonMethod();
SolutionVector& solution = ebosSimulator_.model().solution(/*timeIdx=*/0);
ebosNewtonMethod.update_(/*nextSolution=*/solution,
/*curSolution=*/solution,
/*update=*/dx,
/*resid=*/dx); // the update routines of the black
// oil model do not care about the
// residual
// if the solution is updated, the intensive quantities need to be recalculated
//ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
IntensiveQuantities* dummy = NULL;//only for template spesialization
this->updateIntensiveQuantity(problem,solution, dx, dummy);
//ebosSimulator_.model().invalidateAndUpdateIntensiveQuantitiesSimple<IntensiveQuantities>(problem,
//solution,
// /*timeIdx=*/0);
ebosSimulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
}
/// Return true if output to cout is wanted.

4
opm/simulators/linalg/ISTLSolverEbos.hpp
View File

@@ -532,14 +532,14 @@ namespace Opm
return weights;
}
#if 0
Vector getTrueImpesWeights(int pressureVarIndex,SmallElementContext<TypeTag>& /*elemCtx*/) const
{
Vector weights(rhs_->size());
Amg::getTrueImpesWeights<Evaluation>(pressureVarIndex, weights, simulator_.model());
return weights;
}
#endif
/// Zero out off-diagonal blocks on rows corresponding to overlap cells
/// Diagonal blocks on ovelap rows are set to diag(1.0).

4
opm/simulators/wells/BlackoilWellModel.hpp
View File

@@ -273,8 +273,8 @@ namespace Opm {
void addWellContributions(SparseMatrixAdapter& jacobian) const;
void addReservoirSourceTerms(GlobalEqVector& residual,
SparseMatrixAdapter& jacobian) const;
// void addReservoirSourceTerms(GlobalEqVector& residual,
// SparseMatrixAdapter& jacobian) const;
// called at the beginning of a report step
void beginReportStep(const int time_step);

2
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@@ -1202,6 +1202,7 @@ namespace Opm {
}
}
#if 0
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
@@ -1231,6 +1232,7 @@ namespace Opm {
}
}
}
#endif
template<typename TypeTag>
void

2
opm/simulators/wells/WellInterface.hpp
View File

@@ -288,7 +288,7 @@ public:
const GroupState& group_state,
DeferredLogger& deferred_logger);
const std::vector<RateVector>& connectionRates() const {return connectionRates_;}
// const std::vector<RateVector>& connectionRates() const {return connectionRates_;}
protected:
// simulation parameters