diff --git a/dune/porsol/blackoil/BlackoilFluid.hpp b/dune/porsol/blackoil/BlackoilFluid.hpp
new file mode 100644
index 00000000..d1cd6df2
--- /dev/null
+++ b/dune/porsol/blackoil/BlackoilFluid.hpp
@@ -0,0 +1,585 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_BLACKOILFLUID_HEADER_INCLUDED
+#define OPM_BLACKOILFLUID_HEADER_INCLUDED
+
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+
+namespace Opm
+{
+
+
+ /// Forward declaration for typedef i BlackoilFluid.
+ class BlackoilFluidData;
+
+
+ /// Class responsible for computing all fluid properties from
+ /// face pressures and composition.
+ class BlackoilFluid : public BlackoilDefs
+ {
+ public:
+ typedef FluidStateBlackoil FluidState;
+ typedef BlackoilFluidData FluidData;
+
+ void init(const Dune::EclipseGridParser& parser)
+ {
+ fmi_params_.init(parser);
+ // FluidSystemBlackoil<>::init(parser);
+ pvt_.init(parser);
+ const std::vector& dens = parser.getDENSITY().densities_[0];
+ surface_densities_[Oil] = dens[0];
+ surface_densities_[Water] = dens[1];
+ surface_densities_[Gas] = dens[2];
+ }
+
+ FluidState computeState(PhaseVec phase_pressure, CompVec z) const
+ {
+ FluidState state;
+ state.temperature_ = 300;
+ state.phase_pressure_ = phase_pressure;
+ state.surface_volume_ = z;
+ // FluidSystemBlackoil<>::computeEquilibrium(state); // Sets everything but relperm and mobility.
+ computeEquilibrium(state);
+ FluidMatrixInteractionBlackoil::kr(state.relperm_, fmi_params_, state.saturation_, state.temperature_);
+ FluidMatrixInteractionBlackoil::dkr(state.drelperm_, fmi_params_, state.saturation_, state.temperature_);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ state.mobility_[phase] = state.relperm_[phase]/state.viscosity_[phase];
+ for (int p2 = 0; p2 < numPhases; ++p2) {
+ // Ignoring pressure variation in viscosity for this one.
+ state.dmobility_[phase][p2] = state.drelperm_[phase][p2]/state.viscosity_[phase];
+ }
+ }
+ return state;
+ }
+
+
+ const CompVec& surfaceDensities() const
+ {
+ return surface_densities_;
+ }
+
+ /// \param[in] A state matrix in fortran ordering
+ PhaseVec phaseDensities(const double* A) const
+ {
+ PhaseVec phase_dens(0.0);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ for (int comp = 0; comp < numComponents; ++comp) {
+ phase_dens[phase] += A[numComponents*phase + comp]*surface_densities_[comp];
+ }
+ }
+ return phase_dens;
+ }
+
+
+ // ---- New interface ----
+
+ /// Input: p, z
+ /// Output: B, R
+ template
+ void computeBAndR(States& states) const
+ {
+ const std::vector& p = states.phase_pressure;
+ const std::vector& z = states.surface_volume_density;
+ std::vector& B = states.formation_volume_factor;
+ std::vector& R = states.solution_factor;
+ pvt_.B(p, z, B);
+ pvt_.R(p, z, R);
+ }
+
+ /// Input: p, z
+ /// Output: B, R, mu
+ template
+ void computePvtNoDerivs(States& states) const
+ {
+ computeBAndR(states);
+ const std::vector& p = states.phase_pressure;
+ const std::vector& z = states.surface_volume_density;
+ std::vector& mu = states.viscosity;
+ pvt_.getViscosity(p, z, mu);
+ }
+
+ /// Input: p, z
+ /// Output: B, dB/dp, R, dR/dp, mu
+ template
+ void computePvt(States& states) const
+ {
+ const std::vector& p = states.phase_pressure;
+ const std::vector& z = states.surface_volume_density;
+ std::vector& B = states.formation_volume_factor;
+ std::vector& dB = states.formation_volume_factor_deriv;
+ std::vector& R = states.solution_factor;
+ std::vector& dR = states.solution_factor_deriv;
+ std::vector& mu = states.viscosity;
+ pvt_.dBdp(p, z, B, dB);
+ pvt_.dRdp(p, z, R, dR);
+ pvt_.getViscosity(p, z, mu);
+ }
+
+ /// Input: B, R
+ /// Output: A
+ template
+ void computeStateMatrix(States& states) const
+ {
+ int num = states.formation_volume_factor.size();
+ states.state_matrix.resize(num);
+#pragma omp parallel for
+ for (int i = 0; i < num; ++i) {
+ const PhaseVec& B = states.formation_volume_factor[i];
+ const PhaseVec& R = states.solution_factor[i];
+ PhaseToCompMatrix& At = states.state_matrix[i];
+ // Set the A matrix (A = RB^{-1})
+ // Using A transposed (At) since we really want Fortran ordering:
+ // ultimately that is what the opmpressure C library expects.
+ At = 0.0;
+ At[Aqua][Water] = 1.0/B[Aqua];
+ At[Vapour][Gas] = 1.0/B[Vapour];
+ At[Liquid][Gas] = R[Liquid]/B[Liquid];
+ At[Vapour][Oil] = R[Vapour]/B[Vapour];
+ At[Liquid][Oil] = 1.0/B[Liquid];
+ }
+ }
+
+
+ /// Input: z, B, dB/dp, R, dR/dp
+ /// Output: A, u, sum(u), s, c, cT, ex
+ template
+ void computePvtDepending(States& states) const
+ {
+ int num = states.formation_volume_factor.size();
+ states.state_matrix.resize(num);
+ states.phase_volume_density.resize(num);
+ states.total_phase_volume_density.resize(num);
+ states.saturation.resize(num);
+ states.phase_compressibility.resize(num);
+ states.total_compressibility.resize(num);
+ states.experimental_term.resize(num);
+#pragma omp parallel for
+ for (int i = 0; i < num; ++i) {
+ const CompVec& z = states.surface_volume_density[i];
+ const PhaseVec& B = states.formation_volume_factor[i];
+ const PhaseVec& dB = states.formation_volume_factor_deriv[i];
+ const PhaseVec& R = states.solution_factor[i];
+ const PhaseVec& dR = states.solution_factor_deriv[i];
+ PhaseToCompMatrix& At = states.state_matrix[i];
+ PhaseVec& u = states.phase_volume_density[i];
+ double& tot_phase_vol_dens = states.total_phase_volume_density[i];
+ PhaseVec& s = states.saturation[i];
+ PhaseVec& cp = states.phase_compressibility[i];
+ double& tot_comp = states.total_compressibility[i];
+ double& exp_term = states.experimental_term[i];
+ computeSingleEquilibrium(B, dB, R, dR, z,
+ At, u, tot_phase_vol_dens,
+ s, cp, tot_comp, exp_term);
+ }
+ }
+
+ /// Input: s, mu
+ /// Output: kr, lambda
+ template
+ void computeMobilitiesNoDerivs(States& states) const
+ {
+ int num = states.saturation.size();
+ states.relperm.resize(num);
+ states.mobility.resize(num);
+#pragma omp parallel for
+ for (int i = 0; i < num; ++i) {
+ const CompVec& s = states.saturation[i];
+ const PhaseVec& mu = states.viscosity[i];
+ PhaseVec& kr = states.relperm[i];
+ PhaseVec& lambda = states.mobility[i];
+ FluidMatrixInteractionBlackoil::kr(kr, fmi_params_, s, 300.0);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ lambda[phase] = kr[phase]/mu[phase];
+ }
+
+ }
+ }
+
+ /// Input: s, mu
+ /// Output: kr, dkr/ds, lambda, dlambda/ds
+ template
+ void computeMobilities(States& states) const
+ {
+ int num = states.saturation.size();
+ states.relperm.resize(num);
+ states.relperm_deriv.resize(num);
+ states.mobility.resize(num);
+ states.mobility_deriv.resize(num);
+#pragma omp parallel for
+ for (int i = 0; i < num; ++i) {
+ const CompVec& s = states.saturation[i];
+ const PhaseVec& mu = states.viscosity[i];
+ PhaseVec& kr = states.relperm[i];
+ PhaseJacobian& dkr = states.relperm_deriv[i];
+ PhaseVec& lambda = states.mobility[i];
+ PhaseJacobian& dlambda = states.mobility_deriv[i];
+ FluidMatrixInteractionBlackoil::kr(kr, fmi_params_, s, 300.0);
+ FluidMatrixInteractionBlackoil::dkr(dkr, fmi_params_, s, 300.0);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ lambda[phase] = kr[phase]/mu[phase];
+ for (int p2 = 0; p2 < numPhases; ++p2) {
+ // Ignoring pressure variation in viscosity for this one.
+ dlambda[phase][p2] = dkr[phase][p2]/mu[phase];
+ }
+ }
+
+ }
+ }
+
+
+ private:
+ BlackoilPVT pvt_;
+ FluidMatrixInteractionBlackoilParams fmi_params_;
+ CompVec surface_densities_;
+
+
+ /*!
+ * \brief Assuming the surface volumes and the pressures of all
+ * phases are known, compute everything except relperm and
+ * mobility.
+ */
+ template
+ void computeEquilibrium(FluidState& fluid_state) const
+ {
+ // Get B and R factors.
+ const PhaseVec& p = fluid_state.phase_pressure_;
+ const CompVec& z = fluid_state.surface_volume_;
+ PhaseVec& B = fluid_state.formation_volume_factor_;
+ B[Aqua] = pvt_.B(p[Aqua], z, Aqua);
+ B[Vapour] = pvt_.B(p[Vapour], z, Vapour);
+ B[Liquid] = pvt_.B(p[Liquid], z, Liquid);
+ PhaseVec& R = fluid_state.solution_factor_;
+ R[Aqua] = 0.0;
+ R[Vapour] = pvt_.R(p[Vapour], z, Vapour);
+ R[Liquid] = pvt_.R(p[Liquid], z, Liquid);
+ PhaseVec dB;
+ dB[Aqua] = pvt_.dBdp(p[Aqua], z, Aqua);
+ dB[Vapour] = pvt_.dBdp(p[Vapour], z, Vapour);
+ dB[Liquid] = pvt_.dBdp(p[Liquid], z, Liquid);
+ PhaseVec dR;
+ dR[Aqua] = 0.0;
+ dR[Vapour] = pvt_.dRdp(p[Vapour], z, Vapour);
+ dR[Liquid] = pvt_.dRdp(p[Liquid], z, Liquid);
+
+ // Convenience vars.
+ PhaseToCompMatrix& At = fluid_state.phase_to_comp_;
+ PhaseVec& u = fluid_state.phase_volume_density_;
+ double& tot_phase_vol_dens = fluid_state.total_phase_volume_density_;
+ PhaseVec& s = fluid_state.saturation_;
+ PhaseVec& cp = fluid_state.phase_compressibility_;
+ double& tot_comp = fluid_state.total_compressibility_;
+ double& exp_term = fluid_state.experimental_term_;
+
+ computeSingleEquilibrium(B, dB, R, dR, z,
+ At, u, tot_phase_vol_dens,
+ s, cp, tot_comp, exp_term);
+
+ // Compute viscosities.
+ PhaseVec& mu = fluid_state.viscosity_;
+ mu[Aqua] = pvt_.getViscosity(p[Aqua], z, Aqua);
+ mu[Vapour] = pvt_.getViscosity(p[Vapour], z, Vapour);
+ mu[Liquid] = pvt_.getViscosity(p[Liquid], z, Liquid);
+ }
+
+
+
+ static void computeSingleEquilibrium(const PhaseVec& B,
+ const PhaseVec& dB,
+ const PhaseVec& R,
+ const PhaseVec& dR,
+ const CompVec& z,
+ PhaseToCompMatrix& At,
+ PhaseVec& u,
+ double& tot_phase_vol_dens,
+ PhaseVec& s,
+ PhaseVec& cp,
+ double& tot_comp,
+ double& exp_term)
+ {
+ // Set the A matrix (A = RB^{-1})
+ // Using At since we really want Fortran ordering
+ // (since ultimately that is what the opmpressure
+ // C library expects).
+ // PhaseToCompMatrix& At = fluid_state.phase_to_comp_[i];
+ At = 0.0;
+ At[Aqua][Water] = 1.0/B[Aqua];
+ At[Vapour][Gas] = 1.0/B[Vapour];
+ At[Liquid][Gas] = R[Liquid]/B[Liquid];
+ At[Vapour][Oil] = R[Vapour]/B[Vapour];
+ At[Liquid][Oil] = 1.0/B[Liquid];
+
+ // Update phase volumes. This is the same as multiplying with A^{-1}
+ // PhaseVec& u = fluid_state.phase_volume_density_[i];
+ double detR = 1.0 - R[Vapour]*R[Liquid];
+ u[Aqua] = B[Aqua]*z[Water];
+ u[Vapour] = B[Vapour]*(z[Gas] - R[Liquid]*z[Oil])/detR;
+ u[Liquid] = B[Liquid]*(z[Oil] - R[Vapour]*z[Gas])/detR;
+ tot_phase_vol_dens = u[Aqua] + u[Vapour] + u[Liquid];
+
+ // PhaseVec& s = fluid_state.saturation_[i];
+ for (int phase = 0; phase < numPhases; ++phase) {
+ s[phase] = u[phase]/tot_phase_vol_dens;
+ }
+
+ // Phase compressibilities.
+ // PhaseVec& cp = fluid_state.phase_compressibility_[i];
+ // Set the derivative of the A matrix (A = RB^{-1})
+ PhaseToCompMatrix dAt(0.0);
+ dAt[Aqua][Water] = -dB[Aqua]/(B[Aqua]*B[Aqua]);
+ dAt[Vapour][Gas] = -dB[Vapour]/(B[Vapour]*B[Vapour]);
+ dAt[Liquid][Oil] = -dB[Liquid]/(B[Liquid]*B[Liquid]); // Different order than above.
+ dAt[Liquid][Gas] = dAt[Liquid][Oil]*R[Liquid] + dR[Liquid]/B[Liquid];
+ dAt[Vapour][Oil] = dAt[Vapour][Gas]*R[Vapour] + dR[Vapour]/B[Vapour];
+
+ PhaseToCompMatrix Ait;
+ Dune::FMatrixHelp::invertMatrix(At, Ait);
+
+ PhaseToCompMatrix Ct;
+ Dune::FMatrixHelp::multMatrix(dAt, Ait, Ct);
+
+ cp[Aqua] = Ct[Aqua][Water];
+ cp[Liquid] = Ct[Liquid][Oil] + Ct[Liquid][Gas];
+ cp[Vapour] = Ct[Vapour][Gas] + Ct[Vapour][Oil];
+ tot_comp = cp*s;
+
+ // Experimental term.
+ PhaseVec tmp1, tmp2, tmp3;
+ Ait.mtv(z, tmp1);
+ dAt.mtv(tmp1, tmp2);
+ Ait.mtv(tmp2, tmp3);
+ exp_term = tmp3[Aqua] + tmp3[Liquid] + tmp3[Gas];
+ }
+
+
+ };
+
+
+
+ struct FaceFluidData : public BlackoilDefs
+ {
+ // Canonical state variables.
+ std::vector surface_volume_density; // z
+ std::vector phase_pressure; // p
+
+ // Variables from PVT functions.
+ std::vector formation_volume_factor; // B
+ std::vector solution_factor; // R
+
+ // Variables computed from PVT data.
+ // The A matrices are all in Fortran order (or, equivalently,
+ // we store the transposes).
+ std::vector state_matrix; // A' = (RB^{-1})'
+
+ // Variables computed from saturation.
+ std::vector mobility; // lambda
+ std::vector mobility_deriv; // dlambda/ds
+
+ // Gravity and/or capillary pressure potential differences.
+ std::vector gravity_potential; // (\rho g \delta z)-ish contribution per face
+ };
+
+
+ struct AllFluidData : public BlackoilDefs
+ {
+ // Per-cell data
+ AllFluidStates cell_data;
+ std::vector voldiscr;
+ std::vector relvoldiscr;
+
+ // Per-face data.
+ FaceFluidData face_data;
+
+ public:
+ template
+ void computeNew(const Grid& grid,
+ const Rock& rock,
+ const BlackoilFluid& fluid,
+ const typename Grid::Vector gravity,
+ const std::vector& cell_pressure,
+ const std::vector& face_pressure,
+ const std::vector& cell_z,
+ const CompVec& bdy_z,
+ const double dt)
+ {
+ int num_cells = cell_z.size();
+ ASSERT(num_cells == grid.numCells());
+ const int np = numPhases;
+ const int nc = numComponents;
+ BOOST_STATIC_ASSERT(np == nc);
+ BOOST_STATIC_ASSERT(np == 3);
+
+ // p, z -> B, dB, R, dR, mu, A, dA, u, sum(u), s, c, cT, ex, kr, dkr, lambda, dlambda.
+ cell_data.phase_pressure = cell_pressure;
+ cell_data.surface_volume_density = cell_z;
+ fluid.computePvt(cell_data);
+ fluid.computePvtDepending(cell_data);
+ fluid.computeMobilities(cell_data);
+
+ // Compute volume discrepancies.
+ voldiscr.resize(num_cells);
+ relvoldiscr.resize(num_cells);
+#pragma omp parallel for
+ for (int cell = 0; cell < num_cells; ++cell) {
+ double pv = rock.porosity(cell)*grid.cellVolume(cell);
+ voldiscr[cell] = (cell_data.total_phase_volume_density[cell] - 1.0)*pv/dt;
+ relvoldiscr[cell] = std::fabs(cell_data.total_phase_volume_density[cell] - 1.0);
+ }
+
+
+ // Compute upwinded face properties, including z.
+ computeUpwindProperties(grid, fluid, gravity,
+ cell_pressure, face_pressure,
+ cell_z, bdy_z);
+
+ // Compute state matrices for faces.
+ // p, z -> B, R, A
+ face_data.phase_pressure = face_pressure;
+ fluid.computeBAndR(face_data);
+ fluid.computeStateMatrix(face_data);
+ }
+
+
+ template
+ void computeUpwindProperties(const Grid& grid,
+ const BlackoilFluid& fluid,
+ const typename Grid::Vector gravity,
+ const std::vector& cell_pressure,
+ const std::vector& face_pressure,
+ const std::vector& cell_z,
+ const CompVec& bdy_z)
+ {
+ int num_faces = face_pressure.size();
+ ASSERT(num_faces == grid.numFaces());
+ bool nonzero_gravity = gravity.two_norm() > 0.0;
+ face_data.state_matrix.resize(num_faces);
+ face_data.mobility.resize(num_faces);
+ face_data.mobility_deriv.resize(num_faces);
+ face_data.gravity_potential.resize(num_faces);
+ face_data.surface_volume_density.resize(num_faces);
+#pragma omp parallel for
+ for (int face = 0; face < num_faces; ++face) {
+ // Obtain properties from both sides of the face.
+ typedef typename Grid::Vector Vec;
+ Vec fc = grid.faceCentroid(face);
+ int c[2] = { grid.faceCell(face, 0), grid.faceCell(face, 1) };
+
+ // Get pressures and compute gravity contributions,
+ // to decide upwind directions.
+ PhaseVec phase_p[2];
+ PhaseVec gravcontrib[2];
+ for (int j = 0; j < 2; ++j) {
+ if (c[j] >= 0) {
+ // Pressures
+ phase_p[j] = cell_pressure[c[j]];
+ // Gravity contribution.
+ if (nonzero_gravity) {
+ Vec cdiff = fc;
+ cdiff -= grid.cellCentroid(c[j]);
+ gravcontrib[j] = fluid.phaseDensities(&cell_data.state_matrix[c[j]][0][0]);
+ gravcontrib[j] *= (cdiff*gravity);
+ } else {
+ gravcontrib[j] = 0.0;
+ }
+ } else {
+ // Pressures
+ phase_p[j] = face_pressure[face];
+ // Gravity contribution.
+ gravcontrib[j] = 0.0;
+ }
+ }
+
+ // Gravity contribution:
+ // gravcapf = rho_1*g*(z_12 - z_1) - rho_2*g*(z_12 - z_2)
+ // where _1 and _2 refers to two neigbour cells, z is the
+ // z coordinate of the centroid, and z_12 is the face centroid.
+ // Also compute the potentials.
+ PhaseVec pot[2];
+ for (int phase = 0; phase < numPhases; ++phase) {
+ face_data.gravity_potential[face][phase] = gravcontrib[0][phase] - gravcontrib[1][phase];
+ pot[0][phase] = phase_p[0][phase] + face_data.gravity_potential[face][phase];
+ pot[1][phase] = phase_p[1][phase];
+ }
+
+ // Now we can easily find the upwind direction for every phase,
+ // we can also tell which boundary faces are inflow bdys.
+
+ // Compute face_z, which is averaged from the cells, unless on outflow or noflow bdy.
+ // Get mobilities and derivatives.
+ CompVec face_z(0.0);
+ double face_z_factor = 0.5;
+ PhaseVec phase_mob[2];
+ PhaseJacobian phasemob_deriv[2];
+ for (int j = 0; j < 2; ++j) {
+ if (c[j] >= 0) {
+ face_z += cell_z[c[j]];
+ phase_mob[j] = cell_data.mobility[c[j]];
+ phasemob_deriv[j] = cell_data.mobility_deriv[c[j]];
+ } else if (pot[j][Liquid] > pot[(j+1)%2][Liquid]) {
+ // Inflow boundary.
+ face_z += bdy_z;
+ FluidStateBlackoil bdy_state = fluid.computeState(face_pressure[face], bdy_z);
+ phase_mob[j] = bdy_state.mobility_;
+ phasemob_deriv[j] = bdy_state.dmobility_;
+ } else {
+ // For outflow or noflow boundaries, only cell z is used.
+ face_z_factor = 1.0;
+ // Also, make sure the boundary data are not used for mobilities.
+ pot[j] = -1e100;
+ }
+ }
+ face_z *= face_z_factor;
+ face_data.surface_volume_density[face] = face_z;
+
+ // Computing upwind mobilities and derivatives
+ for (int phase = 0; phase < numPhases; ++phase) {
+ if (pot[0][phase] == pot[1][phase]) {
+ // Average.
+ double aver = 0.5*(phase_mob[0][phase] + phase_mob[1][phase]);
+ face_data.mobility[face][phase] = aver;
+ for (int p2 = 0; p2 < numPhases; ++p2) {
+ face_data.mobility_deriv[face][phase][p2] = phasemob_deriv[0][phase][p2]
+ + phasemob_deriv[1][phase][p2];
+ }
+ } else {
+ // Upwind.
+ int upwind = pot[0][phase] > pot[1][phase] ? 0 : 1;
+ face_data.mobility[face][phase] = phase_mob[upwind][phase];
+ for (int p2 = 0; p2 < numPhases; ++p2) {
+ face_data.mobility_deriv[face][phase][p2] = phasemob_deriv[upwind][phase][p2];
+ }
+ }
+ }
+ }
+ }
+
+
+ };
+
+
+} // namespace Opm
+
+#endif // OPM_BLACKOILFLUID_HEADER_INCLUDED
diff --git a/dune/porsol/blackoil/fluid/BlackoilComponent.hpp b/dune/porsol/blackoil/fluid/BlackoilComponent.hpp
new file mode 100644
index 00000000..61294a68
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/BlackoilComponent.hpp
@@ -0,0 +1,165 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
+#define OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
+
+
+#include
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Components
+
+ * \brief
+ * A component class for the black oil model, intended to be used
+ * for all three components.
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template
+class BlackoilComponent
+{
+public:
+ /*!
+ * \brief A human readable name for the component.
+ */
+ static const char *name()
+ {
+ return "BlackoilComponent";
+ }
+
+ /*!
+ * \brief The molar mass in [kg] of the component.
+ */
+ static Scalar molarMass()
+ { DUNE_THROW(Dune::NotImplemented, "Component::molarMass()"); }
+
+ /*!
+ * \brief Returns the critical temperature in [K] of the component.
+ */
+ static Scalar criticalTemperature()
+ { DUNE_THROW(Dune::NotImplemented, "Component::criticalTemperature()"); }
+
+ /*!
+ * \brief Returns the critical pressure in [Pa] of the component.
+ */
+ static Scalar criticalPressure()
+ { DUNE_THROW(Dune::NotImplemented, "Component::criticalPressure()"); }
+
+ /*!
+ * \brief Returns the temperature in [K] at the component's triple point.
+ */
+ static Scalar tripleTemperature()
+ { DUNE_THROW(Dune::NotImplemented, "Component::tripleTemperature()"); }
+
+ /*!
+ * \brief Returns the pressure in [Pa] at the component's triple point.
+ */
+ static Scalar triplePressure()
+ { DUNE_THROW(Dune::NotImplemented, "Component::triplePressure()"); }
+
+ /*!
+ * \brief The vapor pressure in [Pa] of the component at a given
+ * temperature in [K].
+ *
+ * \param T temperature of the component in [K]
+ */
+ static Scalar vaporPressure(Scalar T)
+ { DUNE_THROW(Dune::NotImplemented, "Component::vaporPressure()"); }
+
+ /*!
+ * \brief The density in [kg/m^3] of the component at a given pressure in [Pa] and temperature in [K].
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static Scalar gasDensity(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::density()"); }
+
+ /*!
+ * \brief The density [kg/m^3] of the liquid component at a given pressure in [Pa] and temperature in [K].
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::density()"); }
+
+ /*!
+ * \brief Specific enthalpy [J/kg] of the pure component in gas.
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::gasEnthalpy()"); }
+
+ /*!
+ * \brief Specific enthalpy [J/kg] of the pure component in liquid.
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::liquidEnthalpy()"); }
+
+ /*!
+ * \brief Specific internal energy [J/kg] of the pure component in gas.
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::gasInternalEnergy()"); }
+
+ /*!
+ * \brief Specific internal energy [J/kg] of pure the pure component in liquid.
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::liquidInternalEnergy()"); }
+
+ /*!
+ * \brief The dynamic viscosity [Pa*s] of the pure component at a given pressure in [Pa] and temperature in [K].
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static Scalar gasViscosity(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::gasViscosity()"); }
+
+ /*!
+ * \brief The dynamic liquid viscosity [Pa*s] of the pure component.
+ *
+ * \param temperature temperature of component in [K]
+ * \param pressure pressure of component in [Pa]
+ */
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ { DUNE_THROW(Dune::NotImplemented, "Component::liquidViscosity()"); }
+
+};
+
+} // end namepace
+
+#endif // OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
diff --git a/dune/porsol/blackoil/fluid/BlackoilDefs.hpp b/dune/porsol/blackoil/fluid/BlackoilDefs.hpp
new file mode 100644
index 00000000..e6e4afcb
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/BlackoilDefs.hpp
@@ -0,0 +1,50 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_BLACKOILDEFS_HEADER_INCLUDED
+#define OPM_BLACKOILDEFS_HEADER_INCLUDED
+
+
+#include
+#include
+#include
+
+namespace Opm
+{
+
+ class BlackoilDefs
+ {
+ public:
+ enum { numComponents = 3 };
+ enum { numPhases = 3 };
+
+ enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 };
+ enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
+
+ typedef double Scalar;
+ typedef Dune::FieldVector CompVec;
+ typedef Dune::FieldVector PhaseVec;
+ BOOST_STATIC_ASSERT(int(numComponents) == int(numPhases));
+ typedef Dune::FieldMatrix PhaseToCompMatrix;
+ typedef Dune::FieldMatrix PhaseJacobian;
+ };
+
+} // namespace Opm
+
+#endif // OPM_BLACKOILDEFS_HEADER_INCLUDED
diff --git a/dune/porsol/blackoil/fluid/BlackoilPVT.cpp b/dune/porsol/blackoil/fluid/BlackoilPVT.cpp
new file mode 100644
index 00000000..c35029f0
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/BlackoilPVT.cpp
@@ -0,0 +1,202 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+
+#include "BlackoilPVT.hpp"
+#include
+#include
+#include "MiscibilityDead.hpp"
+#include "MiscibilityLiveOil.hpp"
+#include "MiscibilityLiveGas.hpp"
+#include "MiscibilityWater.hpp"
+#include
+#include
+
+using namespace Dune;
+
+namespace Opm
+{
+
+
+ void BlackoilPVT::init(const Dune::EclipseGridParser& parser)
+ {
+ typedef std::vector > > table_t;
+ region_number_ = 0;
+
+ // Surface densities. Accounting for different orders in eclipse and our code.
+ if (parser.hasField("DENSITY")) {
+ const int region_number = 0;
+ enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
+ const std::vector& d = parser.getDENSITY().densities_[region_number];
+ densities_[Aqua] = d[ECL_water];
+ densities_[Vapour] = d[ECL_gas];
+ densities_[Liquid] = d[ECL_oil];
+ } else {
+ THROW("Input is missing DENSITY\n");
+ }
+
+ // Water PVT
+ if (parser.hasField("PVTW")) {
+ water_props_.reset(new MiscibilityWater(parser.getPVTW().pvtw_));
+ } else {
+ water_props_.reset(new MiscibilityWater(0.5*Dune::prefix::centi*Dune::unit::Poise)); // Eclipse 100 default
+ }
+
+ // Oil PVT
+ if (parser.hasField("PVDO")) {
+ oil_props_.reset(new MiscibilityDead(parser.getPVDO().pvdo_));
+ } else if (parser.hasField("PVTO")) {
+ oil_props_.reset(new MiscibilityLiveOil(parser.getPVTO().pvto_));
+ } else if (parser.hasField("PVCDO")) {
+ oil_props_.reset(new MiscibilityWater(parser.getPVCDO().pvcdo_));
+ } else {
+ THROW("Input is missing PVDO and PVTO\n");
+ }
+
+ // Gas PVT
+ if (parser.hasField("PVDG")) {
+ gas_props_.reset(new MiscibilityDead(parser.getPVDG().pvdg_));
+ } else if (parser.hasField("PVTG")) {
+ gas_props_.reset(new MiscibilityLiveGas(parser.getPVTG().pvtg_));
+ } else {
+ THROW("Input is missing PVDG and PVTG\n");
+ }
+ }
+
+ BlackoilPVT::CompVec BlackoilPVT::surfaceDensities() const
+ {
+ return densities_;
+ }
+
+ double BlackoilPVT::getViscosity(double press, const CompVec& surfvol, PhaseIndex phase) const
+ {
+ return propsForPhase(phase).getViscosity(region_number_, press, surfvol);
+ }
+
+ double BlackoilPVT::B(double press, const CompVec& surfvol, PhaseIndex phase) const
+ {
+ return propsForPhase(phase).B(region_number_, press, surfvol);
+ }
+
+ double BlackoilPVT::dBdp(double press, const CompVec& surfvol, PhaseIndex phase) const
+ {
+ return propsForPhase(phase).dBdp(region_number_, press, surfvol);
+ }
+
+ double BlackoilPVT::R(double press, const CompVec& surfvol, PhaseIndex phase) const
+ {
+ return propsForPhase(phase).R(region_number_, press, surfvol);
+ }
+
+ double BlackoilPVT::dRdp(double press, const CompVec& surfvol, PhaseIndex phase) const
+ {
+ return propsForPhase(phase).dRdp(region_number_, press, surfvol);
+ }
+
+ const MiscibilityProps& BlackoilPVT::propsForPhase(PhaseIndex phase) const
+ {
+ switch (phase) {
+ case Aqua:
+ return *water_props_;
+ case Liquid:
+ return *oil_props_;
+ case Vapour:
+ return *gas_props_;
+ default:
+ THROW("Unknown phase accessed: " << phase);
+ }
+ }
+
+ void BlackoilPVT::getViscosity(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output) const
+ {
+ int num = pressures.size();
+ output.resize(num);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ propsForPhase(PhaseIndex(phase)).getViscosity(pressures, surfvol, phase, data1_);
+ for (int i = 0; i < num; ++i) {
+ output[i][phase] = data1_[i];
+ }
+ }
+ }
+
+ void BlackoilPVT::B(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output) const
+ {
+ int num = pressures.size();
+ output.resize(num);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ propsForPhase(PhaseIndex(phase)).B(pressures, surfvol, phase, data1_);
+ for (int i = 0; i < num; ++i) {
+ output[i][phase] = data1_[i];
+ }
+ }
+ }
+
+ void BlackoilPVT::dBdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output_B,
+ std::vector& output_dBdp) const
+ {
+ int num = pressures.size();
+ output_B.resize(num);
+ output_dBdp.resize(num);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ propsForPhase(PhaseIndex(phase)).dBdp(pressures, surfvol, phase, data1_, data2_);
+ for (int i = 0; i < num; ++i) {
+ output_B[i][phase] = data1_[i];
+ output_dBdp[i][phase] = data2_[i];
+ }
+ }
+ }
+
+ void BlackoilPVT::R(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output) const
+ {
+ int num = pressures.size();
+ output.resize(num);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ propsForPhase(PhaseIndex(phase)).R(pressures, surfvol, phase, data1_);
+ for (int i = 0; i < num; ++i) {
+ output[i][phase] = data1_[i];
+ }
+ }
+ }
+
+ void BlackoilPVT::dRdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output_R,
+ std::vector& output_dRdp) const
+ {
+ int num = pressures.size();
+ output_R.resize(num);
+ output_dRdp.resize(num);
+ for (int phase = 0; phase < numPhases; ++phase) {
+ propsForPhase(PhaseIndex(phase)).dRdp(pressures, surfvol, phase, data1_, data2_);
+ for (int i = 0; i < num; ++i) {
+ output_R[i][phase] = data1_[i];
+ output_dRdp[i][phase] = data2_[i];
+ }
+ }
+ }
+
+} // namespace Opm
diff --git a/dune/porsol/blackoil/fluid/BlackoilPVT.hpp b/dune/porsol/blackoil/fluid/BlackoilPVT.hpp
new file mode 100644
index 00000000..d5469223
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/BlackoilPVT.hpp
@@ -0,0 +1,88 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_BLACKOILPVT_HEADER_INCLUDED
+#define OPM_BLACKOILPVT_HEADER_INCLUDED
+
+
+#include "MiscibilityProps.hpp"
+#include "BlackoilDefs.hpp"
+#include
+#include
+#include
+
+
+namespace Opm
+{
+ class BlackoilPVT : public BlackoilDefs
+ {
+ public:
+ void init(const Dune::EclipseGridParser& ep);
+
+ double getViscosity(double press,
+ const CompVec& surfvol,
+ PhaseIndex phase) const;
+ CompVec surfaceDensities() const;
+ double B (double press,
+ const CompVec& surfvol,
+ PhaseIndex phase) const;
+ double dBdp(double press,
+ const CompVec& surfvol,
+ PhaseIndex phase) const;
+ double R (double press,
+ const CompVec& surfvol,
+ PhaseIndex phase) const;
+ double dRdp(double press,
+ const CompVec& surfvol,
+ PhaseIndex phase) const;
+
+ void getViscosity(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output) const;
+ void B(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output) const;
+ void dBdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output_B,
+ std::vector& output_dBdp) const;
+ void R(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output) const;
+ void dRdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ std::vector& output_R,
+ std::vector& output_dRdp) const;
+
+ private:
+ int region_number_;
+ const MiscibilityProps& propsForPhase(PhaseIndex phase) const;
+
+ boost::scoped_ptr water_props_;
+ boost::scoped_ptr oil_props_;
+ boost::scoped_ptr gas_props_;
+ CompVec densities_;
+ mutable std::vector data1_;
+ mutable std::vector data2_;
+ };
+
+}
+
+
+#endif // OPM_BLACKOILPVT_HEADER_INCLUDED
diff --git a/dune/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp b/dune/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp
new file mode 100644
index 00000000..cd323f06
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp
@@ -0,0 +1,208 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
+#define OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
+
+#include
+#include
+#include
+#include "BlackoilDefs.hpp"
+#include
+#include
+
+namespace Opm
+{
+
+// Forward declaration needed by associated parameters class.
+template
+class FluidMatrixInteractionBlackoil;
+
+template
+class FluidMatrixInteractionBlackoilParams
+{
+public:
+ typedef ScalarT Scalar;
+ void init(const Dune::EclipseGridParser& ep)
+ {
+ // Extract input data.
+ const Dune::SGOF::table_t& sgof_table = ep.getSGOF().sgof_;
+ const Dune::SWOF::table_t& swof_table = ep.getSWOF().swof_;
+ if (sgof_table.size() != 1 || swof_table.size() != 1) {
+ std::cerr << "We must have exactly one SWOF and one SGOF table (at the moment).\n";
+ throw std::logic_error("Not implemented");
+ }
+ const std::vector& sw = swof_table[0][0];
+ const std::vector& krw = swof_table[0][1];
+ const std::vector& krow = swof_table[0][2];
+ const std::vector& pcow = swof_table[0][3];
+ const std::vector& sg = sgof_table[0][0];
+ const std::vector& krg = sgof_table[0][1];
+ const std::vector& krog = sgof_table[0][2];
+ const std::vector& pcog = sgof_table[0][3];
+
+ // Create tables for krw, krow, krg and krog.
+ int samples = 200;
+ buildUniformMonotoneTable(sw, krw, samples, krw_);
+ buildUniformMonotoneTable(sw, krow, samples, krow_);
+ buildUniformMonotoneTable(sg, krg, samples, krg_);
+ buildUniformMonotoneTable(sg, krog, samples, krog_);
+ krocw_ = krow[0]; // At connate water -> ecl. SWOF
+
+ // Create tables for pcow and pcog.
+ buildUniformMonotoneTable(sw, pcow, samples, pcow_);
+ buildUniformMonotoneTable(sg, pcog, samples, pcog_);
+ }
+
+private:
+ template
+ friend class FluidMatrixInteractionBlackoil;
+
+ Dune::utils::UniformTableLinear krw_;
+ Dune::utils::UniformTableLinear krow_;
+ Dune::utils::UniformTableLinear pcow_;
+ Dune::utils::UniformTableLinear krg_;
+ Dune::utils::UniformTableLinear krog_;
+ Dune::utils::UniformTableLinear pcog_;
+ Scalar krocw_; // = krow_(s_wc)
+};
+
+
+/*!
+ * \ingroup material
+ *
+ * \brief Capillary pressures and relative permeabilities for a black oil system.
+ */
+template >
+class FluidMatrixInteractionBlackoil : public BlackoilDefs
+{
+public:
+ typedef ParamsT Params;
+ typedef typename Params::Scalar Scalar;
+
+ /*!
+ * \brief The linear capillary pressure-saturation curve.
+ *
+ * This material law is linear:
+ * \f[
+ p_C = (1 - \overline{S}_w) (p_{C,max} - p_{C,entry}) + p_{C,entry}
+ \f]
+ *
+ * \param Swe Effective saturation of of the wetting phase \f$\overline{S}_w\f$
+ */
+ template
+ static void pC(pcContainerT &pc,
+ const Params ¶ms,
+ const SatContainerT &saturations,
+ Scalar /*temperature*/)
+ {
+ Scalar sw = saturations[Aqua];
+ Scalar sg = saturations[Vapour];
+ pc[Liquid] = 0.0;
+ pc[Aqua] = params.pcow_(sw);
+ pc[Vapour] = params.pcog_(sg);
+ }
+
+ /*!
+ * \brief The saturation-capillary pressure curve.
+ *
+ * This is the inverse of the capillary pressure-saturation curve:
+ * \f[
+ S_w = 1 - \frac{p_C - p_{C,entry}}{p_{C,max} - p_{C,entry}}
+ \f]
+ *
+ * \param pC Capillary pressure \f$\p_C\f$
+ * \return The effective saturaion of the wetting phase \f$\overline{S}_w\f$
+ */
+ template
+ static void S(SatContainerT &saturations,
+ const Params ¶ms,
+ const pcContainerT &pc,
+ Scalar /*temperature*/)
+ {
+ std::cerr << "FluidMatrixInteractionBlackoil::S() is not implemented yet\n";
+ throw std::logic_error("Not implemented");
+ }
+
+
+ /*!
+ * \brief The relative permeability of all phases.
+ */
+ template
+ static void kr(krContainerT &kr,
+ const Params ¶ms,
+ const SatContainerT &saturations,
+ Scalar /*temperature*/)
+ {
+ // Stone-II relative permeability model.
+ Scalar sw = saturations[Aqua];
+ Scalar sg = saturations[Vapour];
+ Scalar krw = params.krw_(sw);
+ Scalar krg = params.krg_(sg);
+ Scalar krow = params.krow_(sw);
+ Scalar krog = params.krog_(sg);
+ Scalar krocw = params.krocw_;
+ kr[Aqua] = krw;
+ kr[Vapour] = krg;
+ kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
+ if (kr[Liquid] < 0.0) {
+ kr[Liquid] = 0.0;
+ }
+ }
+
+
+ /*!
+ * \brief The saturation derivatives of relative permeability of all phases.
+ */
+ template
+ static void dkr(krContainerT &dkr,
+ const Params ¶ms,
+ const SatContainerT &saturations,
+ Scalar /*temperature*/)
+ {
+ for (int p1 = 0; p1 < numPhases; ++p1) {
+ for (int p2 = 0; p2 < numPhases; ++p2) {
+ dkr[p1][p2] = 0.0;
+ }
+ }
+ // Stone-II relative permeability model.
+ Scalar sw = saturations[Aqua];
+ Scalar sg = saturations[Vapour];
+ Scalar krw = params.krw_(sw);
+ Scalar dkrww = params.krw_.derivative(sw);
+ Scalar krg = params.krg_(sg);
+ Scalar dkrgg = params.krg_.derivative(sg);
+ Scalar krow = params.krow_(sw);
+ Scalar dkrow = params.krow_.derivative(sw);
+ Scalar krog = params.krog_(sg);
+ Scalar dkrog = params.krog_.derivative(sg);
+ Scalar krocw = params.krocw_;
+ dkr[Aqua][Aqua] = dkrww;
+ dkr[Vapour][Vapour] = dkrgg;
+ dkr[Liquid][Aqua] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
+ dkr[Liquid][Vapour] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg);
+ }
+};
+
+} // namespace Opm
+
+
+
+
+#endif // OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
diff --git a/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp b/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp
new file mode 100644
index 00000000..dfc8a17b
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp
@@ -0,0 +1,91 @@
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
+#define OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
+
+
+#include "BlackoilDefs.hpp"
+
+
+namespace Opm
+{
+ /*!
+ * \brief Fluid state for a black oil model.
+ */
+ struct FluidStateBlackoil : public BlackoilDefs
+ {
+ Scalar temperature_;
+ CompVec surface_volume_;
+ PhaseVec phase_pressure_;
+ PhaseVec phase_volume_density_;
+ Scalar total_phase_volume_density_;
+ PhaseVec formation_volume_factor_;
+ PhaseVec solution_factor_;
+ PhaseToCompMatrix phase_to_comp_; // RB^{-1} in Fortran ordering
+ PhaseVec saturation_;
+ PhaseVec phase_compressibility_;
+ Scalar total_compressibility_;
+ Scalar experimental_term_;
+ PhaseVec viscosity_;
+ PhaseVec relperm_;
+ PhaseJacobian drelperm_;
+ PhaseVec mobility_;
+ PhaseJacobian dmobility_;
+ };
+
+
+ /*!
+ * \brief Multiple fluid states for a black oil model.
+ */
+ struct AllFluidStates : public BlackoilDefs
+ {
+ // Canonical state variables.
+ std::vector surface_volume_density; // z
+ std::vector phase_pressure; // p
+
+ // Variables from PVT functions.
+ std::vector formation_volume_factor; // B
+ std::vector formation_volume_factor_deriv; // dB/dp
+ std::vector solution_factor; // R
+ std::vector solution_factor_deriv; // dR/dp
+ std::vector viscosity; // mu
+
+ // Variables computed from PVT data.
+ // The A matrices are all in Fortran order (or, equivalently,
+ // we store the transposes).
+ std::vector state_matrix; // A' = (RB^{-1})'
+ std::vector phase_volume_density; // u
+ std::vector total_phase_volume_density; // sum(u)
+ std::vector saturation; // s = u/sum(u)
+ std::vector phase_compressibility; // c
+ std::vector total_compressibility; // cT
+ std::vector experimental_term; // ex = sum(Ai*dA*Ai*z)
+
+ // Variables computed from saturation.
+ std::vector relperm; // kr
+ std::vector relperm_deriv; // dkr/ds
+ std::vector mobility; // lambda
+ std::vector mobility_deriv; // dlambda/ds
+ };
+
+} // end namespace Opm
+
+
+#endif // OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
diff --git a/dune/porsol/blackoil/fluid/Makefile.am b/dune/porsol/blackoil/fluid/Makefile.am
new file mode 100644
index 00000000..efec558f
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/Makefile.am
@@ -0,0 +1,11 @@
+noinst_LTLIBRARIES = libblackoilfluid_noinst.la
+
+libblackoilfluid_noinst_la_SOURCES = \
+ BlackoilPVT.cpp \
+ MiscibilityDead.cpp \
+ MiscibilityLiveGas.cpp \
+ MiscibilityLiveOil.cpp \
+ MiscibilityProps.cpp
+
+
+include $(top_srcdir)/am/global-rules
diff --git a/dune/porsol/blackoil/fluid/MiscibilityDead.cpp b/dune/porsol/blackoil/fluid/MiscibilityDead.cpp
new file mode 100644
index 00000000..c92b2abf
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/MiscibilityDead.cpp
@@ -0,0 +1,179 @@
+//===========================================================================
+//
+// File: MiscibilityDead.cpp
+//
+// Created: Wed Feb 10 09:06:04 2010
+//
+// Author: Bjørn Spjelkavik
+//
+// Revision: $Id$
+//
+//===========================================================================
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#include
+#include "MiscibilityDead.hpp"
+#include
+#include
+#include
+#include
+#include
+#include
+
+using namespace std;
+using namespace Dune;
+
+namespace Opm
+{
+
+ //------------------------------------------------------------------------
+ // Member functions
+ //-------------------------------------------------------------------------
+
+ /// Constructor
+ MiscibilityDead::MiscibilityDead(const table_t& pvd_table)
+ {
+ const int region_number = 0;
+ if (pvd_table.size() != 1) {
+ THROW("More than one PVT-region");
+ }
+
+ // Copy data
+ const int sz = pvd_table[region_number][0].size();
+ std::vector press(sz);
+ std::vector B_inv(sz);
+ std::vector visc(sz);
+ for (int i = 0; i < sz; ++i) {
+ press[i] = pvd_table[region_number][0][i];
+ B_inv[i] = 1.0 / pvd_table[region_number][1][i];
+ visc[i] = pvd_table[region_number][2][i];
+ }
+ int samples = 1025;
+ buildUniformMonotoneTable(press, B_inv, samples, one_over_B_);
+ buildUniformMonotoneTable(press, visc, samples, viscosity_);
+
+ // Dumping the created tables.
+// static int count = 0;
+// std::ofstream os((std::string("dump-") + boost::lexical_cast(count++)).c_str());
+// os.precision(15);
+// os << "1/B\n\n" << one_over_B_
+// << "\n\nvisc\n\n" << viscosity_ << std::endl;
+ }
+
+ // Destructor
+ MiscibilityDead::~MiscibilityDead()
+ {
+ }
+
+ double MiscibilityDead::getViscosity(int region, double press, const surfvol_t& /*surfvol*/) const
+ {
+ return viscosity_(press);
+ }
+
+ void MiscibilityDead::getViscosity(const std::vector& pressures,
+ const std::vector&,
+ int phase,
+ std::vector& output) const
+ {
+ int num = pressures.size();
+ output.resize(num);
+#pragma omp parallel for
+ for (int i = 0; i < num; ++i) {
+ output[i] = viscosity_(pressures[i][phase]);
+ }
+ }
+
+ double MiscibilityDead::B(int region, double press, const surfvol_t& /*surfvol*/) const
+ {
+ // Interpolate 1/B
+ return 1.0/one_over_B_(press);
+ }
+
+ void MiscibilityDead::B(const std::vector& pressures,
+ const std::vector&,
+ int phase,
+ std::vector& output) const
+ {
+ int num = pressures.size();
+ output.resize(num);
+#pragma omp parallel for
+ for (int i = 0; i < num; ++i) {
+ output[i] = 1.0/one_over_B_(pressures[i][phase]);
+ }
+ }
+
+ double MiscibilityDead::dBdp(int region, double press, const surfvol_t& /*surfvol*/) const
+ {
+ // Interpolate 1/B
+ surfvol_t dummy_surfvol;
+ double Bg = B(region, press, dummy_surfvol);
+ return -Bg*Bg*one_over_B_.derivative(press);
+ }
+
+ void MiscibilityDead::dBdp(const std::vector& pressures,
+ const std::vector& surfvols,
+ int phase,
+ std::vector& output_B,
+ std::vector& output_dBdp) const
+ {
+ B(pressures, surfvols, phase, output_B);
+ int num = pressures.size();
+ output_dBdp.resize(num);
+#pragma omp parallel for
+ for (int i = 0; i < num; ++i) {
+ double Bg = output_B[i];
+ output_dBdp[i] = -Bg*Bg*one_over_B_.derivative(pressures[i][phase]);
+ }
+ }
+
+ double MiscibilityDead::R(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
+ {
+ return 0.0;
+ }
+
+ void MiscibilityDead::R(const std::vector& pressures,
+ const std::vector&,
+ int,
+ std::vector& output) const
+ {
+ int num = pressures.size();
+ output.clear();
+ output.resize(num, 0.0);
+ }
+
+ double MiscibilityDead::dRdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
+ {
+ return 0.0;
+ }
+
+ void MiscibilityDead::dRdp(const std::vector& pressures,
+ const std::vector&,
+ int,
+ std::vector& output_R,
+ std::vector& output_dRdp) const
+ {
+ int num = pressures.size();
+ output_R.clear();
+ output_R.resize(num, 0.0);
+ output_dRdp.clear();
+ output_dRdp.resize(num, 0.0);
+ }
+
+}
diff --git a/dune/porsol/blackoil/fluid/MiscibilityDead.hpp b/dune/porsol/blackoil/fluid/MiscibilityDead.hpp
new file mode 100644
index 00000000..8962bd58
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/MiscibilityDead.hpp
@@ -0,0 +1,89 @@
+//===========================================================================
+//
+// File: MiscibilityDead.hpp
+//
+// Created: Wed Feb 10 09:05:47 2010
+//
+// Author: Bjørn Spjelkavik
+//
+// Revision: $Id$
+//
+//===========================================================================
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef SINTEF_MISCIBILITYDEAD_HEADER
+#define SINTEF_MISCIBILITYDEAD_HEADER
+
+/** Class for immiscible dead oil and dry gas.
+ * Detailed description.
+ */
+
+#include "MiscibilityProps.hpp"
+#include
+
+namespace Opm
+{
+ class MiscibilityDead : public MiscibilityProps
+ {
+ public:
+ typedef std::vector > > table_t;
+
+ MiscibilityDead(const table_t& pvd_table);
+ virtual ~MiscibilityDead();
+
+ virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
+ virtual double B(int region, double press, const surfvol_t& surfvol) const;
+ virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
+ virtual double R(int region, double press, const surfvol_t& surfvol) const;
+ virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
+
+ virtual void getViscosity(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const;
+ virtual void B(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const;
+ virtual void dBdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output_B,
+ std::vector& output_dBdp) const;
+ virtual void R(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const;
+ virtual void dRdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output_R,
+ std::vector& output_dRdp) const;
+
+ private:
+ // PVT properties of dry gas or dead oil
+ Dune::utils::UniformTableLinear one_over_B_;
+ Dune::utils::UniformTableLinear viscosity_;
+ };
+
+}
+
+#endif // SINTEF_MISCIBILITYDEAD_HEADER
+
diff --git a/dune/porsol/blackoil/fluid/MiscibilityLiveGas.cpp b/dune/porsol/blackoil/fluid/MiscibilityLiveGas.cpp
new file mode 100644
index 00000000..3d44239a
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/MiscibilityLiveGas.cpp
@@ -0,0 +1,286 @@
+//===========================================================================
+//
+// File: MiscibilityLiveGas.cpp
+//
+// Created: Wed Feb 10 09:21:53 2010
+//
+// Author: Bjørn Spjelkavik
+//
+// Revision: $Id$
+//
+//===========================================================================
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+
+#include
+#include "MiscibilityLiveGas.hpp"
+#include
+#include
+
+using namespace std;
+using namespace Dune;
+
+namespace Opm
+{
+
+ //------------------------------------------------------------------------
+ // Member functions
+ //-------------------------------------------------------------------------
+
+ /// Constructor
+ MiscibilityLiveGas::MiscibilityLiveGas(const table_t& pvtg)
+ {
+ // GAS, PVTG
+ const int region_number = 0;
+ if (pvtg.size() != 1) {
+ THROW("More than one PVD-region");
+ }
+ saturated_gas_table_.resize(4);
+ const int sz = pvtg[region_number].size();
+ for (int k=0; k<4; ++k) {
+ saturated_gas_table_[k].resize(sz);
+ }
+
+ for (int i=0; i& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const
+ {
+ ASSERT(pressures.size() == surfvol.size());
+ int num = pressures.size();
+ output.resize(num);
+ for (int i = 0; i < num; ++i) {
+ output[i] = miscible_gas(pressures[i][phase], surfvol[i], 2, false);
+ }
+ }
+
+ // Vaporised oil-gas ratio.
+ double MiscibilityLiveGas::R(int /*region*/, double press, const surfvol_t& surfvol) const
+ {
+ if (surfvol[Liquid] == 0.0) {
+ return 0.0;
+ }
+ double R = linearInterpolationExtrap(saturated_gas_table_[0],
+ saturated_gas_table_[3], press);
+ double maxR = surfvol[Liquid]/surfvol[Vapour];
+ if (R < maxR ) { // Saturated case
+ return R;
+ } else {
+ return maxR; // Undersaturated case
+ }
+ }
+
+ void MiscibilityLiveGas::R(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const
+ {
+ ASSERT(pressures.size() == surfvol.size());
+ int num = pressures.size();
+ output.resize(num);
+ for (int i = 0; i < num; ++i) {
+ output[i] = R(0, pressures[i][phase], surfvol[i]);
+ }
+ }
+
+ // Vaporised oil-gas ratio derivative
+ double MiscibilityLiveGas::dRdp(int /*region*/, double press, const surfvol_t& surfvol) const
+ {
+ double R = linearInterpolationExtrap(saturated_gas_table_[0],
+ saturated_gas_table_[3], press);
+ double maxR = surfvol[Liquid]/surfvol[Vapour];
+ if (R < maxR ) { // Saturated case
+ return linearInterpolDerivative(saturated_gas_table_[0],
+ saturated_gas_table_[3],
+ press);
+ } else {
+ return 0.0; // Undersaturated case
+ }
+ }
+
+ void MiscibilityLiveGas::dRdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output_R,
+ std::vector& output_dRdp) const
+ {
+ ASSERT(pressures.size() == surfvol.size());
+ R(pressures, surfvol, phase, output_R);
+ int num = pressures.size();
+ output_dRdp.resize(num);
+ for (int i = 0; i < num; ++i) {
+ output_dRdp[i] = dRdp(0, pressures[i][phase], surfvol[i]); // \TODO Speedup here by using already evaluated R.
+ }
+ }
+
+ double MiscibilityLiveGas::B(int /*region*/, double press, const surfvol_t& surfvol) const
+ {
+ if (surfvol[Vapour] == 0.0) return 1.0; // To handle no-gas case.
+ return miscible_gas(press, surfvol, 1, false);
+ }
+
+ void MiscibilityLiveGas::B(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const
+ {
+ ASSERT(pressures.size() == surfvol.size());
+ int num = pressures.size();
+ output.resize(num);
+ for (int i = 0; i < num; ++i) {
+ output[i] = B(0, pressures[i][phase], surfvol[i]);
+ }
+ }
+
+ double MiscibilityLiveGas::dBdp(int /*region*/, double press, const surfvol_t& surfvol) const
+ {
+ if (surfvol[Vapour] == 0.0) return 0.0; // To handle no-gas case.
+ return miscible_gas(press, surfvol, 1, true);
+ }
+
+ void MiscibilityLiveGas::dBdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output_B,
+ std::vector& output_dBdp) const
+ {
+ ASSERT(pressures.size() == surfvol.size());
+ B(pressures, surfvol, phase, output_B);
+ int num = pressures.size();
+ output_dBdp.resize(num);
+ for (int i = 0; i < num; ++i) {
+ output_dBdp[i] = dBdp(0, pressures[i][phase], surfvol[i]); // \TODO Speedup here by using already evaluated B.
+ }
+ }
+
+ double MiscibilityLiveGas::miscible_gas(double press, const surfvol_t& surfvol, int item,
+ bool deriv) const
+ {
+ int section;
+ double R = linearInterpolationExtrap(saturated_gas_table_[0],
+ saturated_gas_table_[3], press,
+ section);
+ double maxR = surfvol[Liquid]/surfvol[Vapour];
+ if (deriv) {
+ if (R < maxR ) { // Saturated case
+ return linearInterpolDerivative(saturated_gas_table_[0],
+ saturated_gas_table_[item],
+ press);
+ } else { // Undersaturated case
+ int is = section;
+ if (undersat_gas_tables_[is][0].size() < 2) {
+ double val = (saturated_gas_table_[item][is+1]
+ - saturated_gas_table_[item][is]) /
+ (saturated_gas_table_[0][is+1] -
+ saturated_gas_table_[0][is]);
+ return val;
+ }
+ double val1 =
+ linearInterpolationExtrap(undersat_gas_tables_[is][0],
+ undersat_gas_tables_[is][item],
+ maxR);
+ double val2 =
+ linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
+ undersat_gas_tables_[is+1][item],
+ maxR);
+ double val = (val2 - val1)/
+ (saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
+ return val;
+ }
+ } else {
+ if (R < maxR ) { // Saturated case
+ return linearInterpolationExtrap(saturated_gas_table_[0],
+ saturated_gas_table_[item],
+ press);
+ } else { // Undersaturated case
+ int is = section;
+ // Extrapolate from first table section
+ if (is == 0 && press < saturated_gas_table_[0][0]) {
+ return linearInterpolationExtrap(undersat_gas_tables_[0][0],
+ undersat_gas_tables_[0][item],
+ maxR);
+ }
+
+ // Extrapolate from last table section
+ int ltp = saturated_gas_table_[0].size() - 1;
+ if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
+ return linearInterpolationExtrap(undersat_gas_tables_[ltp][0],
+ undersat_gas_tables_[ltp][item],
+ maxR);
+ }
+
+ // Interpolate between table sections
+ double w = (press - saturated_gas_table_[0][is]) /
+ (saturated_gas_table_[0][is+1] -
+ saturated_gas_table_[0][is]);
+ if (undersat_gas_tables_[is][0].size() < 2) {
+ double val = saturated_gas_table_[item][is] +
+ w*(saturated_gas_table_[item][is+1] -
+ saturated_gas_table_[item][is]);
+ return val;
+ }
+ double val1 =
+ linearInterpolationExtrap(undersat_gas_tables_[is][0],
+ undersat_gas_tables_[is][item],
+ maxR);
+ double val2 =
+ linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
+ undersat_gas_tables_[is+1][item],
+ maxR);
+ double val = val1 + w*(val2 - val1);
+ return val;
+ }
+ }
+ }
+
+
+} // namespace Opm
diff --git a/dune/porsol/blackoil/fluid/MiscibilityLiveGas.hpp b/dune/porsol/blackoil/fluid/MiscibilityLiveGas.hpp
new file mode 100644
index 00000000..801649ed
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/MiscibilityLiveGas.hpp
@@ -0,0 +1,92 @@
+//===========================================================================
+//
+// File: MiscibilityLiveGas.hpp
+//
+// Created: Wed Feb 10 09:21:26 2010
+//
+// Author: Bjørn Spjelkavik
+//
+// Revision: $Id$
+//
+//===========================================================================
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#ifndef SINTEF_MISCIBILITYLIVEGAS_HEADER
+#define SINTEF_MISCIBILITYLIVEGAS_HEADER
+
+ /** Class for miscible wet gas.
+ * Detailed description.
+ */
+
+#include "MiscibilityProps.hpp"
+
+namespace Opm
+{
+ class MiscibilityLiveGas : public MiscibilityProps
+ {
+ public:
+ typedef std::vector > > table_t;
+
+ MiscibilityLiveGas(const table_t& pvto);
+ virtual ~MiscibilityLiveGas();
+
+ virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
+ virtual double R(int region, double press, const surfvol_t& surfvol) const;
+ virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
+ virtual double B(int region, double press, const surfvol_t& surfvol) const;
+ virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
+
+ virtual void getViscosity(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const;
+ virtual void B(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const;
+ virtual void dBdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output_B,
+ std::vector& output_dBdp) const;
+ virtual void R(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output) const;
+ virtual void dRdp(const std::vector& pressures,
+ const std::vector& surfvol,
+ int phase,
+ std::vector& output_R,
+ std::vector& output_dRdp) const;
+
+ protected:
+ // item: 1=B 2=mu;
+ double miscible_gas(double press, const surfvol_t& surfvol, int item,
+ bool deriv = false) const;
+ // PVT properties of wet gas (with vaporised oil)
+ std::vector > saturated_gas_table_;
+ std::vector > > undersat_gas_tables_;
+
+ };
+
+}
+
+#endif // SINTEF_MISCIBILITYLIVEGAS_HEADER
+
diff --git a/dune/porsol/blackoil/fluid/MiscibilityLiveOil.cpp b/dune/porsol/blackoil/fluid/MiscibilityLiveOil.cpp
new file mode 100644
index 00000000..4b9fbe33
--- /dev/null
+++ b/dune/porsol/blackoil/fluid/MiscibilityLiveOil.cpp
@@ -0,0 +1,395 @@
+//===========================================================================
+//
+// File: MiscibiltyLiveOil.cpp
+//
+// Created: Wed Feb 10 09:08:25 2010
+//
+// Author: Bjørn Spjelkavik
+//
+// Revision: $Id$
+//
+//===========================================================================
+/*
+ Copyright 2010 SINTEF ICT, Applied Mathematics.
+
+ This file is part of the Open Porous Media project (OPM).
+
+ OPM is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OPM is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OPM. If not, see .
+*/
+
+#include
+#include "MiscibilityLiveOil.hpp"
+#include
+#include
+
+using namespace std;
+using namespace Dune;
+
+namespace Opm
+{
+
+
+ //------------------------------------------------------------------------
+ // Member functions
+ //-------------------------------------------------------------------------
+
+ /// Constructor
+ MiscibilityLiveOil::MiscibilityLiveOil(const table_t& pvto)
+ {
+ // OIL, PVTO
+ const int region_number = 0;
+ if (pvto.size() != 1) {
+ THROW("More than one PVD-region");
+ }
+ saturated_oil_table_.resize(4);
+ const int sz = pvto[region_number].size();
+ for (int k=0; k<4; ++k) {
+ saturated_oil_table_[k].resize(sz);
+ }
+ for (int i=0; i 1) {
+ iPrev = i;
+ continue;
+ }
+
+ bool flagPrev = (iPrev >= 0);
+ bool flagNext = true;
+ if (iNext < i) {
+ iPrev = iNext;
+ flagPrev = true;
+ iNext = i+1;
+ while (undersat_oil_tables_[iNext][0].size() < 2) {
+ ++iNext;
+ }
+ }
+ double slopePrevBinv = 0.0;
+ double slopePrevVisc = 0.0;
+ double slopeNextBinv = 0.0;
+ double slopeNextVisc = 0.0;
+ while (flagPrev || flagNext) {
+ double pressure0 = undersat_oil_tables_[i][0].back();
+ double pressure = 1.0e47;
+ if (flagPrev) {
+ std::vector