ncpprimaryvariables.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_NCP_PRIMARY_VARIABLES_HH
29 #define EWOMS_NCP_PRIMARY_VARIABLES_HH
30 
31 #include "ncpproperties.hh"
32 
35 
36 #include <opm/material/constraintsolvers/NcpFlash.hpp>
37 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
38 #include <opm/material/densead/Math.hpp>
39 
40 #include <dune/common/fvector.hh>
41 
42 namespace Ewoms {
43 
53 template <class TypeTag>
55 {
56  typedef FvBasePrimaryVariables<TypeTag> ParentType;
57 
58  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
59  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
60  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
61  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
62  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
63 
64  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
65  enum { pressure0Idx = Indices::pressure0Idx };
66  enum { saturation0Idx = Indices::saturation0Idx };
67  enum { fugacity0Idx = Indices::fugacity0Idx };
68 
69  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
70  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
71  typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
72 
73  enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
75 
76  typedef Opm::NcpFlash<Scalar, FluidSystem> NcpFlash;
77  typedef Opm::MathToolbox<Evaluation> Toolbox;
78 
79 public:
80  NcpPrimaryVariables() : ParentType()
81  {}
82 
86  NcpPrimaryVariables(Scalar value) : ParentType(value)
87  {}
88 
93  NcpPrimaryVariables(const NcpPrimaryVariables& value) = default;
94  NcpPrimaryVariables& operator=(const NcpPrimaryVariables& value) = default;
95 
99  template <class FluidState>
100  void assignMassConservative(const FluidState& fluidState,
101  const MaterialLawParams& matParams,
102  bool isInEquilibrium = false)
103  {
104  typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
105 
106 #ifndef NDEBUG
107  // make sure the temperature is the same in all fluid phases
108  for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
109  assert(fluidState.temperature(0) == fluidState.temperature(phaseIdx));
110  }
111 #endif // NDEBUG
112 
113  // for the equilibrium case, we don't need complicated
114  // computations.
115  if (isInEquilibrium) {
116  assignNaive(fluidState);
117  return;
118  }
119 
120  // use a flash calculation to calculate a fluid state in
121  // thermodynamic equilibrium
122  typename FluidSystem::template ParameterCache<Scalar> paramCache;
123  Opm::CompositionalFluidState<Scalar, FluidSystem> fsFlash;
124 
125  // use the externally given fluid state as initial value for
126  // the flash calculation
127  fsFlash.assign(fluidState);
128 
129  // calculate the phase densities
130  paramCache.updateAll(fsFlash);
131  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
132  Scalar rho = FluidSystem::density(fsFlash, paramCache, phaseIdx);
133  fsFlash.setDensity(phaseIdx, rho);
134  }
135 
136  // calculate the "global molarities"
137  ComponentVector globalMolarities(0.0);
138  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
139  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
140  globalMolarities[compIdx] +=
141  FsToolbox::value(fsFlash.saturation(phaseIdx))
142  * FsToolbox::value(fsFlash.molarity(phaseIdx, compIdx));
143  }
144  }
145 
146  // run the flash calculation
147  NcpFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
148 
149  // use the result to assign the primary variables
150  assignNaive(fsFlash);
151  }
152 
156  template <class FluidState>
157  void assignNaive(const FluidState& fluidState, unsigned refPhaseIdx = 0)
158  {
159  typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
160 
161  // assign the phase temperatures. this is out-sourced to
162  // the energy module
163  EnergyModule::setPriVarTemperatures(*this, fluidState);
164 
165  // assign fugacities.
166  typename FluidSystem::template ParameterCache<Scalar> paramCache;
167  paramCache.updatePhase(fluidState, refPhaseIdx);
168  Scalar pRef = FsToolbox::value(fluidState.pressure(refPhaseIdx));
169  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
170  // we always compute the fugacities because they are quite exotic quantities
171  // and this easily forgotten to be specified
172  Scalar fugCoeff =
173  FluidSystem::template fugacityCoefficient<FluidState, Scalar>(fluidState,
174  paramCache,
175  refPhaseIdx,
176  compIdx);
177  (*this)[fugacity0Idx + compIdx] =
178  fugCoeff*fluidState.moleFraction(refPhaseIdx, compIdx)*pRef;
179  }
180 
181  // assign pressure of first phase
182  (*this)[pressure0Idx] = FsToolbox::value(fluidState.pressure(/*phaseIdx=*/0));
183 
184  // assign first M - 1 saturations
185  for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
186  (*this)[saturation0Idx + phaseIdx] = FsToolbox::value(fluidState.saturation(phaseIdx));
187  }
188 };
189 
190 } // namespace Ewoms
191 
192 #endif
void assignNaive(const FluidState &fluidState, unsigned refPhaseIdx=0)
Directly retrieve the primary variables from an arbitrary fluid state.
Definition: ncpprimaryvariables.hh:157
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: ncpprimaryvariables.hh:100
Definition: baseauxiliarymodule.hh:37
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
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
NcpPrimaryVariables(Scalar value)
Constructor with assignment from scalar.
Definition: ncpprimaryvariables.hh:86
Represents the primary variables used by the compositional multi-phase NCP model. ...
Definition: ncpprimaryvariables.hh:54
Declares the properties required for the NCP compositional multi-phase model.
Represents the primary variables used by the a model.
Definition: fvbaseprimaryvariables.hh:48
Represents the primary variables used by the a model.