27#ifndef OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
28#define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
35#include <dune/common/fvector.hh>
36#include <dune/common/fmatrix.hh>
37#include <dune/common/version.hh>
48template <
class Scalar>
56 : phaseIdx_(phaseIndex)
64 void set(
unsigned phaseIndex,
unsigned compIndex, Scalar val)
66 phaseIdx_ = phaseIndex;
121template <
class Scalar,
class Flu
idSystem,
class Evaluation = Scalar>
124 static const int numPhases = FluidSystem::numPhases;
125 static const int numComponents = FluidSystem::numComponents;
129 static_assert(numPhases <= numComponents,
130 "This solver requires that the number fluid phases is smaller or equal "
131 "to the number of components");
158 template <
class Flu
idState,
class ParameterCache>
159 static void solve(FluidState& fluidState,
160 ParameterCache& paramCache,
163 unsigned numAuxConstraints,
165 bool setInternalEnergy)
167 static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
168 "The scalar type of the fluid state must be 'Evaluation'");
175 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
176 assert(FluidSystem::isIdealMixture(phaseIdx));
181 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
182 paramCache.updatePhase(fluidState, phaseIdx);
187 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
188 Evaluation fugCoeff = decay<Evaluation>(
189 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
190 fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
196 static const int numEq = numComponents*numPhases;
197 Dune::FieldMatrix<Evaluation, numEq, numEq> M(0.0);
198 Dune::FieldVector<Evaluation, numEq> x(0.0);
199 Dune::FieldVector<Evaluation, numEq> b(0.0);
203 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
204 const Evaluation& entryCol1 =
205 fluidState.fugacityCoefficient(0, compIdx)
206 *fluidState.pressure(0);
207 unsigned col1Idx = compIdx;
209 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
210 unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
211 unsigned col2Idx = phaseIdx*numComponents + compIdx;
213 const Evaluation& entryCol2 =
214 fluidState.fugacityCoefficient(phaseIdx, compIdx)
215 *fluidState.pressure(phaseIdx);
217 M[rowIdx][col1Idx] = entryCol1;
218 M[rowIdx][col2Idx] = -entryCol2;
225 unsigned presentPhases = 0;
226 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
227 if (!(phasePresence& (1 << phaseIdx)))
230 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
234 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
235 unsigned colIdx = phaseIdx*numComponents + compIdx;
237 M[rowIdx][colIdx] = 1.0;
241 assert(presentPhases + numAuxConstraints == numComponents);
244 for (
unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
245 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
246 b[rowIdx] = auxConstraints[auxEqIdx].
value();
248 unsigned colIdx = auxConstraints[auxEqIdx].
phaseIdx()*numComponents + auxConstraints[auxEqIdx].
compIdx();
249 M[rowIdx][colIdx] = 1.0;
254#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
256 Dune::FMatrixPrecision<Scalar>::set_singular_limit(eps);
260 catch (
const Dune::FMatrixError& e) {
261 std::ostringstream oss;
262 oss <<
"Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() <<
"; M="<<M;
272 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
273 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
274 unsigned rowIdx = phaseIdx*numComponents + compIdx;
275 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
277 paramCache.updateComposition(fluidState, phaseIdx);
279 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
280 fluidState.setDensity(phaseIdx, rho);
283 const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
284 fluidState.setViscosity(phaseIdx, mu);
287 if (setInternalEnergy) {
288 const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
289 fluidState.setEnthalpy(phaseIdx, h);
301 template <
class Flu
idState,
class ParameterCache>
302 static void solve(FluidState& fluidState,
303 ParameterCache& paramCache,
305 bool setInternalEnergy)
Provides the opm-material specific exception classes.
Some templates to wrap the valgrind client request macros.
Specifies an auxiliary constraint for the MiscibleMultiPhaseComposition constraint solver.
Definition: MiscibleMultiPhaseComposition.hpp:50
Scalar value() const
Returns value of the mole fraction of the auxiliary constraint.
Definition: MiscibleMultiPhaseComposition.hpp:89
MMPCAuxConstraint(unsigned phaseIndex, unsigned compIndex, Scalar val)
Definition: MiscibleMultiPhaseComposition.hpp:55
unsigned compIdx() const
Returns the index of the component for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:82
void set(unsigned phaseIndex, unsigned compIndex, Scalar val)
Specify the auxiliary constraint.
Definition: MiscibleMultiPhaseComposition.hpp:64
unsigned phaseIdx() const
Returns the index of the fluid phase for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:75
MMPCAuxConstraint()
Definition: MiscibleMultiPhaseComposition.hpp:52
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:123
static void solve(FluidState &fluidState, ParameterCache ¶mCache, bool setViscosity, bool setInternalEnergy)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:302
static void solve(FluidState &fluidState, ParameterCache ¶mCache, int phasePresence, const MMPCAuxConstraint< Evaluation > *auxConstraints, unsigned numAuxConstraints, bool setViscosity, bool setInternalEnergy)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:159
Definition: Exceptions.hpp:46
Definition: Air_Mesitylene.hpp:34
ReturnEval_< Evaluation1, Evaluation2 >::type min(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:346