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.

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

Mixing energy #309

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
8 changes: 7 additions & 1 deletion opm/material/fluidstates/BlackOilFluidState.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -321,7 +321,13 @@ class BlackOilFluidState
* exception!
*/
Scalar internalEnergy(unsigned phaseIdx OPM_UNUSED) const
{ return (*enthalpy_)[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); }
{
//double entFac=0.0;
double entFac=1.0;
// this may be setting which for entFac =1 make this probably E100 compatible
// however this has to be consitent with call to enthalpy
return (*enthalpy_)[phaseIdx] - entFac*pressure(phaseIdx)/density(phaseIdx);
}

//////
// slow methods
Expand Down
102 changes: 75 additions & 27 deletions opm/material/fluidsystems/BlackOilFluidSystem.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -551,7 +551,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
if (enableDissolvedGas()) {
// miscible oil
const LhsEval& Rs = Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
const LhsEval& bo = inverseFormationVolumeFactor<FluidState,LhsEval>(fluidState, phaseIdx, regionIdx);

return
bo*referenceDensity(oilPhaseIdx, regionIdx)
Expand All @@ -560,7 +560,7 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S

// immiscible oil
const LhsEval Rs(0.0);
const auto& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
const auto& bo = inverseFormationVolumeFactor<FluidState,LhsEval>(fluidState, phaseIdx, regionIdx);

return referenceDensity(phaseIdx, regionIdx)*bo;
}
Expand All @@ -569,16 +569,15 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
if (enableVaporizedOil()) {
// miscible gas
const LhsEval& Rv = Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);

const LhsEval& bg = inverseFormationVolumeFactor<FluidState,LhsEval>(fluidState, phaseIdx, regionIdx);
return
bg*referenceDensity(gasPhaseIdx, regionIdx)
+ Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
}

// immiscible gas
const LhsEval Rv(0.0);
const auto& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
const auto& bg = inverseFormationVolumeFactor<FluidState,LhsEval>(fluidState, phaseIdx, regionIdx);
return bg*referenceDensity(phaseIdx, regionIdx);
}

Expand Down Expand Up @@ -655,6 +654,74 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
}

template <class FluidState, class LhsEval = typename FluidState::Scalar>
static LhsEval energy(const FluidState& fluidState,
unsigned phaseIdx,
unsigned regionIdx)
{
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that function should be named either internalEnergy() or enthalpy(). Also, internalEnergy() and enthalpy() need to be consistent, i.e. h_alpha = u_alpha + p_alpha/rho_alpha.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I followed your example and this energy is only used to calculate enthalpy and could probably be made private.

Copy link
Contributor

@andlaus andlaus Oct 18, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah, okay. maybe this can even replace the enthalpy() and/or internalEnergy() functions directly?

assert(0 <= phaseIdx && phaseIdx <= numPhases);
assert(0 <= regionIdx && regionIdx <= numRegions());

//const auto& p = fluidState.pressure(phaseIdx);
//const auto& T = fluidState.temperature(phaseIdx);
const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));
switch (phaseIdx) {
case oilPhaseIdx: {
if (enableDissolvedGas()) {
// miscible oil
const LhsEval& Rs = Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
//const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
const LhsEval& bo = inverseFormationVolumeFactor<FluidState,LhsEval>(fluidState, phaseIdx, regionIdx);
//const LhsEval& bo = Opm::decay<LhsEval>(fluidState.invB(phaseIdx));
const LhsEval& Rv = Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
LhsEval energy =
bo*referenceDensity(oilPhaseIdx, regionIdx)
*oilPvt_->internalEnergy(regionIdx, T, p, Rs)
+ Rs*bo*referenceDensity(gasPhaseIdx, regionIdx)
*gasPvt_->internalEnergy(regionIdx, T, p, Rv);//this a to simple model
return energy/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
}

// immiscible oil
const LhsEval Rs(0.0);
//const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
//const LhsEval& bo = Opm::decay<LhsEval>(fluidState.invB(phaseIdx));
return oilPvt_->internalEnergy(regionIdx, T, p, Rs);
}

case gasPhaseIdx: {
if (enableVaporizedOil()) {
// miscible gas
const LhsEval& Rs = Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& Rv = Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
const LhsEval& bg = inverseFormationVolumeFactor<FluidState,LhsEval>(fluidState, phaseIdx, regionIdx);
//const LhsEval& bg = Opm::decay<LhsEval>(fluidState.invB(phaseIdx));
LhsEval energy =
bg*referenceDensity(gasPhaseIdx, regionIdx)
*gasPvt_->internalEnergy(regionIdx, T, p, Rv)
+ Rv*bg*referenceDensity(oilPhaseIdx, regionIdx)
*oilPvt_->internalEnergy(regionIdx, T, p, Rs);// this is a to simple model
return energy/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
}

// immiscible gas
const LhsEval Rv(0.0);
// const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);

return gasPvt_->internalEnergy(regionIdx, T, p, Rv);

}

case waterPhaseIdx:
return waterPvt_->internalEnergy(regionIdx, T, p);

}

throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
}


/*!
* \brief Returns the formation volume factor \f$B_\alpha\f$ of an "undersaturated"
* fluid phase
Expand Down Expand Up @@ -963,28 +1030,9 @@ class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<S
assert(0 <= regionIdx && regionIdx <= numRegions());

const auto& p = Opm::decay<LhsEval>(fluidState.pressure(phaseIdx));
const auto& T = Opm::decay<LhsEval>(fluidState.temperature(phaseIdx));

switch (phaseIdx) {
case oilPhaseIdx:
return
oilPvt_->internalEnergy(regionIdx, T, p, Opm::BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);

case gasPhaseIdx:
return
gasPvt_->internalEnergy(regionIdx, T, p, Opm::BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);

case waterPhaseIdx:
return
waterPvt_->internalEnergy(regionIdx, T, p)
+ p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);

default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
}

throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
double entFac = 1.0;
return energy<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx)
+ entFac*p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
}

/*!
Expand Down