27 #ifndef OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
28 #define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
32 #include <dune/common/fvector.hh>
33 #include <dune/common/fmatrix.hh>
35 #include <opm/common/Exceptions.hpp>
36 #include <opm/common/ErrorMacros.hpp>
48 template <
class Scalar>
121 template <
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 ¶mCache,
163 unsigned numAuxConstraints,
165 bool setInternalEnergy)
168 static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
169 "The scalar type of the fluid state must be 'Evaluation'");
176 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
177 assert(FluidSystem::isIdealMixture(phaseIdx));
182 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
183 paramCache.updatePhase(fluidState, phaseIdx);
188 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
189 Evaluation fugCoeff = FsToolbox::template toLhs<Evaluation>(
190 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
191 fluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
197 static const int numEq = numComponents*numPhases;
198 Dune::FieldMatrix<Evaluation, numEq, numEq> M(Toolbox::createConstant(0.0));
199 Dune::FieldVector<Evaluation, numEq> x(Toolbox::createConstant(0.0));
200 Dune::FieldVector<Evaluation, numEq> b(Toolbox::createConstant(0.0));
204 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
205 const Evaluation& entryCol1 =
206 fluidState.fugacityCoefficient(0, compIdx)
207 *fluidState.pressure(0);
208 unsigned col1Idx = compIdx;
210 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
211 unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
212 unsigned col2Idx = phaseIdx*numComponents + compIdx;
214 const Evaluation& entryCol2 =
215 fluidState.fugacityCoefficient(phaseIdx, compIdx)
216 *fluidState.pressure(phaseIdx);
218 M[rowIdx][col1Idx] = entryCol1;
219 M[rowIdx][col2Idx] = -entryCol2;
226 unsigned presentPhases = 0;
227 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
228 if (!(phasePresence & (1 << phaseIdx)))
231 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
234 b[rowIdx] = Toolbox::createConstant(1.0);
235 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
236 unsigned colIdx = phaseIdx*numComponents + compIdx;
238 M[rowIdx][colIdx] = Toolbox::createConstant(1.0);
242 assert(presentPhases + numAuxConstraints == numComponents);
245 for (
unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
246 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
247 b[rowIdx] = auxConstraints[auxEqIdx].
value();
249 unsigned colIdx = auxConstraints[auxEqIdx].
phaseIdx()*numComponents + auxConstraints[auxEqIdx].
compIdx();
250 M[rowIdx][colIdx] = 1.0;
255 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-50);
258 catch (
const Dune::FMatrixError &e) {
259 OPM_THROW(NumericalProblem,
260 "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() <<
"; M="<<M);
269 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
270 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
271 unsigned rowIdx = phaseIdx*numComponents + compIdx;
272 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
274 paramCache.updateComposition(fluidState, phaseIdx);
276 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
277 fluidState.setDensity(phaseIdx, rho);
280 const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
281 fluidState.setViscosity(phaseIdx, mu);
284 if (setInternalEnergy) {
285 const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
286 fluidState.setEnthalpy(phaseIdx, h);
298 template <
class Flu
idState,
class ParameterCache>
299 static void solve(FluidState &fluidState,
300 ParameterCache ¶mCache,
302 bool setInternalEnergy)
Definition: Air_Mesitylene.hpp:31
Some templates to wrap the valgrind client request macros.
void set(unsigned phaseIdx, unsigned compIdx, Scalar value)
Specify the auxiliary constraint.
Definition: MiscibleMultiPhaseComposition.hpp:64
MMPCAuxConstraint()
Definition: MiscibleMultiPhaseComposition.hpp:52
Specifies an auxiliary constraint for the MiscibleMultiPhaseComposition constraint solver...
Definition: MiscibleMultiPhaseComposition.hpp:49
unsigned phaseIdx() const
Returns the index of the fluid phase for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:75
Computes the composition of all phases of a N-phase, N-component fluid system assuming that all N pha...
Definition: MiscibleMultiPhaseComposition.hpp:122
unsigned compIdx() const
Returns the index of the component for which the auxiliary constraint is specified.
Definition: MiscibleMultiPhaseComposition.hpp:82
Scalar value() const
Returns value of the mole fraction of the auxiliary constraint.
Definition: MiscibleMultiPhaseComposition.hpp:89
MMPCAuxConstraint(unsigned phaseIdx, unsigned compIdx, Scalar value)
Definition: MiscibleMultiPhaseComposition.hpp:55