pvsprimaryvariables.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_PVS_PRIMARY_VARIABLES_HH
29 #define EWOMS_PVS_PRIMARY_VARIABLES_HH
30 
31 #include "pvsindices.hh"
32 #include "pvsproperties.hh"
33 
36 
37 #include <opm/material/constraintsolvers/NcpFlash.hpp>
38 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
39 #include <opm/common/Valgrind.hpp>
40 #include <opm/common/ErrorMacros.hpp>
41 #include <opm/common/Exceptions.hpp>
42 
43 #include <dune/common/fvector.hh>
44 
45 #include <iostream>
46 
47 namespace Ewoms {
48 
58 template <class TypeTag>
60 {
61  typedef FvBasePrimaryVariables<TypeTag> ParentType;
63  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) Implementation;
64 
65  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
66  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
67  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
68  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
69  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
70  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
71 
72  // primary variable indices
73  enum { pressure0Idx = Indices::pressure0Idx };
74  enum { switch0Idx = Indices::switch0Idx };
75 
76  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
77  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
78  enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
79 
80  typedef typename Opm::MathToolbox<Evaluation> Toolbox;
81  typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
83  typedef Opm::NcpFlash<Scalar, FluidSystem> NcpFlash;
84 
85 public:
86  PvsPrimaryVariables() : ParentType()
87  { Opm::Valgrind::SetDefined(*this); }
88 
92  explicit PvsPrimaryVariables(Scalar value) : ParentType(value)
93  {
94  Opm::Valgrind::CheckDefined(value);
95  Opm::Valgrind::SetDefined(*this);
96 
97  phasePresence_ = 0;
98  }
99 
104  PvsPrimaryVariables(const PvsPrimaryVariables& value) : ParentType(value)
105  {
106  Opm::Valgrind::SetDefined(*this);
107 
108  phasePresence_ = value.phasePresence_;
109  }
110 
114  template <class FluidState>
115  void assignMassConservative(const FluidState& fluidState,
116  const MaterialLawParams& matParams,
117  bool isInEquilibrium = false)
118  {
119 #ifndef NDEBUG
120  // make sure the temperature is the same in all fluid phases
121  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
122  assert(std::abs(fluidState.temperature(0) - fluidState.temperature(phaseIdx)) < 1e-30);
123  }
124 #endif // NDEBUG
125 
126  // for the equilibrium case, we don't need complicated
127  // computations.
128  if (isInEquilibrium) {
129  assignNaive(fluidState);
130  return;
131  }
132 
133  // use a flash calculation to calculate a fluid state in
134  // thermodynamic equilibrium
135  typename FluidSystem::template ParameterCache<Scalar> paramCache;
136  Opm::CompositionalFluidState<Scalar, FluidSystem> fsFlash;
137 
138  // use the externally given fluid state as initial value for
139  // the flash calculation
140  fsFlash.assign(fluidState);
141 
142  // calculate the phase densities
143  paramCache.updateAll(fsFlash);
144  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
145  Scalar rho = FluidSystem::density(fsFlash, paramCache, phaseIdx);
146  fsFlash.setDensity(phaseIdx, rho);
147  }
148  // calculate the "global molarities"
149  ComponentVector globalMolarities(0.0);
150  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
151  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
152  globalMolarities[compIdx] +=
153  fsFlash.saturation(phaseIdx) * fsFlash.molarity(phaseIdx, compIdx);
154  }
155  }
156 
157  // run the flash calculation
158  NcpFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
159 
160  // use the result to assign the primary variables
161  assignNaive(fsFlash);
162  }
163 
168  short phasePresence() const
169  { return phasePresence_; }
170 
177  void setPhasePresence(short value)
178  { phasePresence_ = value; }
179 
187  void setPhasePresent(unsigned phaseIdx, bool yesno = true)
188  {
189  if (yesno)
190  setPhasePresence(phasePresence_ | (1 << phaseIdx));
191  else
192  setPhasePresence(phasePresence_& ~(1 << phaseIdx));
193  }
194 
199  unsigned implicitSaturationIdx() const
200  { return lowestPresentPhaseIdx(); }
201 
210  static bool phaseIsPresent(unsigned phaseIdx, short phasePresence)
211  { return phasePresence& (1 << phaseIdx); }
212 
219  bool phaseIsPresent(unsigned phaseIdx) const
220  { return phasePresence_& (1 << phaseIdx); }
221 
225  unsigned lowestPresentPhaseIdx() const
226  { return static_cast<unsigned>(ffs(phasePresence_) - 1); }
227 
231  ThisType& operator=(const Implementation& value)
232  {
233  ParentType::operator=(value);
234  phasePresence_ = value.phasePresence_;
235 
236  return *this;
237  }
238 
242  ThisType& operator=(Scalar value)
243  {
244  ParentType::operator=(value);
245 
246  phasePresence_ = 0;
247  return *this;
248  }
249 
257  Evaluation explicitSaturationValue(unsigned phaseIdx, unsigned timeIdx) const
258  {
259  if (!phaseIsPresent(phaseIdx) || phaseIdx == lowestPresentPhaseIdx())
260  // non-present phases have saturation 0
261  return 0.0;
262 
263  unsigned varIdx = switch0Idx + phaseIdx - 1;
264  if (timeIdx != 0)
265  Toolbox::createConstant((*this)[varIdx]);
266  return Toolbox::createVariable((*this)[varIdx], varIdx);
267  }
268 
272  template <class FluidState>
273  void assignNaive(const FluidState& fluidState)
274  {
275  typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
276 
277  // assign the phase temperatures. this is out-sourced to
278  // the energy module
279  EnergyModule::setPriVarTemperatures(*this, fluidState);
280 
281  // set the pressure of the first phase
282  (*this)[pressure0Idx] = FsToolbox::value(fluidState.pressure(/*phaseIdx=*/0));
283  Opm::Valgrind::CheckDefined((*this)[pressure0Idx]);
284 
285  // determine the phase presence.
286  phasePresence_ = 0;
287  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
288  // use a NCP condition to determine if the phase is
289  // present or not
290  Scalar a = 1;
291  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
292  a -= FsToolbox::value(fluidState.moleFraction(phaseIdx, compIdx));
293  }
294  Scalar b = FsToolbox::value(fluidState.saturation(phaseIdx));
295 
296  if (b > a)
297  phasePresence_ |= (1 << phaseIdx);
298  }
299 
300  // some phase must be present
301  if (phasePresence_ == 0)
302  OPM_THROW(Opm::NumericalProblem,
303  "Phase state was 0, i.e., no fluid is present");
304 
305  // set the primary variables which correspond to mole
306  // fractions of the present phase which has the lowest index.
307  unsigned lowestPhaseIdx = lowestPresentPhaseIdx();
308  for (unsigned switchIdx = 0; switchIdx < numPhases - 1; ++switchIdx) {
309  unsigned phaseIdx = switchIdx;
310  unsigned compIdx = switchIdx + 1;
311  if (switchIdx >= lowestPhaseIdx)
312  ++phaseIdx;
313 
314  if (phaseIsPresent(phaseIdx)) {
315  (*this)[switch0Idx + switchIdx] = FsToolbox::value(fluidState.saturation(phaseIdx));
316  Opm::Valgrind::CheckDefined((*this)[switch0Idx + switchIdx]);
317  }
318  else {
319  (*this)[switch0Idx + switchIdx] =
320  FsToolbox::value(fluidState.moleFraction(lowestPhaseIdx, compIdx));
321  Opm::Valgrind::CheckDefined((*this)[switch0Idx + switchIdx]);
322  }
323  }
324 
325  // set the mole fractions in of the remaining components in
326  // the phase with the lowest index
327  for (unsigned compIdx = numPhases - 1; compIdx < numComponents - 1; ++compIdx) {
328  (*this)[switch0Idx + compIdx] =
329  FsToolbox::value(fluidState.moleFraction(lowestPhaseIdx, compIdx + 1));
330  Opm::Valgrind::CheckDefined((*this)[switch0Idx + compIdx]);
331  }
332  }
333 
337  void print(std::ostream& os = std::cout) const
338  {
339  os << "(p_" << FluidSystem::phaseName(0) << " = "
340  << this->operator[](pressure0Idx);
341  unsigned lowestPhaseIdx = lowestPresentPhaseIdx();
342  for (unsigned switchIdx = 0; switchIdx < numPhases - 1; ++switchIdx) {
343  unsigned phaseIdx = switchIdx;
344  unsigned compIdx = switchIdx + 1;
345  if (phaseIdx >= lowestPhaseIdx)
346  ++phaseIdx; // skip the saturation of the present
347  // phase with the lowest index
348 
349  if (phaseIsPresent(phaseIdx)) {
350  os << ", S_" << FluidSystem::phaseName(phaseIdx) << " = "
351  << (*this)[switch0Idx + switchIdx];
352  }
353  else {
354  os << ", x_" << FluidSystem::phaseName(lowestPhaseIdx) << "^"
355  << FluidSystem::componentName(compIdx) << " = "
356  << (*this)[switch0Idx + switchIdx];
357  }
358  }
359  for (unsigned compIdx = numPhases - 1; compIdx < numComponents - 1;
360  ++compIdx) {
361  os << ", x_" << FluidSystem::phaseName(lowestPhaseIdx) << "^"
362  << FluidSystem::componentName(compIdx + 1) << " = "
363  << (*this)[switch0Idx + compIdx];
364  }
365  os << ")";
366  os << ", phase presence: " << static_cast<int>(phasePresence_) << std::flush;
367  }
368 
369 private:
370  short phasePresence_;
371 };
372 
373 } // namespace Ewoms
374 
375 #endif
Evaluation explicitSaturationValue(unsigned phaseIdx, unsigned timeIdx) const
Returns an explcitly stored saturation for a given phase.
Definition: pvsprimaryvariables.hh:257
Represents the primary variables used in the primary variable switching compositional model...
Definition: pvsprimaryvariables.hh:59
Definition: baseauxiliarymodule.hh:37
short phasePresence() const
Return the fluid phases which are present in a given control volume.
Definition: pvsprimaryvariables.hh:168
Provides the auxiliary methods required for consideration of the energy equation. ...
Definition: energymodule.hh:59
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
void print(std::ostream &os=std::cout) const
Prints the names of the primary variables and their values.
Definition: pvsprimaryvariables.hh:337
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
unsigned implicitSaturationIdx() const
Returns the index of the phase with&#39;s its saturation is determined by the closure condition of satura...
Definition: pvsprimaryvariables.hh:199
unsigned lowestPresentPhaseIdx() const
Returns the phase with the lowest index that is present.
Definition: pvsprimaryvariables.hh:225
static bool phaseIsPresent(unsigned phaseIdx, short phasePresence)
Returns true iff a phase is present for a given phase presence.
Definition: pvsprimaryvariables.hh:210
The indices for the compositional multi-phase primary variable switching model.
ThisType & operator=(Scalar value)
Assignment operator from a scalar value.
Definition: pvsprimaryvariables.hh:242
void setPhasePresence(short value)
Set which fluid phases are present in a given control volume.
Definition: pvsprimaryvariables.hh:177
PvsPrimaryVariables(const PvsPrimaryVariables &value)
Default constructor.
Definition: pvsprimaryvariables.hh:104
Represents the primary variables used by the a model.
Definition: fvbaseprimaryvariables.hh:48
FvBasePrimaryVariables & operator=(const FvBasePrimaryVariables &value)=default
Assignment from another primary variables object.
bool phaseIsPresent(unsigned phaseIdx) const
Returns true iff a phase is present for the current phase presence.
Definition: pvsprimaryvariables.hh:219
void setPhasePresent(unsigned phaseIdx, bool yesno=true)
Set whether a given indivividual phase should be present or not.
Definition: pvsprimaryvariables.hh:187
Represents the primary variables used by the a model.
Declares the properties required for the compositional multi-phase primary variable switching 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:115
ThisType & operator=(const Implementation &value)
Assignment operator from an other primary variables object.
Definition: pvsprimaryvariables.hh:231
PvsPrimaryVariables(Scalar value)
Constructor with assignment from scalar.
Definition: pvsprimaryvariables.hh:92
void assignNaive(const FluidState &fluidState)
Directly retrieve the primary variables from an arbitrary fluid state.
Definition: pvsprimaryvariables.hh:273