20 #ifndef OPM_BLACKOILFLUID_HEADER_INCLUDED
21 #define OPM_BLACKOILFLUID_HEADER_INCLUDED
24 #include <opm/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp>
25 #include <opm/porsol/blackoil/fluid/FluidStateBlackoil.hpp>
26 #include <opm/porsol/blackoil/fluid/BlackoilPVT.hpp>
27 #include <dune/common/fvector.hh>
36 class BlackoilFluidData;
45 typedef BlackoilFluidData FluidData;
47 void init(
const Opm::Deck& deck)
49 fmi_params_.init(deck);
52 const auto& densityRecord = deck.getKeyword(
"DENSITY").getRecord(0);
53 surface_densities_[Oil] = densityRecord.getItem(
"OIL").getSIDouble(0);
54 surface_densities_[Water] = densityRecord.getItem(
"WATER").getSIDouble(0);
55 surface_densities_[Gas] = densityRecord.getItem(
"GAS").getSIDouble(0);
58 FluidState computeState(PhaseVec phase_pressure, CompVec z)
const
61 state.temperature_ = 300;
62 state.phase_pressure_ = phase_pressure;
63 state.surface_volume_ = z;
65 computeEquilibrium(state);
68 for (
int phase = 0; phase < numPhases; ++phase) {
69 state.mobility_[phase] = state.relperm_[phase]/state.viscosity_[phase];
70 for (
int p2 = 0; p2 < numPhases; ++p2) {
72 state.dmobility_[phase][p2] = state.drelperm_[phase][p2]/state.viscosity_[phase];
79 const CompVec& surfaceDensities()
const
81 return surface_densities_;
87 PhaseVec phase_dens(0.0);
88 for (
int phase = 0; phase < numPhases; ++phase) {
89 for (
int comp = 0; comp < numComponents; ++comp) {
90 phase_dens[phase] += A[numComponents*phase + comp]*surface_densities_[comp];
101 template <
class States>
104 const std::vector<PhaseVec>& p = states.phase_pressure;
105 const std::vector<CompVec>& z = states.surface_volume_density;
106 std::vector<PhaseVec>& B = states.formation_volume_factor;
107 std::vector<PhaseVec>& R = states.solution_factor;
114 template <
class States>
118 const std::vector<PhaseVec>& p = states.phase_pressure;
119 const std::vector<CompVec>& z = states.surface_volume_density;
120 std::vector<PhaseVec>& mu = states.viscosity;
121 pvt_.getViscosity(p, z, mu);
126 template <
class States>
129 const std::vector<PhaseVec>& p = states.phase_pressure;
130 const std::vector<CompVec>& z = states.surface_volume_density;
131 std::vector<PhaseVec>& B = states.formation_volume_factor;
132 std::vector<PhaseVec>& dB = states.formation_volume_factor_deriv;
133 std::vector<PhaseVec>& R = states.solution_factor;
134 std::vector<PhaseVec>& dR = states.solution_factor_deriv;
135 std::vector<PhaseVec>& mu = states.viscosity;
136 pvt_.dBdp(p, z, B, dB);
137 pvt_.dRdp(p, z, R, dR);
138 pvt_.getViscosity(p, z, mu);
143 template <
class States>
146 int num = states.formation_volume_factor.size();
147 states.state_matrix.resize(num);
148 #pragma omp parallel for
149 for (
int i = 0; i < num; ++i) {
150 const PhaseVec& B = states.formation_volume_factor[i];
151 const PhaseVec& R = states.solution_factor[i];
152 PhaseToCompMatrix& At = states.state_matrix[i];
157 At[Aqua][Water] = 1.0/B[Aqua];
158 At[Vapour][Gas] = 1.0/B[Vapour];
159 At[Liquid][Gas] = R[Liquid]/B[Liquid];
160 At[Vapour][Oil] = R[Vapour]/B[Vapour];
161 At[Liquid][Oil] = 1.0/B[Liquid];
168 template <
class States>
171 int num = states.formation_volume_factor.size();
172 states.state_matrix.resize(num);
173 states.phase_volume_density.resize(num);
174 states.total_phase_volume_density.resize(num);
175 states.saturation.resize(num);
176 states.phase_compressibility.resize(num);
177 states.total_compressibility.resize(num);
178 states.experimental_term.resize(num);
179 #pragma omp parallel for
180 for (
int i = 0; i < num; ++i) {
181 const CompVec& z = states.surface_volume_density[i];
182 const PhaseVec& B = states.formation_volume_factor[i];
183 const PhaseVec& dB = states.formation_volume_factor_deriv[i];
184 const PhaseVec& R = states.solution_factor[i];
185 const PhaseVec& dR = states.solution_factor_deriv[i];
186 PhaseToCompMatrix& At = states.state_matrix[i];
187 PhaseVec& u = states.phase_volume_density[i];
188 double& tot_phase_vol_dens = states.total_phase_volume_density[i];
189 PhaseVec& s = states.saturation[i];
190 PhaseVec& cp = states.phase_compressibility[i];
191 double& tot_comp = states.total_compressibility[i];
192 double& exp_term = states.experimental_term[i];
193 computeSingleEquilibrium(B, dB, R, dR, z,
194 At, u, tot_phase_vol_dens,
195 s, cp, tot_comp, exp_term);
201 template <
class States>
204 int num = states.saturation.size();
205 states.relperm.resize(num);
206 states.mobility.resize(num);
207 #pragma omp parallel for
208 for (
int i = 0; i < num; ++i) {
209 const CompVec& s = states.saturation[i];
210 const PhaseVec& mu = states.viscosity[i];
211 PhaseVec& kr = states.relperm[i];
212 PhaseVec& lambda = states.mobility[i];
214 for (
int phase = 0; phase < numPhases; ++phase) {
215 lambda[phase] = kr[phase]/mu[phase];
223 template <
class States>
226 int num = states.saturation.size();
227 states.relperm.resize(num);
228 states.relperm_deriv.resize(num);
229 states.mobility.resize(num);
230 states.mobility_deriv.resize(num);
231 #pragma omp parallel for
232 for (
int i = 0; i < num; ++i) {
233 const CompVec& s = states.saturation[i];
234 const PhaseVec& mu = states.viscosity[i];
235 PhaseVec& kr = states.relperm[i];
236 PhaseJacobian& dkr = states.relperm_deriv[i];
237 PhaseVec& lambda = states.mobility[i];
238 PhaseJacobian& dlambda = states.mobility_deriv[i];
241 for (
int phase = 0; phase < numPhases; ++phase) {
242 lambda[phase] = kr[phase]/mu[phase];
243 for (
int p2 = 0; p2 < numPhases; ++p2) {
245 dlambda[phase][p2] = dkr[phase][p2]/mu[phase];
256 CompVec surface_densities_;
264 template <
class Flu
idState>
265 void computeEquilibrium(FluidState& fluid_state)
const
268 const PhaseVec& p = fluid_state.phase_pressure_;
269 const CompVec& z = fluid_state.surface_volume_;
270 PhaseVec& B = fluid_state.formation_volume_factor_;
271 B[Aqua] = pvt_.B(p[Aqua], z, Aqua);
272 B[Vapour] = pvt_.B(p[Vapour], z, Vapour);
273 B[Liquid] = pvt_.B(p[Liquid], z, Liquid);
274 PhaseVec& R = fluid_state.solution_factor_;
276 R[Vapour] = pvt_.R(p[Vapour], z, Vapour);
277 R[Liquid] = pvt_.R(p[Liquid], z, Liquid);
279 dB[Aqua] = pvt_.dBdp(p[Aqua], z, Aqua);
280 dB[Vapour] = pvt_.dBdp(p[Vapour], z, Vapour);
281 dB[Liquid] = pvt_.dBdp(p[Liquid], z, Liquid);
284 dR[Vapour] = pvt_.dRdp(p[Vapour], z, Vapour);
285 dR[Liquid] = pvt_.dRdp(p[Liquid], z, Liquid);
288 PhaseToCompMatrix& At = fluid_state.phase_to_comp_;
289 PhaseVec& u = fluid_state.phase_volume_density_;
290 double& tot_phase_vol_dens = fluid_state.total_phase_volume_density_;
291 PhaseVec& s = fluid_state.saturation_;
292 PhaseVec& cp = fluid_state.phase_compressibility_;
293 double& tot_comp = fluid_state.total_compressibility_;
294 double& exp_term = fluid_state.experimental_term_;
296 computeSingleEquilibrium(B, dB, R, dR, z,
297 At, u, tot_phase_vol_dens,
298 s, cp, tot_comp, exp_term);
301 PhaseVec& mu = fluid_state.viscosity_;
302 mu[Aqua] = pvt_.getViscosity(p[Aqua], z, Aqua);
303 mu[Vapour] = pvt_.getViscosity(p[Vapour], z, Vapour);
304 mu[Liquid] = pvt_.getViscosity(p[Liquid], z, Liquid);
309 static void computeSingleEquilibrium(
const PhaseVec& B,
314 PhaseToCompMatrix& At,
316 double& tot_phase_vol_dens,
328 At[Aqua][Water] = 1.0/B[Aqua];
329 At[Vapour][Gas] = 1.0/B[Vapour];
330 At[Liquid][Gas] = R[Liquid]/B[Liquid];
331 At[Vapour][Oil] = R[Vapour]/B[Vapour];
332 At[Liquid][Oil] = 1.0/B[Liquid];
336 double detR = 1.0 - R[Vapour]*R[Liquid];
337 u[Aqua] = B[Aqua]*z[Water];
338 u[Vapour] = B[Vapour]*(z[Gas] - R[Liquid]*z[Oil])/detR;
339 u[Liquid] = B[Liquid]*(z[Oil] - R[Vapour]*z[Gas])/detR;
340 tot_phase_vol_dens = u[Aqua] + u[Vapour] + u[Liquid];
343 for (
int phase = 0; phase < numPhases; ++phase) {
344 s[phase] = u[phase]/tot_phase_vol_dens;
350 PhaseToCompMatrix dAt(0.0);
351 dAt[Aqua][Water] = -dB[Aqua]/(B[Aqua]*B[Aqua]);
352 dAt[Vapour][Gas] = -dB[Vapour]/(B[Vapour]*B[Vapour]);
353 dAt[Liquid][Oil] = -dB[Liquid]/(B[Liquid]*B[Liquid]);
354 dAt[Liquid][Gas] = dAt[Liquid][Oil]*R[Liquid] + dR[Liquid]/B[Liquid];
355 dAt[Vapour][Oil] = dAt[Vapour][Gas]*R[Vapour] + dR[Vapour]/B[Vapour];
357 PhaseToCompMatrix Ait;
358 Dune::FMatrixHelp::invertMatrix(At, Ait);
360 PhaseToCompMatrix Ct;
361 Dune::FMatrixHelp::multMatrix(dAt, Ait, Ct);
363 cp[Aqua] = Ct[Aqua][Water];
364 cp[Liquid] = Ct[Liquid][Oil] + Ct[Liquid][Gas];
365 cp[Vapour] = Ct[Vapour][Gas] + Ct[Vapour][Oil];
369 PhaseVec tmp1, tmp2, tmp3;
373 exp_term = tmp3[Aqua] + tmp3[Liquid] + tmp3[Gas];
384 std::vector<CompVec> surface_volume_density;
385 std::vector<PhaseVec> phase_pressure;
388 std::vector<PhaseVec> formation_volume_factor;
389 std::vector<PhaseVec> solution_factor;
394 std::vector<PhaseToCompMatrix> state_matrix;
397 std::vector<PhaseVec> mobility;
398 std::vector<PhaseJacobian> mobility_deriv;
401 std::vector<PhaseVec> gravity_potential;
409 std::vector<double> voldiscr;
410 std::vector<double> relvoldiscr;
416 template <
class Gr
id,
class Rock>
417 void computeNew(
const Grid& grid,
420 const typename Grid::Vector gravity,
421 const std::vector<PhaseVec>& cell_pressure,
422 const std::vector<PhaseVec>& face_pressure,
423 const std::vector<CompVec>& cell_z,
424 const CompVec& bdy_z,
427 int num_cells = cell_z.size();
428 assert(num_cells == grid.numCells());
429 const int np = numPhases;
430 const int nc = numComponents;
431 static_assert(np == nc,
"");
432 static_assert(np == 3,
"");
435 cell_data.phase_pressure = cell_pressure;
436 cell_data.surface_volume_density = cell_z;
442 voldiscr.resize(num_cells);
443 relvoldiscr.resize(num_cells);
444 #pragma omp parallel for
445 for (
int cell = 0; cell < num_cells; ++cell) {
446 double pv = rock.
porosity(cell)*grid.cellVolume(cell);
447 voldiscr[cell] = (cell_data.total_phase_volume_density[cell] - 1.0)*pv/dt;
448 relvoldiscr[cell] = std::fabs(cell_data.total_phase_volume_density[cell] - 1.0);
453 computeUpwindProperties(grid, fluid, gravity,
454 cell_pressure, face_pressure,
459 face_data.phase_pressure = face_pressure;
465 template <
class Gr
id>
466 void computeUpwindProperties(
const Grid& grid,
468 const typename Grid::Vector gravity,
469 const std::vector<PhaseVec>& cell_pressure,
470 const std::vector<PhaseVec>& face_pressure,
471 const std::vector<CompVec>& cell_z,
472 const CompVec& bdy_z)
474 int num_faces = face_pressure.size();
475 assert(num_faces == grid.numFaces());
476 bool nonzero_gravity = gravity.two_norm() > 0.0;
477 face_data.state_matrix.resize(num_faces);
478 face_data.mobility.resize(num_faces);
479 face_data.mobility_deriv.resize(num_faces);
480 face_data.gravity_potential.resize(num_faces);
481 face_data.surface_volume_density.resize(num_faces);
482 #pragma omp parallel for
483 for (
int face = 0; face < num_faces; ++face) {
485 typedef typename Grid::Vector Vec;
486 Vec fc = grid.faceCentroid(face);
487 int c[2] = { grid.faceCell(face, 0), grid.faceCell(face, 1) };
492 PhaseVec gravcontrib[2];
493 for (
int j = 0; j < 2; ++j) {
496 phase_p[j] = cell_pressure[c[j]];
498 if (nonzero_gravity) {
500 cdiff -= grid.cellCentroid(c[j]);
501 gravcontrib[j] = fluid.
phaseDensities(&cell_data.state_matrix[c[j]][0][0]);
502 gravcontrib[j] *= (cdiff*gravity);
504 gravcontrib[j] = 0.0;
508 phase_p[j] = face_pressure[face];
510 gravcontrib[j] = 0.0;
520 for (
int phase = 0; phase < numPhases; ++phase) {
521 face_data.gravity_potential[face][phase] = gravcontrib[0][phase] - gravcontrib[1][phase];
522 pot[0][phase] = phase_p[0][phase] + face_data.gravity_potential[face][phase];
523 pot[1][phase] = phase_p[1][phase];
532 double face_z_factor = 0.5;
533 PhaseVec phase_mob[2];
534 PhaseJacobian phasemob_deriv[2];
535 for (
int j = 0; j < 2; ++j) {
537 face_z += cell_z[c[j]];
538 phase_mob[j] = cell_data.mobility[c[j]];
539 phasemob_deriv[j] = cell_data.mobility_deriv[c[j]];
540 }
else if (pot[j][Liquid] > pot[(j+1)%2][Liquid]) {
544 phase_mob[j] = bdy_state.mobility_;
545 phasemob_deriv[j] = bdy_state.dmobility_;
553 face_z *= face_z_factor;
554 face_data.surface_volume_density[face] = face_z;
557 for (
int phase = 0; phase < numPhases; ++phase) {
558 if (pot[0][phase] == pot[1][phase]) {
560 double aver = 0.5*(phase_mob[0][phase] + phase_mob[1][phase]);
561 face_data.mobility[face][phase] = aver;
562 for (
int p2 = 0; p2 < numPhases; ++p2) {
563 face_data.mobility_deriv[face][phase][p2] = phasemob_deriv[0][phase][p2]
564 + phasemob_deriv[1][phase][p2];
568 int upwind = pot[0][phase] > pot[1][phase] ? 0 : 1;
569 face_data.mobility[face][phase] = phase_mob[upwind][phase];
570 for (
int p2 = 0; p2 < numPhases; ++p2) {
571 face_data.mobility_deriv[face][phase][p2] = phasemob_deriv[upwind][phase][p2];
584 #endif // OPM_BLACKOILFLUID_HEADER_INCLUDED
PhaseVec phaseDensities(const double *A) const
Definition: BlackoilFluid.hpp:85
Definition: BlackoilDefs.hpp:33
Definition: BlackoilFluid.hpp:381
void computeMobilitiesNoDerivs(States &states) const
Input: s, mu Output: kr, lambda.
Definition: BlackoilFluid.hpp:202
void computePvtNoDerivs(States &states) const
Input: p, z Output: B, R, mu.
Definition: BlackoilFluid.hpp:115
static void dkr(krContainerT &dkr, const Params ¶ms, const SatContainerT &saturations, Scalar)
The saturation derivatives of relative permeability of all phases.
Definition: FluidMatrixInteractionBlackoil.hpp:158
Class responsible for computing all fluid properties from face pressures and composition.
Definition: BlackoilFluid.hpp:41
double porosity(int cell_index) const
Read-access to porosity.
Definition: Rock_impl.hpp:76
Definition: BlackoilPVT.hpp:33
A property class for porous media rock.
Definition: ImplicitTransportDefs.hpp:82
Fluid state for a black oil model.
Definition: FluidStateBlackoil.hpp:32
void computeBAndR(States &states) const
Input: p, z Output: B, R.
Definition: BlackoilFluid.hpp:102
Definition: BlackoilFluid.hpp:405
void computeStateMatrix(States &states) const
Input: B, R Output: A.
Definition: BlackoilFluid.hpp:144
Multiple fluid states for a black oil model.
Definition: FluidStateBlackoil.hpp:57
static void kr(krContainerT &kr, const Params ¶ms, const SatContainerT &saturations, Scalar)
The relative permeability of all phases.
Definition: FluidMatrixInteractionBlackoil.hpp:132
void computePvt(States &states) const
Input: p, z Output: B, dB/dp, R, dR/dp, mu.
Definition: BlackoilFluid.hpp:127
void computeMobilities(States &states) const
Input: s, mu Output: kr, dkr/ds, lambda, dlambda/ds.
Definition: BlackoilFluid.hpp:224
void computePvtDepending(States &states) const
Input: z, B, dB/dp, R, dR/dp Output: A, u, sum(u), s, c, cT, ex.
Definition: BlackoilFluid.hpp:169