27 #ifndef OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
28 #define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
33 #include <opm/common/utility/platform_dependent/disable_warnings.h>
35 #include <dune/common/fvector.hh>
36 #include <dune/common/fmatrix.hh>
38 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
41 #include <opm/common/Exceptions.hpp>
42 #include <opm/common/ErrorMacros.hpp>
43 #include <opm/common/Valgrind.hpp>
54 template <
class Scalar>
62 : phaseIdx_(phaseIndex)
70 void set(
unsigned phaseIndex,
unsigned compIndex, Scalar val)
72 phaseIdx_ = phaseIndex;
127 template <
class Scalar,
class Flu
idSystem,
class Evaluation = Scalar>
130 static const int numPhases = FluidSystem::numPhases;
131 static const int numComponents = FluidSystem::numComponents;
135 static_assert(numPhases <= numComponents,
136 "This solver requires that the number fluid phases is smaller or equal "
137 "to the number of components");
164 template <
class Flu
idState,
class ParameterCache>
165 static void solve(FluidState& fluidState,
166 ParameterCache& paramCache,
169 unsigned numAuxConstraints,
171 bool setInternalEnergy)
173 static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
174 "The scalar type of the fluid state must be 'Evaluation'");
181 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
182 assert(FluidSystem::isIdealMixture(phaseIdx));
187 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
188 paramCache.updatePhase(fluidState, phaseIdx);
193 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
194 Evaluation fugCoeff = Opm::decay<Evaluation>(
195 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
196 fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
202 static const int numEq = numComponents*numPhases;
203 Dune::FieldMatrix<Evaluation, numEq, numEq> M(0.0);
204 Dune::FieldVector<Evaluation, numEq> x(0.0);
205 Dune::FieldVector<Evaluation, numEq> b(0.0);
209 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
210 const Evaluation& entryCol1 =
211 fluidState.fugacityCoefficient(0, compIdx)
212 *fluidState.pressure(0);
213 unsigned col1Idx = compIdx;
215 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
216 unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
217 unsigned col2Idx = phaseIdx*numComponents + compIdx;
219 const Evaluation& entryCol2 =
220 fluidState.fugacityCoefficient(phaseIdx, compIdx)
221 *fluidState.pressure(phaseIdx);
223 M[rowIdx][col1Idx] = entryCol1;
224 M[rowIdx][col2Idx] = -entryCol2;
231 unsigned presentPhases = 0;
232 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
233 if (!(phasePresence& (1 << phaseIdx)))
236 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
240 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
241 unsigned colIdx = phaseIdx*numComponents + compIdx;
243 M[rowIdx][colIdx] = 1.0;
247 assert(presentPhases + numAuxConstraints == numComponents);
250 for (
unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
251 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
252 b[rowIdx] = auxConstraints[auxEqIdx].
value();
254 unsigned colIdx = auxConstraints[auxEqIdx].
phaseIdx()*numComponents + auxConstraints[auxEqIdx].
compIdx();
255 M[rowIdx][colIdx] = 1.0;
260 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-50);
263 catch (
const Dune::FMatrixError& e) {
264 OPM_THROW(NumericalProblem,
265 "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() <<
"; M="<<M);
274 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
275 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
276 unsigned rowIdx = phaseIdx*numComponents + compIdx;
277 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
279 paramCache.updateComposition(fluidState, phaseIdx);
281 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
282 fluidState.setDensity(phaseIdx, rho);
285 const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
286 fluidState.setViscosity(phaseIdx, mu);
289 if (setInternalEnergy) {
290 const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
291 fluidState.setEnthalpy(phaseIdx, h);
303 template <
class Flu
idState,
class ParameterCache>
304 static void solve(FluidState& fluidState,
305 ParameterCache& paramCache,
307 bool setInternalEnergy)
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:128
unsigned phaseIdx() const
Returns the index of the fluid phase for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:81
void set(unsigned phaseIndex, unsigned compIndex, Scalar val)
Specify the auxiliary constraint.
Definition: MiscibleMultiPhaseComposition.hpp:70
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:165
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:304
Specifies an auxiliary constraint for the MiscibleMultiPhaseComposition constraint solver...
Definition: MiscibleMultiPhaseComposition.hpp:55
unsigned compIdx() const
Returns the index of the component for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:88
Scalar value() const
Returns value of the mole fraction of the auxiliary constraint.
Definition: MiscibleMultiPhaseComposition.hpp:95