26 #ifndef EWOMS_PVS_PRIMARY_VARIABLES_HH
27 #define EWOMS_PVS_PRIMARY_VARIABLES_HH
32 #include <opm/material/constraintsolvers/NcpFlash.hpp>
33 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
38 #include <dune/common/fvector.hh>
53 template <
class TypeTag>
58 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) Implementation;
60 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
61 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
62 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
63 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
64 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
65 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
68 enum { pressure0Idx = Indices::pressure0Idx };
69 enum { switch0Idx = Indices::switch0Idx };
75 typedef typename Opm::MathToolbox<Evaluation> Toolbox;
76 typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
78 typedef Opm::NcpFlash<Scalar, FluidSystem> NcpFlash;
82 { Valgrind::SetDefined(*
this); }
89 Valgrind::CheckDefined(value);
90 Valgrind::SetDefined(*
this);
101 Valgrind::SetDefined(*
this);
103 phasePresence_ = value.phasePresence_;
109 template <
class Flu
idState>
111 const MaterialLawParams &matParams,
112 bool isInEquilibrium =
false)
116 for (
int phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
117 assert(std::abs(fluidState.temperature(0) - fluidState.temperature(phaseIdx)) < 1e-30);
123 if (isInEquilibrium) {
130 typename FluidSystem::ParameterCache paramCache;
131 Opm::CompositionalFluidState<Scalar, FluidSystem> fsFlash;
135 fsFlash.assign(fluidState);
138 paramCache.updateAll(fsFlash);
139 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
140 fsFlash.setDensity(phaseIdx, FluidSystem::density(fsFlash, paramCache,
145 ComponentVector globalMolarities(0.0);
146 for (
int compIdx = 0; compIdx < numComponents; ++compIdx) {
147 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
148 globalMolarities[compIdx] +=
149 fsFlash.saturation(phaseIdx) * fsFlash.molarity(phaseIdx, compIdx);
155 NcpFlash::template solve<MaterialLaw>(fsFlash, paramCache, matParams,
167 {
return phasePresence_; }
176 { phasePresence_ = value; }
209 {
return phasePresence & (1 << phaseIdx); }
218 {
return phasePresence_ & (1 << phaseIdx); }
224 {
return ffs(phasePresence_) - 1; }
231 ParentType::operator=(value);
232 phasePresence_ = value.phasePresence_;
242 ParentType::operator=(value);
259 return Toolbox::createConstant(0.0);
261 int varIdx = switch0Idx + phaseIdx - 1;
263 Toolbox::createConstant((*
this)[varIdx]);
264 return Toolbox::createVariable((*
this)[varIdx], varIdx);
270 template <
class Flu
idState>
273 typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
277 EnergyModule::setPriVarTemperatures(*
this, fluidState);
280 (*this)[pressure0Idx] = FsToolbox::value(fluidState.pressure(0));
281 Valgrind::CheckDefined((*
this)[pressure0Idx]);
285 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
289 for (
int compIdx = 0; compIdx < numComponents; ++compIdx) {
290 a -= FsToolbox::value(fluidState.moleFraction(phaseIdx, compIdx));
292 Scalar b = FsToolbox::value(fluidState.saturation(phaseIdx));
295 phasePresence_ |= (1 << phaseIdx);
299 if (phasePresence_ == 0)
300 OPM_THROW(Opm::NumericalProblem,
301 "Phase state was 0, i.e., no fluid is present");
306 for (
int switchIdx = 0; switchIdx < numPhases - 1; ++switchIdx) {
307 int phaseIdx = switchIdx;
308 int compIdx = switchIdx + 1;
309 if (switchIdx >= lowestPhaseIdx)
313 (*this)[switch0Idx + switchIdx] = FsToolbox::value(fluidState.saturation(phaseIdx));
314 Valgrind::CheckDefined((*
this)[switch0Idx + switchIdx]);
317 (*this)[switch0Idx + switchIdx] =
318 FsToolbox::value(fluidState.moleFraction(lowestPhaseIdx, compIdx));
319 Valgrind::CheckDefined((*
this)[switch0Idx + switchIdx]);
325 for (
int compIdx = numPhases - 1; compIdx < numComponents - 1;
327 (*this)[switch0Idx + compIdx] =
328 FsToolbox::value(fluidState.moleFraction(lowestPhaseIdx, compIdx + 1));
329 Valgrind::CheckDefined((*
this)[switch0Idx + compIdx]);
336 void print(std::ostream &os = std::cout)
const
338 os <<
"(p_" << FluidSystem::phaseName(0) <<
" = "
339 << this->operator[](pressure0Idx);
341 for (
int switchIdx = 0; switchIdx < numPhases - 1; ++switchIdx) {
342 int phaseIdx = switchIdx;
343 int compIdx = switchIdx + 1;
344 if (phaseIdx >= lowestPhaseIdx)
349 os <<
", S_" << FluidSystem::phaseName(phaseIdx) <<
" = "
350 << (*this)[switch0Idx + switchIdx];
353 os <<
", x_" << FluidSystem::phaseName(lowestPhaseIdx) <<
"^"
354 << FluidSystem::componentName(compIdx) <<
" = "
355 << (*this)[switch0Idx + switchIdx];
358 for (
int compIdx = numPhases - 1; compIdx < numComponents - 1;
360 os <<
", x_" << FluidSystem::phaseName(lowestPhaseIdx) <<
"^"
361 << FluidSystem::componentName(compIdx + 1) <<
" = "
362 << (*this)[switch0Idx + compIdx];
365 os <<
", phase presence: " <<
static_cast<int>(phasePresence_) << std::flush;
369 unsigned char phasePresence_;
Evaluation explicitSaturationValue(int phaseIdx, int timeIdx) const
Returns an explcitly stored saturation for a given phase.
Definition: pvsprimaryvariables.hh:255
bool phaseIsPresent(int phaseIdx) const
Returns true iff a phase is present for the current phase presence.
Definition: pvsprimaryvariables.hh:217
static bool phaseIsPresent(int phaseIdx, short phasePresence)
Returns true iff a phase is present for a given phase presence.
Definition: pvsprimaryvariables.hh:208
Represents the primary variables used by the a model.
void assignMassConservative(const FluidState &fluidState, const MaterialLawParams &matParams, bool isInEquilibrium=false)
Set the primary variables from an arbitrary fluid state in a mass conservative way.
Definition: pvsprimaryvariables.hh:110
Represents the primary variables used in the primary variable switching compositional model...
Definition: pvsprimaryvariables.hh:54
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:468
ThisType & operator=(Scalar value)
Assignment operator from a scalar value.
Definition: pvsprimaryvariables.hh:240
short phasePresence() const
Return the fluid phases which are present in a given control volume.
Definition: pvsprimaryvariables.hh:166
The indices for the compositional multi-phase primary variable switching model.
ThisType & operator=(const Implementation &value)
Assignment operator from an other primary variables object.
Definition: pvsprimaryvariables.hh:229
void assignNaive(const FluidState &fluidState)
Directly retrieve the primary variables from an arbitrary fluid state.
Definition: pvsprimaryvariables.hh:271
PvsPrimaryVariables(const PvsPrimaryVariables &value)
Default constructor.
Definition: pvsprimaryvariables.hh:99
Represents the primary variables used by the a model.
Definition: fvbaseprimaryvariables.hh:41
Provides the auxiliary methods required for consideration of the energy equation. ...
Definition: energymodule.hh:54
void setPhasePresent(int phaseIdx, bool yesno=true)
Set whether a given indivividual phase should be present or not.
Definition: pvsprimaryvariables.hh:185
Definition: baseauxiliarymodule.hh:35
int lowestPresentPhaseIdx() const
Returns the phase with the lowest index that is present.
Definition: pvsprimaryvariables.hh:223
void setPhasePresence(short value)
Set which fluid phases are present in a given control volume.
Definition: pvsprimaryvariables.hh:175
unsigned char implicitSaturationIdx() const
Returns the index of the phase with's its saturation is determined by the closure condition of satura...
Definition: pvsprimaryvariables.hh:197
PvsPrimaryVariables()
Definition: pvsprimaryvariables.hh:81
void print(std::ostream &os=std::cout) const
Prints the names of the primary variables and their values.
Definition: pvsprimaryvariables.hh:336
Declares the properties required for the compositional multi-phase primary variable switching model...
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
PvsPrimaryVariables(Scalar value)
Constructor with assignment from scalar.
Definition: pvsprimaryvariables.hh:87