OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-25 07:07:42 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #6792 from totto82/speedupTEMP

Browse Source 
Speedup temp
This commit is contained in:
Atgeirr Flø Rasmussen authored and GitHub committed 2026-02-18 16:39:02 +01:00
2 parents 858eb9298d a2b4d9a71a
commit b6eaa425a3
41 files changed
+410 -338
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

flow/flow_biofilm.cpp
-2
View File
@@ -91,13 +91,11 @@ private:
// get the energy module type to determine the equation number
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_blackoil_temp.cpp
+1 -2
View File
@@ -49,8 +49,7 @@ template<class TypeTag>
struct EnableDiffusion<TypeTag, TTag::flowBlackoilTempProblem> { static constexpr bool value = false; };
template<class TypeTag>
struct NumDerivatives<TypeTag, TTag::flowBlackoilTempProblem>
{ static constexpr int value = GetPropType<TypeTag, Properties::Indices>::numDerivatives; };
struct AvoidElementContext<TypeTag, TTag::flowBlackoilTempProblem> { static constexpr bool value = true; };
}}
flow/flow_gasoil.cpp
-2
View File
@@ -62,13 +62,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gasoil_energy.cpp
-2
View File
@@ -49,14 +49,12 @@ private:
// get the energy module type to determine the equation number
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gasoildiffuse.cpp
-2
View File
@@ -64,13 +64,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater.cpp
-2
View File
@@ -59,13 +59,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater_brine.cpp
-2
View File
@@ -77,13 +77,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater_dissolution.cpp
-2
View File
@@ -75,13 +75,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater_dissolution_diffuse.cpp
-2
View File
@@ -78,13 +78,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater_energy.cpp
-2
View File
@@ -78,13 +78,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater_saltprec_energy.cpp
-2
View File
@@ -80,13 +80,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater_saltprec_vapwat.cpp
-2
View File
@@ -80,13 +80,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_gaswater_solvent.cpp
-2
View File
@@ -57,13 +57,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_micp.cpp
-2
View File
@@ -60,13 +60,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_oilwater.cpp
-2
View File
@@ -58,14 +58,12 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_oilwater_brine.cpp
-2
View File
@@ -62,13 +62,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_oilwater_polymer.cpp
-2
View File
@@ -48,13 +48,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_oilwater_polymer_injectivity.cpp
-2
View File
@@ -55,13 +55,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilTwoPhaseIndices<0,
0,
2,
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_onephase.cpp
-2
View File
@@ -55,13 +55,11 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flow/flow_onephase_energy.cpp
-2
View File
@@ -50,14 +50,12 @@ private:
using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
public:
using type = Opm::BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
numEnergyVars,
enableSeqImpEnergy,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
flowexperimental/BlackOilEnergyIntensiveQuantitiesGlobalIndex.hpp
+8 -29
View File
@@ -117,17 +117,17 @@ class BlackOilEnergyIntensiveQuantitiesGlobalIndex<TypeTag, EnergyModules::Seque
using ParamCache = typename FluidSystem::template ParameterCache<Evaluation>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
static constexpr unsigned numPhases = FluidSystem::numPhases;
public:
void updateTemperature_([[maybe_unused]] const Problem& problem,
[[maybe_unused]] const PrimaryVariables& priVars,
[[maybe_unused]] unsigned globalSpaceIndex,
[[maybe_unused]] unsigned timeIdx)
{
auto& fs = this->asImp_().fluidState_;
Evaluation T = Evaluation::createVariable(problem.temperature(globalSpaceIndex, timeIdx), Indices::temperatureIdx);
fs.setTemperature(T);
throw std::logic_error("updateTemperature not implemented "
"SequentialImplicitThermal can not be used with"
"global intensive quantites yet.");
}
void updateEnergyQuantities_([[maybe_unused]] const Problem& problem,
@@ -136,32 +136,11 @@ public:
[[maybe_unused]] unsigned timeIdx,
const ParamCache& paramCache)
{
auto& fs = Parent::asImp_().fluidState_;
// compute the specific enthalpy of the fluids, the specific enthalpy of the rock
// and the thermal conductivity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
fs.setEnthalpy(phaseIdx, h);
}
const auto& solidEnergyLawParams = problem().solidEnergyLawParams(globalSpaceIndex, timeIdx);
this->rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIndex, timeIdx);
this->totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
// Retrieve the rock fraction from the problem
// Usually 1 - porosity, but if pvmult is used to modify porosity
// we will apply the same multiplier to the rock fraction
// i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
// multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
this->rockFraction_ = problem.rockFraction(globalSpaceIndex, timeIdx);
throw std::logic_error("updateEnergyQuantities_ not implemented "
"SequentialImplicitThermal can not be used with"
"global intensive quantites yet.");
}
};
template <class TypeTag>
flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp
+2 -2
View File
@@ -133,8 +133,8 @@ class BlackOilIntensiveQuantitiesGlobalIndex
public:
using FluidState = BlackOilFluidState<Evaluation,
FluidSystem,
energyModuleType == EnergyModules::ConstantTemperature,
(energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
energyModuleType != EnergyModules::NoTemperature,
energyModuleType == EnergyModules::FullyImplicitThermal,
compositionSwitchEnabled,
enableVapwat,
enableBrine,
opm/models/blackoil/blackoilenergymodules.hh
+21 -51
View File
@@ -32,6 +32,7 @@
#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/common/quantitycallbacks.hh>
@@ -523,15 +524,14 @@ class BlackOilEnergyIntensiveQuantities<TypeTag, EnergyModules::SequentialImplic
using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
public:
void updateTemperature_(const Problem& problem, unsigned globalDofIdx, unsigned timeIdx)
{
// update the temperature for output (without derivatives)
auto& fs = asImp_().fluidState_;
const Evaluation T = Evaluation::createVariable(problem.temperature(globalDofIdx, timeIdx), Indices::temperatureIdx);
fs.setTemperature(T);
fs.setTemperature(problem.temperature(globalDofIdx, timeIdx));
}
void updateTemperature_(const ElementContext& elemCtx,
@@ -559,56 +559,30 @@ public:
unsigned dofIdx,
unsigned timeIdx)
{
updateEnergyQuantities_(elemCtx.problem(), elemCtx.globalSpaceIndex(dofIdx, timeIdx), timeIdx);
}
void updateEnergyQuantities_(const Problem& problem,
const unsigned globalSpaceIdx,
const unsigned timeIdx)
{
auto& fs = asImp_().fluidState_;
// compute the specific enthalpy of the fluids, the specific enthalpy of the rock
// and the thermal conductivity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const auto& h = FluidSystem::enthalpy(fs, phaseIdx, problem.pvtRegionIndex(globalSpaceIdx));
fs.setEnthalpy(phaseIdx, h);
}
const auto& solidEnergyLawParams = problem.solidEnergyLawParams(globalSpaceIdx, timeIdx);
rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIdx, timeIdx);
totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
// Retrieve the rock fraction from the problem
// Usually 1 - porosity, but if pvmult is used to modify porosity
// we will apply the same multiplier to the rock fraction
// i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
// multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
rockFraction_ = problem.rockFraction(globalSpaceIdx, timeIdx);
}
const Evaluation& rockInternalEnergy() const
{ return rockInternalEnergy_; }
{
throw std::logic_error("Requested the rock internal energy, which is "
"unavailable because energy is not conserved");
}
const Evaluation& totalThermalConductivity() const
{ return totalThermalConductivity_; }
const Scalar& rockFraction() const
{ return rockFraction_; }
{
throw std::logic_error("Requested the total thermal conductivity, which is "
"unavailable because energy is not conserved");
}
protected:
Implementation& asImp_()
{ return *static_cast<Implementation*>(this); }
Evaluation rockInternalEnergy_;
Evaluation totalThermalConductivity_;
Scalar rockFraction_;
};
template <class TypeTag>
@@ -678,10 +652,9 @@ class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::FullyImplicitThe
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
public:
template<class FluidState>
template<class Evaluation, class FluidState, class IntensiveQuantities>
static void updateEnergy(Evaluation& energyFlux,
const unsigned& focusDofIndex,
const unsigned& inIdx,
@@ -786,7 +759,7 @@ public:
updateEnergyBoundary(energyFlux_, inIq, focusDofIdx, inIdx, alpha, boundaryFs);
}
template <class BoundaryFluidState>
template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
static void updateEnergyBoundary(Evaluation& energyFlux,
const IntensiveQuantities& inIq,
unsigned focusDofIndex,
@@ -839,12 +812,11 @@ template <class TypeTag>
class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::ConstantTemperature>
{
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
public:
template<class FluidState>
template<class Evaluation, class FluidState, class IntensiveQuantities>
static void updateEnergy(Evaluation& /*energyFlux*/,
const unsigned& /*focusDofIndex*/,
const unsigned& /*inIdx*/,
@@ -870,7 +842,7 @@ public:
const BoundaryFluidState&)
{}
template <class BoundaryFluidState>
template <class BoundaryFluidState,class IntensiveQuantities>
static void updateEnergyBoundary(Evaluation& /*heatFlux*/,
const IntensiveQuantities& /*inIq*/,
unsigned /*focusDofIndex*/,
@@ -888,12 +860,11 @@ template <class TypeTag>
class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::SequentialImplicitThermal>
{
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
public:
template<class FluidState>
template<class Evaluation, class FluidState, class IntensiveQuantities>
static void updateEnergy(Evaluation& energyFlux,
const unsigned& focusDofIndex,
const unsigned& inIdx,
@@ -932,7 +903,7 @@ public:
const BoundaryFluidState&)
{ }
template <class BoundaryFluidState>
template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
static void updateEnergyBoundary(Evaluation& /*heatFlux*/,
const IntensiveQuantities& /*inIq*/,
unsigned /*focusDofIndex*/,
@@ -951,11 +922,10 @@ class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::NoTemperature>
{
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
template<class FluidState>
template<class Evaluation, class FluidState, class IntensiveQuantities>
static void updateEnergy(Evaluation& /*energyFlux*/,
const unsigned& /*focusDofIndex*/,
const unsigned& /*inIdx*/,
@@ -981,7 +951,7 @@ public:
const BoundaryFluidState&)
{}
template <class BoundaryFluidState>
template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
static void updateEnergyBoundary(Evaluation& /*heatFlux*/,
const IntensiveQuantities& /*inIq*/,
unsigned /*focusDofIndex*/,
opm/models/blackoil/blackoilintensivequantities.hh
+5 -7
View File
@@ -137,8 +137,8 @@ class BlackOilIntensiveQuantities
public:
using FluidState = BlackOilFluidState<Evaluation,
FluidSystem,
energyModuleType == EnergyModules::ConstantTemperature,
(energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
energyModuleType != EnergyModules::NoTemperature,
energyModuleType == EnergyModules::FullyImplicitThermal,
compositionSwitchEnabled,
enableVapwat,
enableBrine,
@@ -147,9 +147,8 @@ public:
Indices::numPhases>;
using ScalarFluidState = BlackOilFluidState<Scalar,
FluidSystem,
energyModuleType == EnergyModules::ConstantTemperature,
(energyModuleType == EnergyModules::FullyImplicitThermal ||
energyModuleType == EnergyModules::SequentialImplicitThermal),
energyModuleType != EnergyModules::NoTemperature,
energyModuleType == EnergyModules::FullyImplicitThermal,
compositionSwitchEnabled,
enableVapwat,
enableBrine,
@@ -633,8 +632,7 @@ public:
if constexpr (enablePolymer) {
asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
}
if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal ||
energyModuleType == EnergyModules::SequentialImplicitThermal) {
if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx);
}
if constexpr (enableFoam) {
opm/models/blackoil/blackoilmodel.hh
-1
View File
@@ -145,7 +145,6 @@ struct Indices<TypeTag, TTag::BlackOilModel>
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal,
getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
opm/models/blackoil/blackoilonephaseindices.hh
+3 -4
View File
@@ -43,7 +43,6 @@ template<unsigned numSolventsV,
unsigned numExtbosV,
unsigned numPolymersV,
unsigned numEnergyV,
bool enableSequentialImplicitThermal,
bool enableFoam,
bool enableBrine,
unsigned PVOffset,
@@ -105,8 +104,8 @@ struct BlackOilOnePhaseIndices
static constexpr int numEq = numPhases + numSolvents + numExtbos + numPolymers +
numEnergy + numFoam + numBrine + numMICPs;
//! For seqential implicit approches we evaluate the intensive quantities with a larger number of derivatives than equations
static constexpr int numDerivatives = numEq + enableSequentialImplicitThermal;
//! We may want to evaluate with more derivatives that equations, but not this time
static constexpr int numDerivatives = numEq;
//////////////////////////////
// Primary variable indices
@@ -184,7 +183,7 @@ struct BlackOilOnePhaseIndices
//! Index of the primary variable for temperature
static constexpr int temperatureIdx =
(enableSequentialImplicitThermal || enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000;
(enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000;
//////////////////////
// Equation indices
opm/models/blackoil/blackoiltwophaseindices.hh
+3 -4
View File
@@ -43,7 +43,6 @@ template<unsigned numSolventsV,
unsigned numExtbosV,
unsigned numPolymersV,
unsigned numEnergyV,
bool enableSequentialImplicitThermal,
bool enableFoam,
bool enableBrine,
int PVOffset,
@@ -105,8 +104,8 @@ struct BlackOilTwoPhaseIndices
static constexpr int numEq = numPhases + numSolvents + numExtbos + numPolymers +
numEnergy + numFoam + numBrine + numBioComp;
//! For seqential implicit approches we evaluate the intensive quantities with a larger number of derivatives than equations
static constexpr int numDerivatives = numEq + enableSequentialImplicitThermal;
//! We may want to evaluate with more derivatives that equations, but not this time
static constexpr int numDerivatives = numEq;
//////////////////////////////
// Primary variable indices
@@ -177,7 +176,7 @@ struct BlackOilTwoPhaseIndices
//! Index of the primary variable for temperature
static constexpr int temperatureIdx =
(enableSequentialImplicitThermal || enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numBioComp + numFoam + numBrine : - 1000;
(enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numBioComp + numFoam + numBrine : - 1000;
//////////////////////
// Equation indices
opm/models/blackoil/blackoilvariableandequationindices.hh
+3 -4
View File
@@ -39,7 +39,6 @@ template<unsigned numSolventsV,
unsigned numExtbosV,
unsigned numPolymersV,
unsigned numEnergyV,
bool enableSequentialImplicitThermal,
bool enableFoam,
bool enableBrine,
int PVOffset,
@@ -93,8 +92,8 @@ struct BlackOilVariableAndEquationIndices
static constexpr int numEq = numPhases + numSolvents + numExtbos + numPolymers +
numEnergy + numFoam + numBrine;
//! For seqential implicit approches we evaluate the intensive quantities with a larger number of derivatives than equations
static constexpr int numDerivatives = numEq + enableSequentialImplicitThermal;
//! We may want to evaluate with more derivatives that equations, but not this time
static constexpr int numDerivatives = numEq;
////////
@@ -160,7 +159,7 @@ struct BlackOilVariableAndEquationIndices
//! Index of the primary variable for temperature
static constexpr int temperatureIdx =
(enableSequentialImplicitThermal || enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
enableFullyImplicitThermal ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
////////
opm/models/discretization/common/tpfalinearizer.hh
+5 -2
View File
@@ -680,7 +680,8 @@ private:
// If DISPERC is in the deck, we initialize the sparse table here as well.
const bool anyFlows = simulator_().problem().eclWriter().outputModule().getFlows().anyFlows();
const auto& blockFlows = simulator_().problem().eclWriter().outputModule().getFlows().blockFlows();
const bool anyFlores = simulator_().problem().eclWriter().outputModule().getFlows().anyFlores();
const bool isTemp = simulator_().vanguard().eclState().getSimulationConfig().isTemp();
const bool anyFlores = simulator_().problem().eclWriter().outputModule().getFlows().anyFlores() || isTemp;
const bool dispersionActive = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
if (((!(anyFlows || !blockFlows.empty()) || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && (!dispersionActive && !enableBioeffects)) {
return;
@@ -800,7 +801,9 @@ public:
OPM_TIMEBLOCK(updateFlows);
const bool enableFlows = simulator_().problem().eclWriter().outputModule().getFlows().hasFlows();
const auto& blockFlows = simulator_().problem().eclWriter().outputModule().getFlows().blockFlows();
const bool enableFlores = simulator_().problem().eclWriter().outputModule().getFlows().hasFlores();
// We reuse the fluxes in the TEMP option
const bool isTemp = simulator_().vanguard().eclState().getSimulationConfig().isTemp();
const bool enableFlores = simulator_().problem().eclWriter().outputModule().getFlows().hasFlores() || isTemp;
if (!enableFlows && !enableFlores && blockFlows.empty()) {
return;
}
opm/simulators/aquifers/AquiferAnalytical.hpp
+2 -2
View File
@@ -78,8 +78,8 @@ public:
using FluidState = BlackOilFluidState<Eval,
FluidSystem,
energyModuleType == EnergyModules::ConstantTemperature,
(energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
energyModuleType != EnergyModules::NoTemperature,
energyModuleType == EnergyModules::FullyImplicitThermal,
BlackoilIndices::gasEnabled,
enableVapwat,
enableBrine,
opm/simulators/flow/EquilInitializer.hpp
+3 -3
View File
@@ -88,8 +88,8 @@ public:
// internal energy!
using ScalarFluidState = BlackOilFluidState<Scalar,
FluidSystem,
energyModuleType == EnergyModules::ConstantTemperature,
(energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
energyModuleType != EnergyModules::NoTemperature,
energyModuleType == EnergyModules::FullyImplicitThermal,
enableDissolution,
enableVapwat,
enableBrine,
@@ -173,7 +173,7 @@ public:
const auto& rho = FluidSystem::density(fluidState, phaseIdx, regionIdx);
fluidState.setDensity(phaseIdx, rho);
if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal) {
if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, regionIdx);
fluidState.setEnthalpy(phaseIdx, h);
}
opm/simulators/flow/FlowProblemBlackoil.hpp
+1 -1
View File
@@ -857,7 +857,7 @@ public:
const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
fluidState.setDensity(phaseIdx, rho);
if constexpr (energyModuleType == EnergyModules::SequentialImplicitThermal || energyModuleType == EnergyModules::FullyImplicitThermal) {
if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
fluidState.setEnthalpy(phaseIdx, h);
}
opm/simulators/flow/GenericTemperatureModel.hpp
+3 -2
View File
@@ -51,7 +51,9 @@ template<class Grid, class GridView, class DofMapper, class Stencil, class Fluid
class GenericTemperatureModel
{
public:
using EnergyMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
// the jacobian matrix
using MatrixBlockTemp = MatrixBlock<Scalar, 1, 1>;
using EnergyMatrix = Dune::BCRSMatrix<MatrixBlockTemp>;
using EnergyVector = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
static constexpr int dimWorld = Grid::dimensionworld;
@@ -85,7 +87,6 @@ protected:
const DofMapper& dofMapper_;
EnergyVector energyVector_;
std::unique_ptr<EnergyMatrix> energyMatrix_;
std::vector<Scalar> temperature_;
std::vector<Scalar> energy_rates_;
bool doTemp_{false};
opm/simulators/flow/GenericTemperatureModel_impl.hpp
+3 -39
View File
@@ -115,48 +115,12 @@ void GenericTemperatureModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>
doInit(std::size_t numGridDof)
{
doTemp_ = eclState_.getSimulationConfig().isTemp();
temperature_.resize(numGridDof);
energyVector_.resize(numGridDof);
// allocate matrix for storing the Jacobian of the temperature residual
energyMatrix_ = std::make_unique<EnergyMatrix>(numGridDof, numGridDof, EnergyMatrix::random);
// find the sparsity pattern of the temperature matrix
using NeighborSet = std::set<unsigned>;
std::vector<NeighborSet> neighbors(numGridDof);
Stencil stencil(gridView_, dofMapper_);
for (const auto& elem : elements(gridView_)) {
stencil.update(elem);
for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
neighbors[myIdx].insert(neighborIdx);
}
}
}
// allocate space for the rows of the matrix
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
energyMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
}
energyMatrix_->endrowsizes();
// fill the rows with indices. each degree of freedom talks to
// all of its neighbors. (it also talks to itself since
// degrees of freedom are sometimes quite egocentric.)
for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
for (; nIt != nEndIt; ++nIt) {
energyMatrix_->addindex(dofIdx, *nIt);
}
}
energyMatrix_->endindices();
if (!doTemp_)
return;
energyVector_.resize(numGridDof);
maxTempChange_ = Parameters::Get<Parameters::MaxTemperatureChange<Scalar>>();
}
opm/simulators/flow/TemperatureModel.hpp
+294 -87
View File
@@ -37,6 +37,8 @@
#include <opm/simulators/flow/GenericTemperatureModel.hpp>
#include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
#include <opm/simulators/aquifers/AquiferGridUtils.hpp>
#include <opm/models/discretization/common/tpfalinearizer.hh>
#include <opm/simulators/linalg/istlsparsematrixadapter.hh>
#include <algorithm>
#include <cassert>
@@ -58,6 +60,110 @@ namespace Opm {
template<typename Scalar, typename IndexTraits> class WellState;
template <class TypeTag>
class BlackOilEnergyIntensiveQuantitiesTemp
{
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
enum { enableDisgasInWater = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
public:
using EvaluationTemp = DenseAd::Evaluation<Scalar, 1>;
using FluidStateTemp = BlackOilFluidState<EvaluationTemp,
FluidSystem,
true, //store temperature
true, // store enthalpy
compositionSwitchEnabled,
enableVapwat,
enableBrine,
enableSaltPrecipitation,
enableDisgasInWater,
Indices::numPhases>;
void updateTemperature_(const Problem& problem, unsigned globalDofIdx, unsigned timeIdx)
{
const EvaluationTemp T = EvaluationTemp::createVariable(problem.temperature(globalDofIdx, timeIdx), 0);
fluidState_.setTemperature(T);
}
void updateEnergyQuantities_(const Problem& problem,
const unsigned globalSpaceIdx,
const unsigned timeIdx)
{
// compute the specific enthalpy of the fluids, the specific enthalpy of the rock
// and the thermal conductivity coefficients
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const auto& h = FluidSystem::enthalpy(fluidState_, phaseIdx, problem.pvtRegionIndex(globalSpaceIdx));
fluidState_.setEnthalpy(phaseIdx, h);
}
const auto& solidEnergyLawParams = problem.solidEnergyLawParams(globalSpaceIdx, timeIdx);
rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fluidState_);
const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIdx, timeIdx);
totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fluidState_);
// Retrieve the rock fraction from the problem
// Usually 1 - porosity, but if pvmult is used to modify porosity
// we will apply the same multiplier to the rock fraction
// i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
// multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
rockFraction_ = problem.rockFraction(globalSpaceIdx, timeIdx);
}
const EvaluationTemp& rockInternalEnergy() const
{ return rockInternalEnergy_; }
const EvaluationTemp& totalThermalConductivity() const
{ return totalThermalConductivity_; }
const Scalar& rockFraction() const
{ return rockFraction_; }
const FluidStateTemp& fluidStateTemp() const
{ return fluidState_; }
template <class FluidState>
void setFluidState(const FluidState& fs)
{
// copy the needed part of the fluid state
fluidState_.setPvtRegionIndex(fs.pvtRegionIndex());
fluidState_.setRs(getValue(fs.Rs()));
fluidState_.setRv(getValue(fs.Rv()));
for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
fluidState_.setPressure(phaseIdx, getValue(fs.pressure(phaseIdx)));
fluidState_.setDensity(phaseIdx, getValue(fs.density(phaseIdx)));
fluidState_.setSaturation(phaseIdx, getValue(fs.saturation(phaseIdx)));
fluidState_.setInvB(phaseIdx, getValue(fs.invB(phaseIdx)));
}
}
protected:
EvaluationTemp rockInternalEnergy_;
EvaluationTemp totalThermalConductivity_;
FluidStateTemp fluidState_;
Scalar rockFraction_;
};
/*!
* \ingroup BlackOilSimulator
*
@@ -84,24 +190,38 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
using Stencil = GetPropType<TypeTag, Properties::Stencil>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using TemperatureEvaluation = DenseAd::Evaluation<Scalar,1>;
using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using EnergyModule = BlackOilEnergyModule<TypeTag>;
using IndexTraits = typename FluidSystem::IndexTraitsType;
using WellStateType = WellState<Scalar, IndexTraits>;
using IntensiveQuantitiesTemp = BlackOilEnergyIntensiveQuantitiesTemp<TypeTag>;
using FluidStateTemp = typename IntensiveQuantitiesTemp::FluidStateTemp;
using Evaluation = typename IntensiveQuantitiesTemp::EvaluationTemp;
using EnergyMatrix = typename BaseType::EnergyMatrix;
using EnergyVector = typename BaseType::EnergyVector;
using MatrixBlockTemp = typename BaseType::MatrixBlockTemp;
using SpareMatrixEnergyAdapter = typename Linear::IstlSparseMatrixAdapter<MatrixBlockTemp>;
//using ResidualNBInfo = typename LocalResidual::ResidualNBInfo;
struct ResidualNBInfo
{
double faceArea;
double inAlpha;
double outAlpha;
};
using NeighborInfoCPU = NeighborInfoStruct<ResidualNBInfo, MatrixBlockTemp>;
enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
enum { numPhases = FluidSystem::numPhases };
enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
static constexpr int temperatureIdx = 0;
public:
explicit TemperatureModel(Simulator& simulator)
@@ -134,6 +254,48 @@ public:
// find and store the overlap cells
const auto& elemMapper = simulator_.model().elementMapper();
detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
// set the scaling factor
scalingFactor_ = getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
// find the sparsity pattern of the temperature matrix
using NeighborSet = std::set<unsigned>;
std::vector<NeighborSet> neighbors(numCells);
Stencil stencil(this->gridView_, this->dofMapper_);
neighborInfo_.reserve(numCells, 6 * numCells);
std::vector<NeighborInfoCPU> loc_nbinfo;
for (const auto& elem : elements(this->gridView_)) {
stencil.update(elem);
for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
const unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
loc_nbinfo.resize(stencil.numDof() - 1); // Do not include the primary dof in neighborInfo_
for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
const unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
neighbors[myIdx].insert(neighborIdx);
if (dofIdx > 0) {
const auto scvfIdx = dofIdx - 1;
const auto& scvf = stencil.interiorFace(scvfIdx);
const Scalar area = scvf.area();
Scalar inAlpha = simulator_.problem().thermalHalfTransmissibility(myIdx, neighborIdx);
Scalar outAlpha = simulator_.problem().thermalHalfTransmissibility(neighborIdx, myIdx);
ResidualNBInfo nbinfo{area, inAlpha, outAlpha};
loc_nbinfo[dofIdx - 1] = NeighborInfoCPU{neighborIdx, nbinfo, nullptr};
}
}
neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
}
}
// allocate matrix for storing the Jacobian of the temperature residual
energyMatrix_ = std::make_unique<SpareMatrixEnergyAdapter>(simulator_);
diagMatAddress_.resize(numCells);
energyMatrix_->reserve(neighbors);
for (unsigned globI = 0; globI < numCells; globI++) {
const auto& nbInfos = neighborInfo_[globI];
diagMatAddress_[globI] = energyMatrix_->blockAddress(globI, globI);
for (auto& nbInfo : nbInfos) {
nbInfo.matBlockAddress = energyMatrix_->blockAddress(nbInfo.neighbor, globI);
}
}
}
void beginTimeStep()
@@ -142,10 +304,13 @@ public:
return;
}
// We copy the intensive quantities here to make it possible to update them
// We use the specialized intensive quantities here with only the temperature derivative
const unsigned int numCells = simulator_.model().numTotalDof();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned globI = 0; globI < numCells; ++globI) {
intQuants_[globI] = simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0);
intQuants_[globI].setFluidState(simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0).fluidState());
intQuants_[globI].updateTemperature_(simulator_.problem(), globI, /*timeIdx*/ 0);
intQuants_[globI].updateEnergyQuantities_(simulator_.problem(), globI, /*timeIdx*/ 0);
}
@@ -164,10 +329,13 @@ public:
return;
}
// We copy the intensive quantities here to make it possible to update them
// We use the specialized intensive quantities here with only the temperature derivative
const unsigned int numCells = simulator_.model().numTotalDof();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned globI = 0; globI < numCells; ++globI) {
intQuants_[globI] = simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0);
intQuants_[globI].setFluidState(simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0).fluidState());
intQuants_[globI].updateTemperature_(simulator_.problem(), globI, /*timeIdx*/ 0);
intQuants_[globI].updateEnergyQuantities_(simulator_.problem(), globI, /*timeIdx*/ 0);
}
@@ -179,6 +347,13 @@ public:
auto& ws = wellState.well(wellID);
ws.energy_rate = this->energy_rates_[wellID];
}
// Update well temperature
const auto& wellPtrs = simulator_.problem().wellModel().localNonshutWells();
for (const auto& wellPtr : wellPtrs) {
auto& ws = wellState.well(wellPtr->name());
this->computeWellTemperature(*wellPtr, ws);
}
}
/*!
@@ -241,13 +416,16 @@ protected:
{
const unsigned int numCells = simulator_.model().numTotalDof();
EnergyVector dx(numCells);
bool conv = this->linearSolve_(*this->energyMatrix_, dx, this->energyVector_);
bool conv = this->linearSolve_(this->energyMatrix_->istlMatrix(), dx, this->energyVector_);
if (!conv) {
if (simulator_.gridView().comm().rank() == 0) {
OpmLog::warning("Temp model: Linear solver did not converge. Temperature values not updated.");
}
}
else {
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned globI = 0; globI < numCells; ++globI) {
this->temperature_[globI] -= std::clamp(dx[globI][0], -this->maxTempChange_, this->maxTempChange_);
intQuants_[globI].updateTemperature_(simulator_.problem(), globI, /*timeIdx*/ 0);
@@ -313,9 +491,9 @@ protected:
void computeStorageTerm(unsigned globI, LhsEval& storage)
{
const auto& intQuants = intQuants_[globI];
const auto& poro = decay<LhsEval>(intQuants.porosity());
const auto& poro = getValue(simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0).porosity());
// accumulate the internal energy of the fluids
const auto& fs = intQuants.fluidState();
const auto& fs = intQuants.fluidStateTemp();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@@ -332,20 +510,15 @@ protected:
const Scalar rockFraction = intQuants.rockFraction();
const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
storage += rockFraction*uRock;
storage*= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
storage*= scalingFactor_;
}
template < class ResidualNBInfo>
void computeFluxTerm(unsigned globI, unsigned globJ,
const ResidualNBInfo& res_nbinfo,
template <class RateVector>
void computeFluxTerm(const FluidStateTemp& fsIn,
const FluidStateTemp& fsEx,
const RateVector& darcyFlux,
Evaluation& flux)
{
const IntensiveQuantities& intQuantsIn = intQuants_[globI];
const IntensiveQuantities& intQuantsEx = intQuants_[globJ];
RateVector tmp(0.0); //not used
RateVector darcyFlux(0.0);
LocalResidual::computeFlux(tmp, darcyFlux, globI, globJ, intQuantsIn, intQuantsEx,
res_nbinfo, simulator_.problem().moduleParams());
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@@ -353,14 +526,12 @@ protected:
const unsigned activeCompIdx =
FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
bool inIsUp = darcyFlux[activeCompIdx] > 0;
const IntensiveQuantities& up = inIsUp ? intQuantsIn : intQuantsEx;
const auto& fs = up.fluidState();
const auto& fs = inIsUp ? fsIn : fsEx;
if (inIsUp) {
flux += fs.enthalpy(phaseIdx)
* fs.density(phaseIdx)
* darcyFlux[activeCompIdx];
* getValue(darcyFlux[activeCompIdx]);
}
else {
flux += getValue(fs.enthalpy(phaseIdx))
@@ -368,18 +539,15 @@ protected:
* getValue(darcyFlux[activeCompIdx]);
}
}
flux *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
flux *= scalingFactor_;
}
template < class ResidualNBInfo>
void computeHeatFluxTerm(unsigned globI, unsigned globJ,
void computeHeatFluxTerm(const IntensiveQuantitiesTemp& intQuantsIn,
const IntensiveQuantitiesTemp& intQuantsEx,
const ResidualNBInfo& res_nbinfo,
Evaluation& heatFlux)
{
const IntensiveQuantities& intQuantsIn = intQuants_[globI];
const IntensiveQuantities& intQuantsEx = intQuants_[globJ];
const Scalar inAlpha = simulator_.problem().thermalHalfTransmissibility(globI, globJ);
const Scalar outAlpha = simulator_.problem().thermalHalfTransmissibility(globJ, globI);
short interiorDofIdx = 0; // NB
short exteriorDofIdx = 1; // NB
EnergyModule::ExtensiveQuantities::updateEnergy(heatFlux,
@@ -388,55 +556,79 @@ protected:
exteriorDofIdx,
intQuantsIn,
intQuantsEx,
intQuantsIn.fluidState(),
intQuantsEx.fluidState(),
inAlpha,
outAlpha,
intQuantsIn.fluidStateTemp(),
intQuantsEx.fluidStateTemp(),
res_nbinfo.inAlpha,
res_nbinfo.outAlpha,
res_nbinfo.faceArea);
heatFlux *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>()*res_nbinfo.faceArea;
heatFlux *= scalingFactor_*res_nbinfo.faceArea;
}
void assembleEquations()
{
this->energyVector_ = 0.0;
(*this->energyMatrix_) = 0.0;
Scalar dt = simulator_.timeStepSize();
const unsigned int numCells = simulator_.model().numTotalDof();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned globI = 0; globI < numCells; ++globI) {
this->energyVector_[globI] = 0.0;
energyMatrix_->clearRow(globI, 0.0);
}
Scalar dt = simulator_.timeStepSize();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned globI = 0; globI < numCells; ++globI) {
MatrixBlockTemp bMat;
Scalar volume = simulator_.model().dofTotalVolume(globI);
Scalar storefac = volume / dt;
Evaluation storage = 0.0;
computeStorageTerm(globI, storage);
this->energyVector_[globI] += storefac * ( getValue(storage) - storage1_[globI] );
(*this->energyMatrix_)[globI][globI][0][0] += storefac * storage.derivative(Indices::temperatureIdx);
bMat[0][0] = storefac * storage.derivative(temperatureIdx);
*diagMatAddress_[globI] += bMat;
}
const auto& neighborInfo = simulator_.model().linearizer().getNeighborInfo();
#ifdef _OPENMP
#pragma omp parallel for
#endif
const auto& floresInfo = this->simulator_.problem().model().linearizer().getFloresInfo();
const bool enableDriftCompensation = Parameters::Get<Parameters::EnableDriftCompensationTemp>();
const auto& problem = simulator_.problem();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (unsigned globI = 0; globI < numCells; ++globI) {
const auto& nbInfos = neighborInfo[globI];
const auto& nbInfos = neighborInfo_[globI];
const auto& floresInfos = floresInfo[globI];
int loc = 0;
const auto& intQuantsIn = intQuants_[globI];
MatrixBlockTemp bMat;
for (const auto& nbInfo : nbInfos) {
unsigned globJ = nbInfo.neighbor;
const auto& intQuantsEx = intQuants_[globJ];
assert(globJ != globI);
const auto& darcyflux = floresInfos[loc].flow;
// compute convective flux
Evaluation flux = 0.0;
computeFluxTerm(globI, globJ, nbInfo.res_nbinfo, flux);
this->energyVector_[globI] += getValue(flux);
(*this->energyMatrix_)[globI][globI][0][0] += flux.derivative(Indices::temperatureIdx);
(*this->energyMatrix_)[globJ][globI][0][0] -= flux.derivative(Indices::temperatureIdx);
computeFluxTerm(intQuantsIn.fluidStateTemp(), intQuantsEx.fluidStateTemp(), darcyflux, flux);
// compute conductive flux
Evaluation heatFlux = 0.0;
computeHeatFluxTerm(globI, globJ, nbInfo.res_nbinfo, heatFlux);
computeHeatFluxTerm(intQuantsIn, intQuantsEx, nbInfo.res_nbinfo, heatFlux);
heatFlux += flux;
this->energyVector_[globI] += getValue(heatFlux);
(*this->energyMatrix_)[globI][globI][0][0] += heatFlux.derivative(Indices::temperatureIdx);
(*this->energyMatrix_)[globJ][globI][0][0] -= heatFlux.derivative(Indices::temperatureIdx);
bMat[0][0] = heatFlux.derivative(temperatureIdx);
*diagMatAddress_[globI] += bMat;
bMat *= -1.0;
//SparseAdapter syntax: jacobian_->addToBlock(globJ, globI, bMat);
*nbInfo.matBlockAddress += bMat;
loc++;
}
if (enableDriftCompensation) {
auto dofDriftRate = problem.drift()[globI]/dt;
const auto& fs = intQuantsIn.fluidStateTemp();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
const unsigned activeCompIdx =
FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
auto drift_hrate = dofDriftRate[activeCompIdx]*getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx)) / getValue(fs.invB(phaseIdx));
this->energyVector_[globI] -= drift_hrate*scalingFactor_;
}
}
}
@@ -446,33 +638,19 @@ protected:
this->assembleEquationWell(*wellPtr);
}
const auto& problem = simulator_.problem();
bool enableDriftCompensation = Parameters::Get<Parameters::EnableDriftCompensationTemp>();
if (enableDriftCompensation) {
for (unsigned globalDofIdx = 0; globalDofIdx < numCells; ++globalDofIdx) {
auto dofDriftRate = problem.drift()[globalDofIdx]/dt;
const auto& fs = intQuants_[globalDofIdx].fluidState();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
const unsigned activeCompIdx =
FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
auto drift_hrate = dofDriftRate[activeCompIdx]*getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx)) / getValue(fs.invB(phaseIdx));
this->energyVector_[globalDofIdx] -= drift_hrate*getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
}
}
}
if (simulator_.gridView().comm().size() > 1) {
// Set dirichlet conditions for overlapping cells
// loop over precalculated overlap rows and columns
for (const auto row : overlapRows_) {
// Zero out row.
(*this->energyMatrix_)[row] = 0.0;
//diagonal block set to diag(1.0).
(*this->energyMatrix_)[row][row][0][0] = 1.0;
}
}
// For standard ilu0 we need to set the overlapping cells to identity. But since we use "paroverilu0"
// here this is not needed.
// if (simulator_.gridView().comm().size() > 1) {
// // Set dirichlet conditions for overlapping cells
// // loop over precalculated overlap rows and columns
// for (const auto row : overlapRows_) {
// // Zero out row.
// (*this->energyMatrix_)[row] = 0.0;
//
// //diagonal block set to diag(1.0).
// (*this->energyMatrix_)[row][row][0][0] = 1.0;
// }
//}
}
template<class Well>
@@ -482,9 +660,10 @@ protected:
std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value();
const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
this->energy_rates_[well_index] = 0.0;
MatrixBlockTemp bMat;
for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
const auto globI = ws.perf_data.cell_index[i];
auto fs = intQuants_[globI].fluidState(); //copy to make it possible to change the temp in the injector
auto fs = intQuants_[globI].fluidStateTemp(); //copy to make it possible to change the temp in the injector
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@@ -513,18 +692,46 @@ protected:
rate *= fs.enthalpy(phaseIdx) * getValue(fs.density(phaseIdx)) / getValue(fs.invB(phaseIdx));
}
this->energy_rates_[well_index] += getValue(rate);
rate *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
rate *= scalingFactor_;
this->energyVector_[globI] -= getValue(rate);
(*this->energyMatrix_)[globI][globI][0][0] -= rate.derivative(Indices::temperatureIdx);
bMat[0][0] = -rate.derivative(temperatureIdx);
*diagMatAddress_[globI] += bMat;
}
}
}
template<class Well, class SingleWellState>
void computeWellTemperature(const Well& well, SingleWellState& ws)
{
if (well.isInjector()) {
if (ws.status != WellStatus::STOP) {
ws.temperature = well.wellEcl().inj_temperature();
return;
}
}
const int np = simulator_.problem().wellModel().wellState().numPhases();
std::array<Scalar,2> weighted{0.0,0.0};
auto& [weighted_temperature, total_weight] = weighted;
for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
const auto globI = ws.perf_data.cell_index[i];
const auto& fs = intQuants_[globI].fluidStateTemp();
Scalar weight_factor = simulator_.problem().wellModel().computeTemperatureWeightFactor(i, np, fs, ws);
total_weight += weight_factor;
weighted_temperature += weight_factor * fs.temperature(/*phaseIdx*/0).value();
}
//simulator_.gridView().comm().sum(weighted.data(), 2);
ws.temperature = weighted_temperature / total_weight;
}
const Simulator& simulator_;
EnergyVector storage1_;
std::vector<IntensiveQuantities> intQuants_;
std::vector<IntensiveQuantitiesTemp> intQuants_;
SparseTable<NeighborInfoCPU> neighborInfo_{};
std::vector<MatrixBlockTemp*> diagMatAddress_{};
std::unique_ptr<SpareMatrixEnergyAdapter> energyMatrix_;
std::vector<int> overlapRows_;
std::vector<int> interiorRows_;
Scalar scalingFactor_{1.0};
};
// need for the old linearizer
opm/simulators/utils/InstantiationIndicesMacros.hpp
+24 -25
View File
@@ -63,30 +63,29 @@
template class CLASS<FS<T>,__VA_ARGS__>;
#define INSTANTIATE_TYPE_INDICES(CLASS, T) \
INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,false,0u,1u,5u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,true,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,false,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,false,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,true,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,false,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,true,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<1u,0u,0u,0u,false,false,false,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,0u,2u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,true,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,true,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,false,true,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,false,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,1u,0u,0u,false,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,1u,0u,false,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,false,true,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,false,true,false,0u,0u>)
INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,true,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<1u,0u,0u,0u,false,false,0u,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,2u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,true,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,true,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,1u,0u,0u,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,1u,0u,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,false,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,true,0u,0u>) \
INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,true,false,0u,0u>)
#endif
opm/simulators/wells/BlackoilWellModel.hpp
+22
View File
@@ -460,6 +460,28 @@ template<class Scalar> class WellContributions;
const ModelParameters& param() const
{ return param_; }
template<class FluidState, class SingleWellState>
static Scalar computeTemperatureWeightFactor(const int perf_index, const int np, const FluidState& fs, const SingleWellState& ws)
{
const auto& perf_phase_rate = ws.perf_data.phase_rates;
// we only have one temperature pr cell any phaseIdx will do
Scalar cellTemperatures = fs.temperature(/*phaseIdx*/0).value();
Scalar weight_factor = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
Scalar cellInternalEnergy = fs.enthalpy(phaseIdx).value() -
fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
Scalar cellBinv = fs.invB(phaseIdx).value();
Scalar cellDensity = fs.density(phaseIdx).value();
Scalar perfPhaseRate = perf_phase_rate[perf_index*np + phaseIdx];
weight_factor += cellDensity * (perfPhaseRate / cellBinv) * (cellInternalEnergy / cellTemperatures);
}
return (std::abs(weight_factor) + 1e-13);
}
protected:
Simulator& simulator_;
opm/simulators/wells/BlackoilWellModel_impl.hpp
+5 -31
View File
@@ -2173,13 +2173,8 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
computeWellTemperature()
{
if constexpr (energyModuleType_ == EnergyModules::FullyImplicitThermal ||
energyModuleType_ == EnergyModules::SequentialImplicitThermal) {
int np = this->numPhases();
Scalar cellInternalEnergy;
Scalar cellBinv;
Scalar cellDensity;
Scalar perfPhaseRate;
if constexpr (energyModuleType_ == EnergyModules::FullyImplicitThermal) {
const int np = this->numPhases();
const int nw = this->numLocalWells();
for (auto wellID = 0*nw; wellID < nw; ++wellID) {
const Well& well = this->wells_ecl_[wellID];
@@ -2190,38 +2185,17 @@ namespace Opm {
continue;
}
}
std::array<Scalar,2> weighted{0.0,0.0};
auto& [weighted_temperature, total_weight] = weighted;
auto& well_info = this->local_parallel_well_info_[wellID].get();
auto& perf_data = ws.perf_data;
auto& perf_phase_rate = perf_data.phase_rates;
const auto& well_info = this->local_parallel_well_info_[wellID].get();
using int_type = decltype(this->well_perf_data_[wellID].size());
for (int_type perf = 0, end_perf = this->well_perf_data_[wellID].size(); perf < end_perf; ++perf) {
const int cell_idx = this->well_perf_data_[wellID][perf].cell_index;
const auto& intQuants = simulator_.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// we on only have one temperature pr cell any phaseIdx will do
Scalar cellTemperatures = fs.temperature(/*phaseIdx*/0).value();
Scalar weight_factor = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
cellInternalEnergy = fs.enthalpy(phaseIdx).value() -
fs.pressure(phaseIdx).value() / fs.density(phaseIdx).value();
cellBinv = fs.invB(phaseIdx).value();
cellDensity = fs.density(phaseIdx).value();
perfPhaseRate = perf_phase_rate[perf*np + phaseIdx];
weight_factor += cellDensity * perfPhaseRate / cellBinv * cellInternalEnergy / cellTemperatures;
}
weight_factor = std::abs(weight_factor) + 1e-13;
Scalar weight_factor = computeTemperatureWeightFactor(perf, np, fs, ws);
total_weight += weight_factor;
weighted_temperature += weight_factor * cellTemperatures;
weighted_temperature += weight_factor * fs.temperature(/*phaseIdx*/0).value();
}
well_info.communication().sum(weighted.data(), 2);
this->wellState().well(wellID).temperature = weighted_temperature / total_weight;
opm/simulators/wells/WellInterface.hpp
+2 -2
View File
@@ -129,8 +129,8 @@ public:
// For the conversion between the surface volume rate and reservoir voidage rate
using FluidState = BlackOilFluidState<Eval,
FluidSystem,
energyModuleType == EnergyModules::ConstantTemperature,
(energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
energyModuleType != EnergyModules::NoTemperature,
energyModuleType == EnergyModules::FullyImplicitThermal,
Indices::compositionSwitchIdx >= 0,
has_watVapor,
has_brine,
python/simulators/PyGasWaterSimulator.cpp
-1
View File
@@ -46,7 +46,6 @@ namespace Opm {
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal,
getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,
python/simulators/PyOnePhaseSimulator.cpp
-1
View File
@@ -45,7 +45,6 @@ namespace Opm::Properties {
getPropValue<TypeTag, Properties::EnableExtbo>(),
getPropValue<TypeTag, Properties::EnablePolymer>(),
getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal,
getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal,
getPropValue<TypeTag, Properties::EnableFoam>(),
getPropValue<TypeTag, Properties::EnableBrine>(),
/*PVOffset=*/0,