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 add tpfa diffusion

Browse Source 
fixed tpfa for diffusion module

fixed type error
This commit is contained in:
hnil
2023-08-23 11:50:52 +02:00
committed by Kai Bao
parent 7ec77ada7e
commit 7acecc7e57
3 changed files with 84 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

78
opm/models/blackoil/blackoildiffusionmodule.hh
View File

@@ -47,6 +47,13 @@ namespace Opm {
template <class TypeTag, bool enableDiffusion>
class BlackOilDiffusionModule;
template <class TypeTag, bool enableDiffusion>
class BlackOilDiffusionExtensiveQuantities;
template <class TypeTag, bool enableDiffusion>
class NumericalConstants;
/*!
* \copydoc Opm::BlackOilDiffusionModule
*/
@@ -95,6 +102,7 @@ class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/true>
using Toolbox = MathToolbox<Evaluation>;
public:
using ExtensiveQuantities = BlackOilDiffusionExtensiveQuantities<TypeTag,true>;
/*!
* \brief Register all run-time parameters for the diffusion module.
*/
@@ -112,10 +120,21 @@ public:
// Only work if diffusion is enabled run-time by DIFFUSE in the deck
if(!FluidSystem::enableDiffusion())
return;
const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
const auto& diffusivity = extQuants.diffusivity();
const auto& effectiveDiffusionCoefficient = extQuants.effectiveDiffusionCoefficient();
addDiffusiveFlux(flux, fluidStateI, fluidStateJ, diffusivity, effectiveDiffusionCoefficient);
}
template<class FluidState,class EvaluationArray>
static void addDiffusiveFlux(RateVector& flux,
const FluidState& fluidStateI,
const FluidState& fluidStateJ,
const Evaluation& diffusivity,
const EvaluationArray& effectiveDiffusionCoefficient)
{
unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
@@ -127,7 +146,7 @@ public:
continue;
}
// no diffusion in gas phase in water + gas system.
// no diffusion in gas phase in water + gas system.
if (FluidSystem::gasPhaseIdx == phaseIdx && !FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
continue;
}
@@ -165,8 +184,9 @@ public:
- bSAvg
* convFactor
* diffR
* extQuants.diffusivity()
* extQuants.effectiveDiffusionCoefficient(phaseIdx, solventCompIdx);
* diffusivity
* effectiveDiffusionCoefficient[phaseIdx][solventCompIdx];
// mass flux of solute component (gas in oil or oil in gas)
unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx);
@@ -174,8 +194,8 @@ public:
bSAvg
* diffR
* convFactor
* extQuants.diffusivity()
* extQuants.effectiveDiffusionCoefficient(phaseIdx, soluteCompIdx);
* diffusivity
* effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx];
}
}
@@ -271,7 +291,7 @@ class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
enum { numPhases = FluidSystem::numPhases };
enum { numComponents = FluidSystem::numComponents };
@@ -339,9 +359,17 @@ protected:
// Only work if diffusion is enabled run-time by DIFFUSE in the deck
if(!FluidSystem::enableDiffusion())
return;
using Toolbox = MathToolbox<Evaluation>;
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
update_(fluidState, paramCache, intQuants, timeIdx);
}
template<class FluidState>
void update_(FluidState& fluidState,
typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
const IntensiveQuantities& intQuants,
unsigned timeIdx){
using Toolbox = MathToolbox<Evaluation>;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@@ -353,8 +381,10 @@ protected:
}
// Based on Millington, R. J., & Quirk, J. P. (1961).
//constexpr auto& numconst = GetPropValue<TypeTag, Properties::NumericalConstants>;
constexpr double myeps = 0.0001; //numconst.blackoildiffusionmoduleeps;
const Evaluation& base =
Toolbox::max(0.0001,
Toolbox::max(myeps, //0.0001,
intQuants.porosity()
* intQuants.fluidState().saturation(phaseIdx));
tortuosity_[phaseIdx] =
@@ -436,6 +466,7 @@ public:
throw std::logic_error("The method effectiveDiffusionCoefficient() "
"does not make sense if diffusion is disabled.");
}
};
/*!
@@ -450,6 +481,7 @@ class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
using GridView = GetPropType<TypeTag, Properties::GridView>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Toolbox = MathToolbox<Evaluation>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
enum { dimWorld = GridView::dimensionworld };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
@@ -457,7 +489,8 @@ class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
using DimVector = Dune::FieldVector<Scalar, dimWorld>;
using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
public:
using EvaluationArray = Evaluation[numPhases][numComponents];
protected:
/*!
* \brief Update the quantities required to calculate
@@ -477,10 +510,17 @@ protected:
const Scalar diffusivity = elemCtx.problem().diffusivity(elemCtx, face.interiorIndex(), face.exteriorIndex());
const Scalar faceArea = face.area();
diffusivity_ = diffusivity / faceArea; //opm-models expects pr area flux
diffusivity_ = diffusivity / faceArea;
update_(diffusivity_, effectiveDiffusionCoefficient_, intQuantsInside, intQuantsOutside);
Valgrind::CheckDefined(diffusivity_);
}
public:
static void update_(Scalar& diffusivity,
EvaluationArray& effectiveDiffusionCoefficient,
const IntensiveQuantities& intQuantsInside,
const IntensiveQuantities& intQuantsOutside){
//opm-models expects pr area flux
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@@ -492,17 +532,17 @@ protected:
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
// use the arithmetic average for the effective
// diffusion coefficients.
effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
effectiveDiffusionCoefficient[phaseIdx][compIdx] =
(intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
+
intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
/ 2;
Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
//Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
}
}
}
protected:
template <class Context, class FluidState>
void updateBoundary_(const Context&,
unsigned,
@@ -532,9 +572,13 @@ public:
const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
{ return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
const auto& effectiveDiffusionCoefficient() const{
return effectiveDiffusionCoefficient_;
}
private:
Scalar diffusivity_;
Evaluation effectiveDiffusionCoefficient_[numPhases][numComponents];
EvaluationArray effectiveDiffusionCoefficient_;
};
} // namespace Opm

20
opm/models/blackoil/blackoillocalresidualtpfa.hh
View File

@@ -118,6 +118,7 @@ public:
double Vex;
double inAlpha;
double outAlpha;
double diffusivity;
};
/*!
* \copydoc FvBaseLocalResidual::computeStorage
@@ -277,7 +278,7 @@ public:
// 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.)
Scalar g = problem.gravity()[dimWorld - 1];
const Scalar g = problem.gravity()[dimWorld - 1];
const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
@@ -294,8 +295,9 @@ public:
// for thermal harmonic mean of half trans
const Scalar inAlpha = problem.thermalHalfTransmissibility(globalIndexIn, globalIndexEx);
const Scalar outAlpha = problem.thermalHalfTransmissibility(globalIndexEx, globalIndexIn);
const Scalar diffusivity = problem.diffusivity(globalIndexEx, globalIndexIn);
const ResidualNBInfo res_nbinfo {trans, faceArea, thpres, distZ * g, facedir, Vin, Vex, inAlpha, outAlpha};
const ResidualNBInfo res_nbinfo {trans, faceArea, thpres, distZ * g, facedir, Vin, Vex, inAlpha, outAlpha, diffusivity};
calculateFluxes_(flux,
darcy,
@@ -324,6 +326,7 @@ public:
const FaceDir::DirEnum facedir = nbInfo.faceDirection;
const Scalar inAlpha = nbInfo.inAlpha;
const Scalar outAlpha = nbInfo.outAlpha;
const Scalar diffusivity = nbInfo.diffusivity;
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
@@ -440,9 +443,18 @@ public:
// BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
// deal with diffusion (if present)
static_assert(!enableDiffusion, "Relevant computeFlux() method must be implemented for this module before enabling.");
// DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
//static_assert(!enableDiffusion, "Relevant computeFlux() method must be implemented for this module before enabling.");
if constexpr(enableDiffusion){
typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
Scalar tmpdiffusivity = diffusivity/faceArea;
DiffusionModule::ExtensiveQuantities::update_(tmpdiffusivity,effectiveDiffusionCoefficient,intQuantsIn,intQuantsEx);
DiffusionModule::addDiffusiveFlux(flux,
intQuantsIn.fluidState(),
intQuantsEx.fluidState(),
tmpdiffusivity,
effectiveDiffusionCoefficient);
}
// deal with micp (if present)
static_assert(!enableMICP, "Relevant computeFlux() method must be implemented for this module before enabling.");
// MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);

8
opm/models/discretization/common/tpfalinearizer.hh
View File

@@ -112,6 +112,7 @@ class TpfaLinearizer
static const bool linearizeNonLocalElements = getPropValue<TypeTag, Properties::LinearizeNonLocalElements>();
static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
static const bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
// copying the linearizer is not a good idea
TpfaLinearizer(const TpfaLinearizer&);
//! \endcond
@@ -449,14 +450,19 @@ private:
Scalar inAlpha {0.};
Scalar outAlpha {0.};
FaceDirection dirId = FaceDirection::Unknown;
Scalar diffusivity {0.};
if constexpr(enableEnergy){
inAlpha = problem_().thermalHalfTransmissibility(myIdx, neighborIdx);
outAlpha = problem_().thermalHalfTransmissibility(neighborIdx, myIdx);
}
if constexpr(enableDiffusion){
diffusivity = problem_().diffusivity(myIdx, neighborIdx);
}
if (materialLawManager->hasDirectionalRelperms()) {
dirId = scvf.faceDirFromDirId();
}
loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, {trans, area, thpres, dZg, dirId, Vin, Vex, inAlpha, outAlpha}, nullptr};
loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, {trans, area, thpres, dZg, dirId, Vin, Vex, inAlpha, outAlpha, diffusivity}, nullptr};
}
}
neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());