My Project
Loading...
Searching...
No Matches
MiscibleMultiPhaseComposition.hpp
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*/
27#ifndef OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
28#define OPM_MISCIBLE_MULTIPHASE_COMPOSITION_HPP
29
31
33
35
36#include <dune/common/fvector.hh>
37#include <dune/common/fmatrix.hh>
38
39namespace Opm {
40
48template <class Scalar>
50{
51public:
53 {}
54
55 MMPCAuxConstraint(unsigned phaseIndex, unsigned compIndex, Scalar val)
56 : phaseIdx_(phaseIndex)
57 , compIdx_(compIndex)
58 , value_(val)
59 {}
60
64 void set(unsigned phaseIndex, unsigned compIndex, Scalar val)
65 {
66 phaseIdx_ = phaseIndex;
67 compIdx_ = compIndex;
68 value_ = val;
69 }
70
75 unsigned phaseIdx() const
76 { return phaseIdx_; }
77
82 unsigned compIdx() const
83 { return compIdx_; }
84
89 Scalar value() const
90 { return value_; }
91
92private:
93 unsigned phaseIdx_;
94 unsigned compIdx_;
95 Scalar value_;
96};
97
121template <class Scalar, class FluidSystem, class Evaluation = Scalar>
123{
124 static const int numPhases = FluidSystem::numPhases;
125 static const int numComponents = FluidSystem::numComponents;
126
128
129 static_assert(numPhases <= numComponents,
130 "This solver requires that the number fluid phases is smaller or equal "
131 "to the number of components");
132
133
134public:
158 template <class FluidState, class ParameterCache>
159 static void solve(FluidState& fluidState,
160 ParameterCache& paramCache,
161 int phasePresence,
162 const MMPCAuxConstraint<Evaluation>* auxConstraints,
163 unsigned numAuxConstraints,
164 bool setViscosity,
165 bool setInternalEnergy)
166 {
167 static_assert(std::is_same<typename FluidState::Scalar, Evaluation>::value,
168 "The scalar type of the fluid state must be 'Evaluation'");
169
170#ifndef NDEBUG
171 // currently this solver can only handle fluid systems which
172 // assume ideal mixtures of all fluids. TODO: relax this
173 // (requires solving a non-linear system of equations, i.e. using
174 // newton method.)
175 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
176 assert(FluidSystem::isIdealMixture(phaseIdx));
177 }
178#endif
179
180 // compute all fugacity coefficients
181 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
182 paramCache.updatePhase(fluidState, phaseIdx);
183
184 // since we assume ideal mixtures, the fugacity
185 // coefficients of the components cannot depend on
186 // composition, i.e. the parameters in the cache are valid
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);
191 }
192 }
193
194 // create the linear system of equations which defines the
195 // mole fractions
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);
200
201 // assemble the equations expressing the fact that the
202 // fugacities of each component are equal in all phases
203 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
204 const Evaluation& entryCol1 =
205 fluidState.fugacityCoefficient(/*phaseIdx=*/0, compIdx)
206 *fluidState.pressure(/*phaseIdx=*/0);
207 unsigned col1Idx = compIdx;
208
209 for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
210 unsigned rowIdx = (phaseIdx - 1)*numComponents + compIdx;
211 unsigned col2Idx = phaseIdx*numComponents + compIdx;
212
213 const Evaluation& entryCol2 =
214 fluidState.fugacityCoefficient(phaseIdx, compIdx)
215 *fluidState.pressure(phaseIdx);
216
217 M[rowIdx][col1Idx] = entryCol1;
218 M[rowIdx][col2Idx] = -entryCol2;
219 }
220 }
221
222 // assemble the equations expressing the assumption that the
223 // sum of all mole fractions in each phase must be 1 for the
224 // phases present.
225 unsigned presentPhases = 0;
226 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
227 if (!(phasePresence& (1 << phaseIdx)))
228 continue;
229
230 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases;
231 presentPhases += 1;
232
233 b[rowIdx] = 1.0;
234 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
235 unsigned colIdx = phaseIdx*numComponents + compIdx;
236
237 M[rowIdx][colIdx] = 1.0;
238 }
239 }
240
241 assert(presentPhases + numAuxConstraints == numComponents);
242
243 // incorperate the auxiliary equations, i.e., the explicitly given mole fractions
244 for (unsigned auxEqIdx = 0; auxEqIdx < numAuxConstraints; ++auxEqIdx) {
245 unsigned rowIdx = numComponents*(numPhases - 1) + presentPhases + auxEqIdx;
246 b[rowIdx] = auxConstraints[auxEqIdx].value();
247
248 unsigned colIdx = auxConstraints[auxEqIdx].phaseIdx()*numComponents + auxConstraints[auxEqIdx].compIdx();
249 M[rowIdx][colIdx] = 1.0;
250 }
251
252 // solve for all mole fractions
253 try {
254 M.solve(x, b);
255 }
256 catch (const Dune::FMatrixError& e) {
257 std::ostringstream oss;
258 oss << "Numerical problem in MiscibleMultiPhaseComposition::solve(): " << e.what() << "; M="<<M;
259 throw NumericalProblem(oss.str());
260 }
261 catch (...) {
262 throw;
263 }
264
265
266 // set all mole fractions and the additional quantities in
267 // the fluid state
268 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
269 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
270 unsigned rowIdx = phaseIdx*numComponents + compIdx;
271 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
272 }
273 paramCache.updateComposition(fluidState, phaseIdx);
274
275 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
276 fluidState.setDensity(phaseIdx, rho);
277
278 if (setViscosity) {
279 const Evaluation& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
280 fluidState.setViscosity(phaseIdx, mu);
281 }
282
283 if (setInternalEnergy) {
284 const Evaluation& h = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
285 fluidState.setEnthalpy(phaseIdx, h);
286 }
287 }
288 }
289
297 template <class FluidState, class ParameterCache>
298 static void solve(FluidState& fluidState,
299 ParameterCache& paramCache,
300 bool setViscosity,
301 bool setInternalEnergy)
302 {
303 solve(fluidState,
304 paramCache,
305 /*phasePresence=*/0xffffff,
306 /*numAuxConstraints=*/0,
307 /*auxConstraints=*/0,
308 setViscosity,
309 setInternalEnergy);
310 }
311};
312
313} // namespace Opm
314
315#endif
Provides the OPM specific exception classes.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
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
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
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 &paramCache, 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:298
static void solve(FluidState &fluidState, ParameterCache &paramCache, 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:40
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30
Definition MathToolbox.hpp:50