Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Extend support for KRNUM to solvent module #751

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Extend support for KRNUM to solvent module #751

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 23 additions & 12 deletions opm/models/blackoil/blackoilintensivequantities.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,11 +39,11 @@
#include "blackoilmicpmodules.hh"
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/models/common/directionalmobility.hh>

#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/utility/CopyablePtr.hpp>

#include <dune/common/fmatrix.hh>

#include <cstring>
Expand Down Expand Up @@ -421,22 +421,15 @@ public:
}
paramCache.updateAll(fluidState_);

auto mobilities = getMobilityVector();
// compute the phase densities and transform the phase permeabilities into mobilities
int nmobilities = 1;
std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
if (dirMob_) {
for (int i=0; i<3; i++) {
nmobilities += 1;
mobilities.push_back(&(dirMob_->getArray(i)));
}
}
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx))
continue;
const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
fluidState_.setInvB(phaseIdx, b);
const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
for (int i = 0; i<nmobilities; i++) {
for (std::size_t i = 0; i<mobilities.size(); i++) {
if (enableExtbo && phaseIdx == oilPhaseIdx) {
(*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
}
Expand Down Expand Up @@ -563,6 +556,20 @@ public:
*/
const FluidState& fluidState() const
{ return fluidState_; }

/*!
* Returns a vector of pointer to the mobilities. This includes the directional mobilities if they have
* been enabled.
*/
std::vector<std::array<Evaluation,numPhases>*> getMobilityVector() {
std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
if (dirMob_) {
for (int i=0; i<3; i++) {
mobilities.push_back(&(dirMob_->getArray(i)));
}
}
return mobilities;
}

/*!
* \copydoc ImmiscibleIntensiveQuantities::mobility
Expand Down Expand Up @@ -658,8 +665,12 @@ private:

// Instead of writing a custom copy constructor and a custom assignment operator just to handle
// the dirMob_ unique ptr member variable when copying BlackOilIntensiveQuantites (see for example
// updateIntensitiveQuantities_() in fvbaseelementcontext.hh for a copy example) we write the below
// custom wrapper class CopyablePtr which wraps the unique ptr and makes it copyable.
// updateIntensitiveQuantities_() in fvbaseelementcontext.hh for a copy example) we use the
// wrapper class
//
// DirectionalMobilityPtr (aka Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>)
//
// which wraps the unique ptr and makes it copyable.
//
// The advantage of this approach is that we avoid having to call all the base class copy constructors and
// assignment operators explicitly (which is needed when writing the custom copy constructor and assignment
Expand Down
141 changes: 110 additions & 31 deletions opm/models/blackoil/blackoilsolventmodules.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
#include <opm/models/common/quantitycallbacks.hh>

#include <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
#include <opm/utility/CopyablePtr.hpp>

#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
Expand All @@ -48,6 +49,8 @@
#include <opm/input/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
#endif

#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>

#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/Exceptions.hpp>

Expand Down Expand Up @@ -90,7 +93,6 @@ class BlackOilSolventModule
static constexpr unsigned numPhases = FluidSystem::numPhases;
static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();


public:
#if HAVE_ECL_INPUT
/*!
Expand Down Expand Up @@ -669,6 +671,27 @@ class BlackOilSolventIntensiveQuantities
static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
static constexpr double cutOff = 1e-12;

struct SolventDirectionalMobility {
SolventDirectionalMobility(const SolventDirectionalMobility& other)
: mobX_{other.mobX_}, mobY_{other.mobY_}, mobZ_{other.mobZ_} {}
SolventDirectionalMobility(const Evaluation& mX, const Evaluation& mY, const Evaluation& mZ)
: mobX_{mX}, mobY_{mY}, mobZ_{mZ} {}
SolventDirectionalMobility() : mobX_{0.0}, mobY_{0.0}, mobZ_{0.0} {}

Evaluation& getMob(int idx) {
switch(idx) {
case 0:
return mobX_;
case 1:
return mobY_;
case 2:
return mobZ_;
default:
throw std::runtime_error("Unexpected mobility array index");
}
}
Evaluation mobX_, mobY_, mobZ_;
};

public:
/*!
Expand Down Expand Up @@ -711,8 +734,16 @@ public:
auto& fs = asImp_().fluidState_;
fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_);

solventMobility_ = 0.0;

const auto& problem = elemCtx.problem();
auto mobilities = asImp_().getMobilityVector();
//if (problem.materialLawManager()->hasDirectionalRelperms()) {
if (mobilities.size() > 1) {
solventDirectionalMobilities_ = std::make_unique<SolventDirectionalMobility>();
}
auto solventMobilities = getSolventMobilities();
for (int i=0; i<solventMobilities.size(); i++) {
*solventMobilities[i] = 0.0;
}
// apply a cut-off. Don't waste calculations if no solvent
if (solventSaturation().value() < cutOff)
return;
Expand All @@ -724,7 +755,6 @@ public:
const Evaluation& pgImisc = fs.pressure(gasPhaseIdx);

// compute capillary pressure for miscible fluid
const auto& problem = elemCtx.problem();
const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
Evaluation pgMisc = 0.0;
std::array<Evaluation, numPhases> pC;
Expand Down Expand Up @@ -790,26 +820,28 @@ public:
const Evaluation mkrgt = msfnKrsg.eval(F_totalGas, /*extrapolate=*/true) * sof2Krn.eval(oilGasSolventSat, /*extrapolate=*/true);
const Evaluation mkro = msfnKro.eval(F_totalGas, /*extrapolate=*/true) * sof2Krn.eval(oilGasSolventSat, /*extrapolate=*/true);

Evaluation& kro = asImp_().mobility_[oilPhaseIdx];
Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
for (std::size_t i = 0; i<mobilities.size(); i++) {
Evaluation& kro = (*mobilities[i])[oilPhaseIdx];
Evaluation& krg = (*mobilities[i])[gasPhaseIdx];

// combine immiscible and miscible part of the relperm
krg *= (1.0 - miscibility);
krg += miscibility * mkrgt;
kro *= (1.0 - miscibility);
kro += miscibility * mkro;
// combine immiscible and miscible part of the relperm
krg *= (1.0 - miscibility);
krg += miscibility * mkrgt;
kro *= (1.0 - miscibility);
kro += miscibility * mkro;
}
}



// compute the mobility of the solvent "phase" and modify the gas phase
const auto& ssfnKrg = SolventModule::ssfnKrg(elemCtx, dofIdx, timeIdx);
const auto& ssfnKrs = SolventModule::ssfnKrs(elemCtx, dofIdx, timeIdx);

Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
solventMobility_ = krg * ssfnKrs.eval(Fsolgas, /*extrapolate=*/true);
krg *= ssfnKrg.eval(Fhydgas, /*extrapolate=*/true);

for (std::size_t i = 0; i<mobilities.size(); i++) {
Evaluation& krg = (*mobilities[i])[gasPhaseIdx];
*solventMobilities[i] = krg * ssfnKrs.eval(Fsolgas, /*extrapolate=*/true);
krg *= ssfnKrg.eval(Fhydgas, /*extrapolate=*/true);
}
}

/*!
Expand Down Expand Up @@ -837,9 +869,10 @@ public:

effectiveProperties(elemCtx, scvIdx, timeIdx);

solventMobility_ /= solventViscosity_;


auto solventMobilities = getSolventMobilities();
for (std::size_t i=0; i<solventMobilities.size(); i++) {
*solventMobilities[i] /= solventViscosity_;
}
}

const Evaluation& solventSaturation() const
Expand All @@ -851,9 +884,40 @@ public:
const Evaluation& solventViscosity() const
{ return solventViscosity_; }

std::vector<Evaluation *> getSolventMobilities()
{
std::vector<Evaluation *> solventMobilities = {&solventMobility_};
if (solventDirectionalMobilities_) {
for (int i=0; i<3; i++) {
solventMobilities.push_back(&(solventDirectionalMobilities_->getMob(i)));
}
}
return solventMobilities;
}

const Evaluation& solventMobility() const
{ return solventMobility_; }

const Evaluation& solventMobility(FaceDir::DirEnum facedir) const
{
using Dir = FaceDir::DirEnum;
if (solventDirectionalMobilities_) {
switch(facedir) {
case Dir::XPlus:
return solventDirectionalMobilities_->mobX_;
case Dir::YPlus:
return solventDirectionalMobilities_->mobY_;
case Dir::ZPlus:
return solventDirectionalMobilities_->mobZ_;
default:
throw std::runtime_error("Unexpected face direction");
}
}
else {
return solventMobility_;
}

}
const Evaluation& solventInverseFormationVolumeFactor() const
{ return solventInvFormationVolumeFactor_; }

Expand Down Expand Up @@ -1039,6 +1103,7 @@ protected:
Evaluation solventDensity_;
Evaluation solventViscosity_;
Evaluation solventMobility_;
Utility::CopyablePtr<SolventDirectionalMobility> solventDirectionalMobilities_;
Evaluation solventInvFormationVolumeFactor_;

Scalar solventRefDensity_;
Expand Down Expand Up @@ -1129,6 +1194,7 @@ public:
const auto& gradCalc = elemCtx.gradientCalculator();
PressureCallback<TypeTag> pressureCallback(elemCtx);

const auto& problem = elemCtx.problem();
const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
const auto& faceNormal = scvf.normal();

Expand All @@ -1148,8 +1214,8 @@ public:
if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) {
// estimate the gravitational acceleration at a given SCV face
// using the arithmetic mean
const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
const auto& gIn = problem.gravity(elemCtx, i, timeIdx);
const auto& gEx = problem.gravity(elemCtx, j, timeIdx);

const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
Expand Down Expand Up @@ -1208,15 +1274,20 @@ public:
}

const auto& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
const auto& materialLawManager = problem.materialLawManager();
FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
if (materialLawManager->hasDirectionalRelperms()) {
facedir = scvf.faceDirFromDirId(); // direction (X, Y, or Z) of the face
}

// this is also slightly hacky because it assumes that the derivative of the
// flux between two DOFs only depends on the primary variables in the
// upstream direction. For non-TPFA flux approximation schemes, this is not
// true...
if (solventUpstreamDofIdx_ == i)
solventVolumeFlux_ = solventPGradNormal*up.solventMobility();
solventVolumeFlux_ = solventPGradNormal*up.solventMobility(facedir);
else
solventVolumeFlux_ = solventPGradNormal*scalarValue(up.solventMobility());
solventVolumeFlux_ = solventPGradNormal*scalarValue(up.solventMobility(facedir));
}

/*!
Expand All @@ -1230,6 +1301,7 @@ public:
unsigned timeIdx)
{
const ExtensiveQuantities& extQuants = asImp_();
const auto& problem = elemCtx.problem();

unsigned interiorDofIdx = extQuants.interiorIndex();
unsigned exteriorDofIdx = extQuants.exteriorIndex();
Expand All @@ -1241,12 +1313,12 @@ public:
unsigned I = elemCtx.globalSpaceIndex(interiorDofIdx, timeIdx);
unsigned J = elemCtx.globalSpaceIndex(exteriorDofIdx, timeIdx);

Scalar thpres = elemCtx.problem().thresholdPressure(I, J);
Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
Scalar thpres = problem.thresholdPressure(I, J);
Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
Scalar g = problem.gravity()[dimWorld - 1];

Scalar zIn = elemCtx.problem().dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
Scalar zEx = elemCtx.problem().dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
Scalar distZ = zIn - zEx;

const Evaluation& rhoIn = intQuantsIn.solventDensity();
Expand Down Expand Up @@ -1285,16 +1357,23 @@ public:
return;
}

Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
Scalar faceArea = scvf.area();
const auto& materialLawManager = problem.materialLawManager();
FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
if (materialLawManager->hasDirectionalRelperms()) {
facedir = scvf.faceDirFromDirId(); // direction (X, Y, or Z) of the face
}

const IntensiveQuantities& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
if (solventUpstreamDofIdx_ == interiorDofIdx)
solventVolumeFlux_ =
up.solventMobility()
up.solventMobility(facedir)
*(-trans/faceArea)
*pressureDiffSolvent;
else
solventVolumeFlux_ =
scalarValue(up.solventMobility())
scalarValue(up.solventMobility(facedir))
*(-trans/faceArea)
*pressureDiffSolvent;
}
Expand Down