BlackoilEnergyModule: only instance code when fully implicit

use the dummy overload otherwise
This commit is contained in:
Arne Morten Kvarving committed 2026-06-05 15:16:38 +02:00
1 parent f6c253a3f9
commit 842fa44a47
6 files changed
+134 -135

No files matched your search

+65 -93
View File
@@ -54,13 +54,16 @@
namespace Opm {
template <class TypeTag, EnergyModules activeModule>
class BlackOilEnergyModule;
/*!
* \ingroup BlackOil
* \brief Contains the high level supplements required to extend the black oil
* model by energy.
*/
template <class TypeTag, EnergyModules activeModule = getPropValue<TypeTag, Properties::EnergyModuleType>()>
class BlackOilEnergyModule
template <class TypeTag>
class BlackOilEnergyModule<TypeTag, EnergyModules::FullyImplicitThermal>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
@@ -77,7 +80,7 @@ class BlackOilEnergyModule
static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;
static constexpr unsigned enableFullyImplicitThermal = (activeModule == EnergyModules::FullyImplicitThermal);
static constexpr unsigned enableFullyImplicitThermal = true;
static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
static constexpr unsigned numPhases = FluidSystem::numPhases;
@@ -89,9 +92,7 @@ public:
*/
static void registerParameters()
{
if constexpr (enableFullyImplicitThermal) {
VtkBlackOilEnergyModule<TypeTag>::registerParameters();
}
VtkBlackOilEnergyModule<TypeTag>::registerParameters();
}
/*!
@@ -100,19 +101,12 @@ public:
static void registerOutputModules(Model& model,
Simulator& simulator)
{
if constexpr (enableFullyImplicitThermal) {
model.addOutputModule(std::make_unique<VtkBlackOilEnergyModule<TypeTag>>(simulator));
}
model.addOutputModule(std::make_unique<VtkBlackOilEnergyModule<TypeTag>>(simulator));
}
static OPM_HOST_DEVICE bool primaryVarApplies(unsigned pvIdx)
{
if constexpr (enableFullyImplicitThermal) {
return pvIdx == temperatureIdx;
}
else {
return false;
}
return pvIdx == temperatureIdx;
}
static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
@@ -132,12 +126,7 @@ public:
static OPM_HOST_DEVICE bool eqApplies(unsigned eqIdx)
{
if constexpr (enableFullyImplicitThermal) {
return eqIdx == contiEnergyEqIdx;
}
else {
return false;
}
return eqIdx == contiEnergyEqIdx;
}
static std::string eqName([[maybe_unused]] unsigned eqIdx)
@@ -160,73 +149,66 @@ public:
const IntensiveQuantities& intQuants)
{
using LhsEval = typename StorageType::value_type;
const FluidSystem& fsys = intQuants.getFluidSystem();
if constexpr (enableFullyImplicitThermal) {
const FluidSystem& fsys = intQuants.getFluidSystem();
const auto& poro = decay<LhsEval>(intQuants.porosity());
const auto& poro = decay<LhsEval>(intQuants.porosity());
// accumulate the internal energy of the fluids
const auto& fs = intQuants.fluidState();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!fsys.phaseIsActive(phaseIdx)) {
continue;
}
const auto& u = decay<LhsEval>(fs.internalEnergy(phaseIdx));
const auto& S = decay<LhsEval>(fs.saturation(phaseIdx));
const auto& rho = decay<LhsEval>(fs.density(phaseIdx));
storage[contiEnergyEqIdx] += poro*S*u*rho;
// accumulate the internal energy of the fluids
const auto& fs = intQuants.fluidState();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!fsys.phaseIsActive(phaseIdx)) {
continue;
}
// add the internal energy of the rock
const Scalar rockFraction = intQuants.rockFraction();
const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
storage[contiEnergyEqIdx] += rockFraction * uRock;
storage[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
const auto& u = decay<LhsEval>(fs.internalEnergy(phaseIdx));
const auto& S = decay<LhsEval>(fs.saturation(phaseIdx));
const auto& rho = decay<LhsEval>(fs.density(phaseIdx));
storage[contiEnergyEqIdx] += poro*S*u*rho;
}
// add the internal energy of the rock
const Scalar rockFraction = intQuants.rockFraction();
const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
storage[contiEnergyEqIdx] += rockFraction * uRock;
storage[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
}
OPM_HOST_DEVICE static void computeFlux([[maybe_unused]] RateVector& flux,
[[maybe_unused]] const ElementContext& elemCtx,
[[maybe_unused]] unsigned scvfIdx,
[[maybe_unused]] unsigned timeIdx)
OPM_HOST_DEVICE static void computeFlux(RateVector& flux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx)
{
if constexpr (enableFullyImplicitThermal) {
flux[contiEnergyEqIdx] = 0.0;
flux[contiEnergyEqIdx] = 0.0;
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
const unsigned focusIdx = elemCtx.focusDofIndex();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
if (upIdx == focusIdx) {
addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
}
else {
addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
}
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
const unsigned focusIdx = elemCtx.focusDofIndex();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
// diffusive energy flux
flux[contiEnergyEqIdx] += extQuants.energyFlux();
flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
if (upIdx == focusIdx) {
addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
}
else {
addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
}
}
// diffusive energy flux
flux[contiEnergyEqIdx] += extQuants.energyFlux();
flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
}
template<class RateVectorT>
OPM_HOST_DEVICE static void addHeatFlux(RateVectorT& flux,
const Evaluation& heatFlux)
{
if constexpr (enableFullyImplicitThermal) {
// diffusive energy flux
flux[contiEnergyEqIdx] += heatFlux;
flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
}
// diffusive energy flux
flux[contiEnergyEqIdx] += heatFlux;
flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
}
template <class UpEval, class RateVectorT, class Eval, class FluidState>
@@ -262,9 +244,7 @@ public:
OPM_HOST_DEVICE static void addToEnthalpyRate(RateVector& flux,
const Evaluation& hRate)
{
if constexpr (enableFullyImplicitThermal) {
flux[contiEnergyEqIdx] += hRate;
}
flux[contiEnergyEqIdx] += hRate;
}
/*!
@@ -274,9 +254,7 @@ public:
OPM_HOST_DEVICE static void assignPrimaryVars(PrimaryVariables& priVars,
const FluidState& fluidState)
{
if constexpr (enableFullyImplicitThermal) {
priVars[temperatureIdx] = getValue(fluidState.temperature(/*phaseIdx=*/0));
}
priVars[temperatureIdx] = getValue(fluidState.temperature(/*phaseIdx=*/0));
}
/*!
@@ -286,10 +264,8 @@ public:
const PrimaryVariables& oldPv,
const EqVector& delta)
{
if constexpr (enableFullyImplicitThermal) {
// do a plain unchopped Newton update
newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
}
// do a plain unchopped Newton update
newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
}
/*!
@@ -316,26 +292,22 @@ public:
template <class DofEntity>
static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
{
if constexpr (enableFullyImplicitThermal) {
const unsigned dofIdx = model.dofMapper().index(dof);
const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
outstream << priVars[temperatureIdx];
}
const unsigned dofIdx = model.dofMapper().index(dof);
const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
outstream << priVars[temperatureIdx];
}
template <class DofEntity>
static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
{
if constexpr (enableFullyImplicitThermal) {
const unsigned dofIdx = model.dofMapper().index(dof);
PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
const unsigned dofIdx = model.dofMapper().index(dof);
PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
instream >> priVars0[temperatureIdx];
instream >> priVars0[temperatureIdx];
// set the primary variables for the beginning of the current time step.
priVars1 = priVars0[temperatureIdx];
}
// set the primary variables for the beginning of the current time step.
priVars1 = priVars0[temperatureIdx];
}
};
+9 -5
View File
@@ -83,8 +83,8 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
static constexpr bool enableFullyImplicitThermal =
getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
@@ -94,7 +94,7 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
using BrineModule = BlackOilBrineModule<TypeTag, enableBrine>;
using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
using EnergyModule = BlackOilEnergyModule<TypeTag>;
using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;
@@ -191,7 +191,9 @@ public:
}
// deal with energy (if present)
EnergyModule::addStorage(storage, intQuants);
if constexpr (enableFullyImplicitThermal) {
EnergyModule::addStorage(storage, intQuants);
}
// deal with foam (if present)
if constexpr (enableFoam) {
@@ -255,7 +257,9 @@ public:
}
// deal with energy (if present)
EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
if constexpr (enableFullyImplicitThermal) {
EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
}
// deal with foam (if present)
if constexpr (enableFoam) {
@@ -101,9 +101,10 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
static constexpr bool enableFullyImplicitThermal
= (getPropValue<TypeTag, Properties::EnergyModuleType>()
== EnergyModules::FullyImplicitThermal);
static constexpr EnergyModules energyModuleType =
getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr bool enableFullyImplicitThermal =
energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
@@ -121,7 +122,7 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
using EnergyModule = BlackOilEnergyModule<TypeTag>;
using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
@@ -233,7 +234,9 @@ public:
}
// deal with energy (if present)
EnergyModule::addStorage(storage, intQuants);
if constexpr (enableFullyImplicitThermal) {
EnergyModule::addStorage(storage, intQuants);
}
// deal with foam (if present)
if constexpr (enableFoam) {
+31 -17
View File
@@ -358,12 +358,14 @@ private:
static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool waterEnabled = Indices::waterEnabled;
using BioeffectsModule = BlackOilBioeffectsModule<TypeTag, enableBioeffects>;
using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
using EnergyModule = BlackOilEnergyModule<TypeTag>;
using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;
using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
@@ -393,7 +395,9 @@ public:
if constexpr (enablePolymer) {
PolymerModule::registerParameters();
}
EnergyModule::registerParameters();
if constexpr (enableFullyImplicitThermal) {
EnergyModule::registerParameters();
}
if constexpr (enableDiffusion) {
DiffusionModule::registerParameters();
}
@@ -448,12 +452,13 @@ public:
}
}
if (EnergyModule::primaryVarApplies(pvIdx)) {
return EnergyModule::primaryVarName(pvIdx);
}
else {
throw std::logic_error("Invalid primary variable index");
if constexpr (enableFullyImplicitThermal) {
if (EnergyModule::primaryVarApplies(pvIdx)) {
return EnergyModule::primaryVarName(pvIdx);
}
}
throw std::logic_error("Invalid primary variable index");
}
/*!
@@ -485,12 +490,13 @@ public:
}
}
if (EnergyModule::eqApplies(eqIdx)) {
return EnergyModule::eqName(eqIdx);
}
else {
throw std::logic_error("Invalid equation index");
if constexpr (enableFullyImplicitThermal) {
if (EnergyModule::eqApplies(eqIdx)) {
return EnergyModule::eqName(eqIdx);
}
}
throw std::logic_error("Invalid equation index");
}
/*!
@@ -537,8 +543,10 @@ public:
}
// deal with primary variables stemming from the energy module
if (EnergyModule::primaryVarApplies(pvIdx)) {
return EnergyModule::primaryVarWeight(pvIdx);
if constexpr (enableFullyImplicitThermal) {
if (EnergyModule::primaryVarApplies(pvIdx)) {
return EnergyModule::primaryVarWeight(pvIdx);
}
}
// if the primary variable is either the gas saturation, Rs or Rv
@@ -612,7 +620,9 @@ public:
if constexpr (enablePolymer) {
PolymerModule::serializeEntity(asImp_(), outstream, dof);
}
EnergyModule::serializeEntity(asImp_(), outstream, dof);
if constexpr (enableFullyImplicitThermal) {
EnergyModule::serializeEntity(asImp_(), outstream, dof);
}
}
/*!
@@ -664,7 +674,9 @@ public:
if constexpr (enablePolymer) {
PolymerModule::deserializeEntity(asImp_(), instream, dof);
}
EnergyModule::deserializeEntity(asImp_(), instream, dof);
if constexpr (enableFullyImplicitThermal) {
EnergyModule::deserializeEntity(asImp_(), instream, dof);
}
using PVM_G = typename PrimaryVariables::GasMeaning;
using PVM_W = typename PrimaryVariables::WaterMeaning;
@@ -731,7 +743,9 @@ protected:
if constexpr (enablePolymer) {
PolymerModule::registerOutputModules(asImp_(), this->simulator_);
}
EnergyModule::registerOutputModules(asImp_(), this->simulator_);
if constexpr (enableFullyImplicitThermal) {
EnergyModule::registerOutputModules(asImp_(), this->simulator_);
}
if constexpr (enableBioeffects) {
BioeffectsModule::registerOutputModules(asImp_(), this->simulator_);
}
@@ -102,16 +102,18 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag, VectorTy
// component indices from the fluid system
enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr bool enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();
static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
enum { enableMICP = Indices::enableMICP };
enum { gasCompIdx = FluidSystem::gasCompIdx };
enum { waterCompIdx = FluidSystem::waterCompIdx };
@@ -339,7 +341,9 @@ public:
const bool precipitatedSaltPresent = enableSaltPrecipitation ? saltSaturation > 0.0 : false;
const bool oneActivePhases = fluidState.fluidSystem().numActivePhases() == 1;
// deal with the primary variables for the energy extension
EnergyModule::assignPrimaryVars(*this, fluidState);
if constexpr (enableFullyImplicitThermal) {
EnergyModule::assignPrimaryVars(*this, fluidState);
}
// Determine the meaning of the pressure primary variables
// Depending on the phases present, this variable is either interpreted as the
+15 -13
View File
@@ -196,7 +196,8 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
using Indices = GetPropType<TypeTag, Properties::Indices>;
using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using EnergyModule = BlackOilEnergyModule<TypeTag>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
using IndexTraits = typename FluidSystem::IndexTraitsType;
using WellStateType = WellState<Scalar, IndexTraits>;
@@ -554,7 +555,7 @@ protected:
flux *= scalingFactor_;
}
template < class ResidualNBInfo>
template <class ResidualNBInfo>
void computeHeatFluxTerm(const IntensiveQuantitiesTemp& intQuantsIn,
const IntensiveQuantitiesTemp& intQuantsEx,
const ResidualNBInfo& res_nbinfo,
@@ -562,17 +563,18 @@ protected:
{
short interiorDofIdx = 0; // NB
short exteriorDofIdx = 1; // NB
EnergyModule::ExtensiveQuantities::updateEnergy(heatFlux,
interiorDofIdx, // focusDofIndex,
interiorDofIdx,
exteriorDofIdx,
intQuantsIn,
intQuantsEx,
intQuantsIn.fluidStateTemp(),
intQuantsEx.fluidStateTemp(),
res_nbinfo.inAlpha,
res_nbinfo.outAlpha,
res_nbinfo.faceArea);
BlackOilEnergyExtensiveQuantities<TypeTag, energyModuleType>::
updateEnergy(heatFlux,
interiorDofIdx, // focusDofIndex,
interiorDofIdx,
exteriorDofIdx,
intQuantsIn,
intQuantsEx,
intQuantsIn.fluidStateTemp(),
intQuantsEx.fluidStateTemp(),
res_nbinfo.inAlpha,
res_nbinfo.outAlpha,
res_nbinfo.faceArea);
heatFlux *= scalingFactor_*res_nbinfo.faceArea;
}