OilfieldToolsNet/opm-upscaling
4
0
Fork 0
mirror of https://github.com/OPM/opm-upscaling.git synced 2025-02-25 18:45:23 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #228 from atgeirr/more-removals

Browse Source 
Deleting more files not needed after the impes sims were removed.
This commit is contained in:
Arne Morten Kvarving 2017-11-16 14:45:28 +01:00 committed by GitHub
commit fb90035d0e
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
27 changed files with 0 additions and 32425 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
CMakeLists_files.cmake
View File

@ -32,11 +32,6 @@ list (APPEND TEST_SOURCE_FILES
)
list(APPEND MAIN_SOURCE_FILES
opm/porsol/blackoil/fluid/BlackoilPVT.cpp
opm/porsol/blackoil/fluid/MiscibilityDead.cpp
opm/porsol/blackoil/fluid/MiscibilityLiveGas.cpp
opm/porsol/blackoil/fluid/MiscibilityLiveOil.cpp
opm/porsol/blackoil/fluid/MiscibilityProps.cpp
opm/porsol/common/blas_lapack.cpp
opm/porsol/common/BoundaryPeriodicity.cpp
opm/porsol/common/ImplicitTransportDefs.cpp
@ -91,7 +86,6 @@ list (APPEND TEST_DATA_FILES
list (APPEND EXAMPLE_SOURCE_FILES
examples/aniso_implicitcap_test.cpp
examples/aniso_simulator_test.cpp
examples/co2_blackoil_pvt.cpp
examples/cpchop.cpp
examples/cpchop_depthtrend.cpp
examples/cpregularize.cpp
@ -153,23 +147,6 @@ list (APPEND PROGRAM_SOURCE_FILES
# originally generated with the command:
# find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort
list (APPEND PUBLIC_HEADER_FILES
opm/porsol/blackoil/BlackoilFluid.hpp
opm/porsol/blackoil/BlackoilInitialization.hpp
opm/porsol/blackoil/BlackoilSimulator.hpp
opm/porsol/blackoil/BlackoilWells.hpp
opm/porsol/blackoil/co2fluid/benchmark3co2tables.hh
opm/porsol/blackoil/co2fluid/BlackoilCo2PVT.hpp
opm/porsol/blackoil/ComponentTransport.hpp
opm/porsol/blackoil/fluid/BlackoilComponent.hpp
opm/porsol/blackoil/fluid/BlackoilDefs.hpp
opm/porsol/blackoil/fluid/BlackoilPVT.hpp
opm/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp
opm/porsol/blackoil/fluid/FluidStateBlackoil.hpp
opm/porsol/blackoil/fluid/MiscibilityDead.hpp
opm/porsol/blackoil/fluid/MiscibilityLiveGas.hpp
opm/porsol/blackoil/fluid/MiscibilityLiveOil.hpp
opm/porsol/blackoil/fluid/MiscibilityProps.hpp
opm/porsol/blackoil/fluid/MiscibilityWater.hpp
opm/porsol/common/BCRSMatrixBlockAssembler.hpp
opm/porsol/common/blas_lapack.hpp
opm/porsol/common/BoundaryConditions.hpp
@ -211,8 +188,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/porsol/mimetic/IncompFlowSolverHybrid.hpp
opm/porsol/mimetic/MimeticIPAnisoRelpermEvaluator.hpp
opm/porsol/mimetic/MimeticIPEvaluator.hpp
opm/porsol/mimetic/TpfaCompressibleAssembler.hpp
opm/porsol/mimetic/TpfaCompressible.hpp
opm/upscaling/ParserAdditions.hpp
opm/upscaling/SinglePhaseUpscaler.hpp
opm/upscaling/SteadyStateUpscaler.hpp

50
examples/co2_blackoil_pvt.cpp
View File

@ -1,50 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/utility/StopWatch.hpp>
#include <opm/porsol/blackoil/co2fluid/BlackoilCo2PVT.hpp>
#include <iostream>
int main(int argc, char** argv)
try
{
double temperature = 300.;
if (argc == 2) {
temperature = std::atof(argv[1]);
}
Opm::BlackoilCo2PVT boPvt;
Opm::time::StopWatch clock;
clock.start();
boPvt.init( Opm::Deck{} );
boPvt.generateBlackOilTables(temperature);
clock.stop();
std::cout << "\n\nInitialisation and table generation - clock time (secs): "
<< clock.secsSinceStart() << std::endl;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

584
opm/porsol/blackoil/BlackoilFluid.hpp
View File

@ -1,584 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILFLUID_HEADER_INCLUDED
#define OPM_BLACKOILFLUID_HEADER_INCLUDED
#include <opm/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp>
#include <opm/porsol/blackoil/fluid/FluidStateBlackoil.hpp>
#include <opm/porsol/blackoil/fluid/BlackoilPVT.hpp>
#include <dune/common/fvector.hh>
#include <vector>
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 Opm::Deck& deck)
{
fmi_params_.init(deck);
// FluidSystemBlackoil<>::init(parser);
pvt_.init(deck);
const auto& densityRecord = deck.getKeyword("DENSITY").getRecord(0);
surface_densities_[Oil] = densityRecord.getItem("OIL").getSIDouble(0);
surface_densities_[Water] = densityRecord.getItem("WATER").getSIDouble(0);
surface_densities_[Gas] = densityRecord.getItem("GAS").getSIDouble(0);
}
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<double>::kr(state.relperm_, fmi_params_, state.saturation_, state.temperature_);
FluidMatrixInteractionBlackoil<double>::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 <class States>
void computeBAndR(States& states) const
{
const std::vector<PhaseVec>& p = states.phase_pressure;
const std::vector<CompVec>& z = states.surface_volume_density;
std::vector<PhaseVec>& B = states.formation_volume_factor;
std::vector<PhaseVec>& R = states.solution_factor;
pvt_.B(p, z, B);
pvt_.R(p, z, R);
}
/// Input: p, z
/// Output: B, R, mu
template <class States>
void computePvtNoDerivs(States& states) const
{
computeBAndR(states);
const std::vector<PhaseVec>& p = states.phase_pressure;
const std::vector<CompVec>& z = states.surface_volume_density;
std::vector<PhaseVec>& mu = states.viscosity;
pvt_.getViscosity(p, z, mu);
}
/// Input: p, z
/// Output: B, dB/dp, R, dR/dp, mu
template <class States>
void computePvt(States& states) const
{
const std::vector<PhaseVec>& p = states.phase_pressure;
const std::vector<CompVec>& z = states.surface_volume_density;
std::vector<PhaseVec>& B = states.formation_volume_factor;
std::vector<PhaseVec>& dB = states.formation_volume_factor_deriv;
std::vector<PhaseVec>& R = states.solution_factor;
std::vector<PhaseVec>& dR = states.solution_factor_deriv;
std::vector<PhaseVec>& 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 <class States>
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 <class States>
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 <class States>
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<double>::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 <class States>
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<double>::kr(kr, fmi_params_, s, 300.0);
FluidMatrixInteractionBlackoil<double>::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<double> 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 <class FluidState>
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<CompVec> surface_volume_density; // z
std::vector<PhaseVec> phase_pressure; // p
// Variables from PVT functions.
std::vector<PhaseVec> formation_volume_factor; // B
std::vector<PhaseVec> solution_factor; // R
// Variables computed from PVT data.
// The A matrices are all in Fortran order (or, equivalently,
// we store the transposes).
std::vector<PhaseToCompMatrix> state_matrix; // A' = (RB^{-1})'
// Variables computed from saturation.
std::vector<PhaseVec> mobility; // lambda
std::vector<PhaseJacobian> mobility_deriv; // dlambda/ds
// Gravity and/or capillary pressure potential differences.
std::vector<PhaseVec> gravity_potential; // (\rho g \delta z)-ish contribution per face
};
struct AllFluidData : public BlackoilDefs
{
// Per-cell data
AllFluidStates cell_data;
std::vector<double> voldiscr;
std::vector<double> relvoldiscr;
// Per-face data.
FaceFluidData face_data;
public:
template <class Grid, class Rock>
void computeNew(const Grid& grid,
const Rock& rock,
const BlackoilFluid& fluid,
const typename Grid::Vector gravity,
const std::vector<PhaseVec>& cell_pressure,
const std::vector<PhaseVec>& face_pressure,
const std::vector<CompVec>& 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;
static_assert(np == nc, "");
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 <class Grid>
void computeUpwindProperties(const Grid& grid,
const BlackoilFluid& fluid,
const typename Grid::Vector gravity,
const std::vector<PhaseVec>& cell_pressure,
const std::vector<PhaseVec>& face_pressure,
const std::vector<CompVec>& 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

330
opm/porsol/blackoil/BlackoilInitialization.hpp
View File

@ -1,330 +0,0 @@
/*
Copyright 2011 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILINITIALIZATION_HEADER_INCLUDED
#define OPM_BLACKOILINITIALIZATION_HEADER_INCLUDED
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <iostream>
namespace Opm
{
/// Base class for initialization codes.
template <class Simulator>
class BlackoilInitialization
{
public:
typedef typename Simulator::State State;
typedef typename Simulator::Grid Grid;
typedef typename Simulator::Fluid Fluid;
virtual void init(const Opm::ParameterGroup& param,
const Grid& grid,
const Fluid& fluid,
typename Grid::Vector gravity,
State& simstate) = 0;
};
/// Initialize basic cases.
template <class Simulator>
class BasicInitialization : public BlackoilInitialization<Simulator>
{
public:
typedef typename Simulator::State State;
typedef typename Simulator::Grid Grid;
typedef typename Simulator::Fluid Fluid;
virtual void init(const Opm::ParameterGroup& param,
const Grid& grid,
const Fluid& fluid,
typename Grid::Vector gravity,
State& simstate)
{
typedef typename Fluid::CompVec CompVec;
typedef typename Fluid::PhaseVec PhaseVec;
if (param.getDefault("heterogenous_initial_mix", false)) {
CompVec init_oil(0.0);
init_oil[Fluid::Oil] = 1.0;
CompVec init_water(0.0);
init_water[Fluid::Water] = 1.0;
simstate.cell_z_.resize(grid.numCells());
std::fill(simstate.cell_z_.begin(), simstate.cell_z_.begin() + simstate.cell_z_.size()/2, init_oil);
std::fill(simstate.cell_z_.begin() + simstate.cell_z_.size()/2, simstate.cell_z_.end(), init_water);
OPM_MESSAGE("******* Assuming zero capillary pressures *******");
PhaseVec init_p(100.0*Opm::unit::barsa);
simstate.cell_pressure_.resize(grid.numCells(), init_p);
// if (gravity.two_norm() != 0.0) {
// double ref_gravpot = grid.cellCentroid(0)*gravity;
// double rho = init_z*fluid_.surfaceDensities(); // Assuming incompressible, and constant initial z.
// for (int cell = 1; cell < grid.numCells(); ++cell) {
// double press = rho*(grid.cellCentroid(cell)*gravity - ref_gravpot) + simstate.cell_pressure_[0][0];
// simstate.cell_pressure_[cell] = PhaseVec(press);
// }
// }
} else if (param.getDefault("unstable_initial_mix", false)) {
CompVec init_oil(0.0);
init_oil[Fluid::Oil] = 1.0;
init_oil[Fluid::Gas] = 0.0;
CompVec init_water(0.0);
init_water[Fluid::Water] = 1.0;
CompVec init_gas(0.0);
init_gas[Fluid::Gas] = 150.0;
simstate.cell_z_.resize(grid.numCells());
std::fill(simstate.cell_z_.begin(),
simstate.cell_z_.begin() + simstate.cell_z_.size()/3,
init_water);
std::fill(simstate.cell_z_.begin() + simstate.cell_z_.size()/3,
simstate.cell_z_.begin() + 2*(simstate.cell_z_.size()/3),
init_oil);
std::fill(simstate.cell_z_.begin() + 2*(simstate.cell_z_.size()/3),
simstate.cell_z_.end(),
init_gas);
OPM_MESSAGE("******* Assuming zero capillary pressures *******");
PhaseVec init_p(100.0*Opm::unit::barsa);
simstate.cell_pressure_.resize(grid.numCells(), init_p);
if (gravity.two_norm() != 0.0) {
typename Fluid::FluidState state = fluid.computeState(simstate.cell_pressure_[0], simstate.cell_z_[0]);
simstate.cell_z_[0] *= 1.0/state.total_phase_volume_density_;
for (int cell = 1; cell < grid.numCells(); ++cell) {
double fluid_vol_dens;
int cnt =0;
do {
double rho = 0.5*((simstate.cell_z_[cell]+simstate.cell_z_[cell-1])*fluid.surfaceDensities());
double press = rho*((grid.cellCentroid(cell) - grid.cellCentroid(cell-1))*gravity) + simstate.cell_pressure_[cell-1][0];
simstate.cell_pressure_[cell] = PhaseVec(press);
state = fluid.computeState(simstate.cell_pressure_[cell], simstate.cell_z_[cell]);
fluid_vol_dens = state.total_phase_volume_density_;
simstate.cell_z_[cell] *= 1.0/fluid_vol_dens;
++cnt;
} while (std::fabs((fluid_vol_dens-1.0)) > 1.0e-8 && cnt < 10);
}
} else {
std::cout << "---- Exit - BlackoilSimulator.hpp: No gravity, no fun ... ----" << std::endl;
exit(-1);
}
} else if (param.getDefault("CO2-injection", false)) {
CompVec init_water(0.0);
// Initially water filled (use Oil-component for water in order
// to utilise blackoil mechanisms for brine-co2 interaction)
init_water[Fluid::Oil] = 1.0;
simstate.cell_z_.resize(grid.numCells());
std::fill(simstate.cell_z_.begin(),simstate.cell_z_.end(),init_water);
double datum_pressure_barsa = param.getDefault<double>("datum_pressure", 200.0);
double datum_pressure = Opm::unit::convert::from(datum_pressure_barsa, Opm::unit::barsa);
PhaseVec init_p(datum_pressure);
simstate.cell_pressure_.resize(grid.numCells(), init_p);
// Simple initial condition based on "incompressibility"-assumption
double zMin = grid.cellCentroid(0)[2];
for (int cell = 1; cell < grid.numCells(); ++cell) {
if (grid.cellCentroid(cell)[2] < zMin)
zMin = grid.cellCentroid(cell)[2];
}
typename Fluid::FluidState state = fluid.computeState(init_p, init_water);
simstate.cell_z_[0] *= 1.0/state.total_phase_volume_density_;
double density = (init_water*fluid.surfaceDensities())/state.total_phase_volume_density_;
for (int cell = 0; cell < grid.numCells(); ++cell) {
double pressure(datum_pressure + (grid.cellCentroid(cell)[2] - zMin)*gravity[2]*density);
simstate.cell_pressure_[cell] = PhaseVec(pressure);
state = fluid.computeState(simstate.cell_pressure_[cell], simstate.cell_z_[cell]);
simstate.cell_z_[cell] *= 1.0/state.total_phase_volume_density_;
}
} else {
CompVec init_z(0.0);
double initial_mixture_gas = param.getDefault("initial_mixture_gas", 0.0);
double initial_mixture_oil = param.getDefault("initial_mixture_oil", 1.0);
double initial_mixture_water = param.getDefault("initial_mixture_water", 0.0);
init_z[Fluid::Water] = initial_mixture_water;
init_z[Fluid::Gas] = initial_mixture_gas;
init_z[Fluid::Oil] = initial_mixture_oil;
simstate.cell_z_.resize(grid.numCells(), init_z);
OPM_MESSAGE("******* Assuming zero capillary pressures *******");
PhaseVec init_p(param.getDefault("initial_pressure", 100.0*Opm::unit::barsa));
simstate.cell_pressure_.resize(grid.numCells(), init_p);
if (gravity.two_norm() != 0.0) {
double ref_gravpot = grid.cellCentroid(0)*gravity;
double rho = init_z*fluid.surfaceDensities(); // Assuming incompressible, and constant initial z.
for (int cell = 1; cell < grid.numCells(); ++cell) {
double press = rho*(grid.cellCentroid(cell)*gravity - ref_gravpot) + simstate.cell_pressure_[0][0];
simstate.cell_pressure_[cell] = PhaseVec(press);
}
}
}
}
};
/// Initialize SPE9 type case.
template <class Simulator>
class SPE9Initialization : public BlackoilInitialization<Simulator>
{
public:
typedef typename Simulator::State State;
typedef typename Simulator::Grid Grid;
typedef typename Simulator::Fluid Fluid;
virtual void init(const Opm::ParameterGroup& param,
const Grid& grid,
const Fluid& fluid,
typename Grid::Vector gravity,
State& simstate)
{
typedef typename Fluid::CompVec CompVec;
double zeroDepth = param.getDefault("zero_depth", 2743.2);
int nx = param.getDefault<int>("nx", 24);
int ny = param.getDefault<int>("ny", 25);
int nz = param.getDefault<int>("nz", 15);
// double datum_depth = param.getDefault<double>("datum_depth", 2753.87) - zeroDepth;
double datum_pressure_barsa = param.getDefault<double>("datum_pressure", 248.22);
double datum_pressure = Opm::unit::convert::from(datum_pressure_barsa, Opm::unit::barsa);
double wo_contact_depth = param.getDefault<double>("wo_contact_depth", 3032.76) - zeroDepth;
double go_contact_depth = param.getDefault<double>("go_contact_depth", 2682.24) - zeroDepth;
double connate_water_saturation = param.getDefault<double>("connate_water_saturation", 0.151090);
double residual_oil_saturation = param.getDefault<double>("residual_oil_saturation", 0.118510);
double initial_mixture_gas = param.getDefault("initial_mixture_gas", 247.43);
double initial_mixture_oil = param.getDefault("initial_mixture_oil", 1.0);
// Initial fluid state
CompVec oil_sample(0.0);
oil_sample[Fluid::Oil] = initial_mixture_oil;
oil_sample[Fluid::Gas] = initial_mixture_gas;
CompVec water_sample(0.0);
water_sample[Fluid::Water] = 1.0;
simstate.cell_z_.resize(grid.numCells());
simstate.cell_pressure_.resize(grid.numCells());
// Datum -cell
// For now, assume that datum_depth corresponds the centroid of cell 0 (reasonable approx)
simstate.cell_pressure_[0] = datum_pressure;
typename Fluid::FluidState state = fluid.computeState(simstate.cell_pressure_[0],oil_sample);
simstate.cell_z_[0] = oil_sample;
simstate.cell_z_[0] *= (1.0-connate_water_saturation)/state.total_phase_volume_density_;
state = fluid.computeState(simstate.cell_pressure_[0],water_sample);
simstate.cell_z_[0][Fluid::Water] = water_sample[Fluid::Water];
simstate.cell_z_[0][Fluid::Water] *= connate_water_saturation/state.total_phase_volume_density_;
// Rest of the cells -- NOTE: Assume uniform cell properties in y-direction
for (int i=0; i<nx; ++i) {
int k0=i*nz;
for (int k=0; k<nz; ++k) {
int kk=k0+k;
if (i>0 && k==0) {
computeCellState(grid, fluid, gravity,
kk, kk-nz, wo_contact_depth, go_contact_depth, connate_water_saturation,
residual_oil_saturation, simstate);
} else if (k>0) {
computeCellState(grid, fluid, gravity,
kk, kk-1, wo_contact_depth, go_contact_depth, connate_water_saturation,
residual_oil_saturation, simstate);
}
// Copy cell properties to y-layers
for (int j=1; j<ny; ++j) {
int jj = j*nx*nz + kk;
simstate.cell_z_[jj] = simstate.cell_z_[kk];
simstate.cell_pressure_[jj] = simstate.cell_pressure_[kk];
}
}
}
}
bool computeCellState(const Grid& grid,
const Fluid& fluid,
typename Grid::Vector gravity,
int iCell,
int iRef,
double wo_contact_depth,
double /* go_contact_depth */,
double connate_water_saturation,
double residual_oil_saturation,
State& simstate)
{
typedef typename Fluid::PhaseVec PhaseVec;
const int maxCnt = 30;
const double eps = 1.0e-8;
simstate.cell_z_[iCell] = simstate.cell_z_[iRef];
bool below_wo_contact = false;
if (grid.cellCentroid(iCell)[2] > wo_contact_depth)
below_wo_contact = true;
double gZ = (grid.cellCentroid(iCell) - grid.cellCentroid(iRef))*gravity;
double fluid_vol_dens;
int cnt =0;
do {
double rho = 0.5*(simstate.cell_z_[iCell]*fluid.surfaceDensities()
+ simstate.cell_z_[iRef]*fluid.surfaceDensities());
double press = rho*gZ + simstate.cell_pressure_[iRef][0];
simstate.cell_pressure_[iCell] = PhaseVec(press);
typename Fluid::FluidState state = fluid.computeState(simstate.cell_pressure_[iCell], simstate.cell_z_[iCell]);
fluid_vol_dens = state.total_phase_volume_density_;
double oil_vol_dens = state.phase_volume_density_[Fluid::Liquid]
+ state.phase_volume_density_[Fluid::Vapour];
double wat_vol_dens = state.phase_volume_density_[Fluid::Aqua];
if (below_wo_contact) {
simstate.cell_z_[iCell][Fluid::Oil] *= residual_oil_saturation/oil_vol_dens;
simstate.cell_z_[iCell][Fluid::Gas] *= residual_oil_saturation/oil_vol_dens;
simstate.cell_z_[iCell][Fluid::Water] *= (1.0-residual_oil_saturation)/wat_vol_dens;
} else {
simstate.cell_z_[iCell][Fluid::Oil] *= (1.0-connate_water_saturation)/oil_vol_dens;
simstate.cell_z_[iCell][Fluid::Gas] *= (1.0-connate_water_saturation)/oil_vol_dens;
simstate.cell_z_[iCell][Fluid::Water] *= connate_water_saturation/wat_vol_dens;
}
++cnt;
} while (std::fabs(fluid_vol_dens-1.0) > eps && cnt < maxCnt);
if (cnt == maxCnt) {
std::cout << "z_cell_[" << iCell << "]: " << simstate.cell_z_[iCell]
<< " pressure: " << simstate.cell_pressure_[iCell][Fluid::Liquid]
<< " cnt: " << cnt
<< " eps: " << std::fabs(fluid_vol_dens-1.0) << std::endl;
}
return (cnt < maxCnt);
}
};
} // namespace Opm
#endif // OPM_BLACKOILINITIALIZATION_HEADER_INCLUDED

587
opm/porsol/blackoil/BlackoilSimulator.hpp
View File

@ -1,587 +0,0 @@
/*
Copyright 2011 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILSIMULATOR_HEADER_INCLUDED
#define OPM_BLACKOILSIMULATOR_HEADER_INCLUDED
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <boost/filesystem/convenience.hpp>
#include <boost/lexical_cast.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/porsol/common/BoundaryConditions.hpp>
#include <opm/porsol/blackoil/BlackoilInitialization.hpp>
#include <opm/porsol/common/SimulatorUtilities.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <iostream>
#include <string>
#include <vector>
#include <numeric>
namespace Opm
{
template<class GridT, class Rock, class FluidT, class Wells, class FlowSolver, class TransportSolver>
class BlackoilSimulator
{
public:
void init(const Opm::ParameterGroup& param);
void simulate();
typedef GridT Grid;
typedef FluidT Fluid;
typedef typename Fluid::CompVec CompVec;
typedef typename Fluid::PhaseVec PhaseVec;
struct State
{
std::vector<PhaseVec> cell_pressure_;
std::vector<PhaseVec> face_pressure_;
std::vector<double> well_bhp_pressure_;
std::vector<double> well_perf_pressure_;
std::vector<double> well_perf_flux_;
std::vector<CompVec> cell_z_;
};
private:
Grid grid_;
Rock rock_;
Fluid fluid_;
Wells wells_;
FlowSolver flow_solver_;
TransportSolver transport_solver_;
Opm::BasicBoundaryConditions<true, false> flow_bc_;
typename Grid::Vector gravity_;
std::vector<double> src_;
State state_;
PhaseVec bdy_pressure_;
CompVec bdy_z_;
double total_time_;
double initial_stepsize_;
bool increase_stepsize_;
double stepsize_increase_factor_;
double minimum_stepsize_;
double maximum_stepsize_;
std::vector<double> report_times_;
bool do_impes_;
bool ignore_impes_stability_;
std::string output_dir_;
int output_interval_;
static void output(const Grid& grid,
const Fluid& fluid,
const State& simstate,
const std::vector<double>& face_flux,
const int step,
const std::string& filebase);
};
// Method implementations below.
template<class Grid, class Rock, class Fluid, class Wells, class FlowSolver, class TransportSolver>
void
BlackoilSimulator<Grid, Rock, Fluid, Wells, FlowSolver, TransportSolver>::
init(const Opm::ParameterGroup& param)
{
using namespace Opm;
std::string fileformat = param.getDefault<std::string>("fileformat", "cartesian");
if (fileformat == "eclipse") {
Opm::ParseContext parseContext;
Opm::Parser parser;
auto deck = parser.parseFile(param.get<std::string>("filename") , parseContext);
Opm::EclipseGrid inputGrid( deck );
double z_tolerance = param.getDefault<double>("z_tolerance", 0.0);
bool periodic_extension = param.getDefault<bool>("periodic_extension", false);
bool turn_normals = param.getDefault<bool>("turn_normals", false);
grid_.processEclipseFormat(inputGrid, z_tolerance, periodic_extension, turn_normals);
double perm_threshold_md = param.getDefault("perm_threshold_md", 0.0);
double perm_threshold = Opm::unit::convert::from(perm_threshold_md, Opm::prefix::milli*Opm::unit::darcy);
rock_.init(deck, grid_.globalCell(), perm_threshold);
fluid_.init(deck);
wells_.init(deck, grid_, rock_);
} else if (fileformat == "cartesian") {
std::array<int, 3> dims = {{ param.getDefault<int>("nx", 1),
param.getDefault<int>("ny", 1),
param.getDefault<int>("nz", 1) }};
std::array<double, 3> cellsz = {{ param.getDefault<double>("dx", 1.0),
param.getDefault<double>("dy", 1.0),
param.getDefault<double>("dz", 1.0) }};
grid_.createCartesian(dims, cellsz);
double default_poro = param.getDefault("default_poro", 1.0);
double default_perm_md = param.getDefault("default_perm_md", 100.0);
double default_perm = unit::convert::from(default_perm_md, prefix::milli*unit::darcy);
OPM_MESSAGE("Warning: For generated cartesian grids, we use uniform rock properties.");
rock_.init(grid_.size(0), default_poro, default_perm);
Opm::ParseContext parseContext;
Opm::Parser parser;
auto deck = parser.parseFile(param.get<std::string>("filename") , parseContext);
fluid_.init(deck);
wells_.init(deck, grid_, rock_);
} else {
OPM_THROW(std::runtime_error, "Unknown file format string: " << fileformat);
}
flow_solver_.init(param);
transport_solver_.init(param);
if (param.has("timestep_file")) {
std::ifstream is(param.get<std::string>("timestep_file").c_str());
std::istream_iterator<double> beg(is);
std::istream_iterator<double> end;
report_times_.assign(beg, end);
// File contains deltas, we want accumulated times.
std::partial_sum(report_times_.begin(), report_times_.end(), report_times_.begin());
assert(!report_times_.empty());
total_time_ = report_times_.back();
initial_stepsize_ = report_times_.front();
} else {
total_time_ = param.getDefault("total_time", 30*unit::day);
initial_stepsize_ = param.getDefault("initial_stepsize", 1.0*unit::day);
increase_stepsize_ = param.getDefault("increase_stepsize", false);
if (increase_stepsize_) {
stepsize_increase_factor_ = param.getDefault("stepsize_increase_factor", 1.5);
maximum_stepsize_ = param.getDefault("maximum_stepsize", 1.0*unit::day);
} else {
stepsize_increase_factor_ = 1.0;
maximum_stepsize_ = 1e100;
}
}
minimum_stepsize_ = param.getDefault("minimum_stepsize", 0.0);
do_impes_ = param.getDefault("do_impes", false);
if (do_impes_) {
ignore_impes_stability_ = param.getDefault("ignore_impes_stability", false);
}
output_dir_ = param.getDefault<std::string>("output_dir", "output");
output_interval_ = param.getDefault("output_interval", 1);
// Boundary conditions.
typedef Opm::FlowBC BC;
flow_bc_.resize(7);
bool bdy_dirichlet = param.getDefault("bdy_dirichlet", false);
if (bdy_dirichlet) {
flow_bc_.flowCond(1) = BC(BC::Dirichlet, param.get<double>("bdy_pressure_left"));
flow_bc_.flowCond(2) = BC(BC::Dirichlet, param.get<double>("bdy_pressure_right"));
} else if (param.getDefault("lateral_dirichlet", false)) {
flow_bc_.flowCond(1) = BC(BC::Dirichlet, -17.0); // Use a negative value to instruct flow solver
flow_bc_.flowCond(2) = BC(BC::Dirichlet, -17.0); // to use initial face pressures (hydrostatic)
flow_bc_.flowCond(3) = BC(BC::Dirichlet, -17.0); // as boundary conditions.
flow_bc_.flowCond(4) = BC(BC::Dirichlet, -17.0);
}
// Gravity.
gravity_ = 0.0;
if (param.has("gravity")) {
std::string g = param.get<std::string>("gravity");
if (g == "standard") {
gravity_[2] = Opm::unit::gravity;
} else {
gravity_[2] = boost::lexical_cast<double>(g);
}
}
// Initial state.
if (param.getDefault("spe9_init", false)) {
SPE9Initialization<BlackoilSimulator> initializer;
initializer.init(param, grid_, fluid_, gravity_, state_);
} else {
BasicInitialization<BlackoilSimulator> initializer;
initializer.init(param, grid_, fluid_, gravity_, state_);
}
// Write initial state to std::cout
/*
for (int cell = 0; cell < grid_.numCells(); ++cell) {
std::cout.precision(2);
std::cout << std::fixed << std::showpoint;
std::cout << std::setw(5) << cell << std::setw(12) << grid_.cellCentroid(cell)[0]
<< std::setw(12) << grid_.cellCentroid(cell)[1]
<< std::setw(12) << grid_.cellCentroid(cell)[2]
<< std::setw(20) << state_.cell_pressure_[cell][0]
<< std::setw(15) << state_.cell_z_[cell][0]
<< std::setw(15) << state_.cell_z_[cell][1]
<< std::setw(15) << state_.cell_z_[cell][2]
<< std::endl;
if ((cell+1)%nz == 0) {
std::cout << "------------------------------------------------------------------------------------------------------------------" << std::endl;
}
}
*/
bdy_z_ = flow_solver_.inflowMixture();
bdy_pressure_ = 300.0*Opm::unit::barsa;
// PhaseVec bdy_pressure_(100.0*Opm::unit::barsa); // WELLS
// Rescale z values so that pore volume is filled exactly
// (to get zero initial volume discrepancy).
for (int cell = 0; cell < grid_.numCells(); ++cell) {
typename Fluid::FluidState state = fluid_.computeState(state_.cell_pressure_[cell], state_.cell_z_[cell]);
double fluid_vol_dens = state.total_phase_volume_density_;
state_.cell_z_[cell] *= 1.0/fluid_vol_dens;
}
int num_faces = grid_.numFaces();
state_.face_pressure_.resize(num_faces);
for (int face = 0; face < num_faces; ++face) {
int bid = grid_.boundaryId(face);
if (flow_bc_.flowCond(bid).isDirichlet() && flow_bc_.flowCond(bid).pressure() >= 0.0) {
state_.face_pressure_[face] = flow_bc_.flowCond(bid).pressure();
} else {
int c[2] = { grid_.faceCell(face, 0), grid_.faceCell(face, 1) };
state_.face_pressure_[face] = 0.0;
int num = 0;
for (int j = 0; j < 2; ++j) {
if (c[j] >= 0) {
state_.face_pressure_[face] += state_.cell_pressure_[c[j]];
++num;
}
}
state_.face_pressure_[face] /= double(num);
}
}
// Flow solver setup.
flow_solver_.setup(grid_, rock_, fluid_, wells_, gravity_, flow_bc_, &(state_.face_pressure_));
// Transport solver setup.
transport_solver_.setup(grid_, rock_, fluid_, wells_, flow_solver_.faceTransmissibilities(), gravity_);
// Simple source terms.
src_.resize(grid_.numCells(), 0.0);
// Set initial well perforation pressures equal to cell pressures,
// and perforation fluxes equal to zero.
// Set initial well bhp values to the target if bhp well, or to
// first perforation pressure if not.
state_.well_perf_pressure_.clear();
for (int well = 0; well < wells_.numWells(); ++well) {
int num_perf = wells_.numPerforations(well);
for (int perf = 0; perf < num_perf; ++perf) {
int cell = wells_.wellCell(well, perf);
state_.well_perf_pressure_.push_back(state_.cell_pressure_[cell][Fluid::Liquid]);
}
if (wells_.control(well) == Wells::Pressure) {
state_.well_bhp_pressure_.push_back(wells_.target(well));
} else {
int cell = wells_.wellCell(well, 0);
state_.well_bhp_pressure_.push_back(state_.cell_pressure_[cell][Fluid::Liquid]);
}
}
state_.well_perf_flux_.clear();
state_.well_perf_flux_.resize(state_.well_perf_pressure_.size(), 0.0);
wells_.update(grid_.numCells(), state_.well_perf_pressure_, state_.well_perf_flux_);
// Check for unused parameters (potential typos).
if (param.anyUnused()) {
std::cout << "***** WARNING: Unused parameters: *****\n";
param.displayUsage();
}
// Write parameters used to file, ensuring directory exists.
std::string paramfilename = output_dir_ + "/simulator-parameters.param";
boost::filesystem::path fpath(paramfilename);
if (fpath.has_branch_path()) {
create_directories(fpath.branch_path());
}
param.writeParam(paramfilename);
}
template<class Grid, class Rock, class Fluid, class Wells, class FlowSolver, class TransportSolver>
void
BlackoilSimulator<Grid, Rock, Fluid, Wells, FlowSolver, TransportSolver>::
simulate()
{
double voldisclimit = flow_solver_.volumeDiscrepancyLimit();
double stepsize = initial_stepsize_;
double current_time = 0.0;
int step = 0;
std::vector<double> face_flux;
State start_state;
std::string output_name = output_dir_ + "/" + "blackoil-output";
while (current_time < total_time_) {
start_state = state_;
// Do not run past total_time_.
if (current_time + stepsize > total_time_) {
stepsize = total_time_ - current_time;
}
std::cout << "\n\n================ Simulation step number " << step
<< " ==============="
<< "\n Current time (days) " << Opm::unit::convert::to(current_time, Opm::unit::day)
<< "\n Current stepsize (days) " << Opm::unit::convert::to(stepsize, Opm::unit::day)
<< "\n Total time (days) " << Opm::unit::convert::to(total_time_, Opm::unit::day)
<< "\n" << std::endl;
// Solve flow system.
enum FlowSolver::ReturnCode result
= flow_solver_.solve(state_.cell_pressure_, state_.face_pressure_, state_.cell_z_, face_flux,
state_.well_bhp_pressure_,
state_.well_perf_pressure_, state_.well_perf_flux_, src_, stepsize);
// Check if the flow solver succeeded.
if (result == FlowSolver::VolumeDiscrepancyTooLarge) {
OPM_THROW(std::runtime_error, "Flow solver refused to run due to too large volume discrepancy.");
} else if (result == FlowSolver::FailedToConverge) {
std::cout << "********* Nonlinear convergence failure: Shortening (pressure) stepsize, redoing step number " << step <<" **********" << std::endl;
stepsize *= 0.5;
state_ = start_state;
wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
continue;
}
assert(result == FlowSolver::SolveOk);
// Update wells with new perforation pressures and fluxes.
wells_.update(grid_.numCells(), state_.well_perf_pressure_, state_.well_perf_flux_);
// Transport and check volume discrepancy.
bool voldisc_ok = true;
if (!do_impes_) {
double actual_computed_time
= transport_solver_.transport(bdy_pressure_, bdy_z_,
face_flux, state_.cell_pressure_, state_.face_pressure_,
stepsize, voldisclimit, state_.cell_z_);
voldisc_ok = (actual_computed_time == stepsize);
if (voldisc_ok) {
// Just for output.
flow_solver_.volumeDiscrepancyAcceptable(state_.cell_pressure_, state_.face_pressure_,
state_.well_perf_pressure_, state_.cell_z_, stepsize);
}
} else {
// First check IMPES stepsize.
double max_dt = ignore_impes_stability_ ? 1e100 : flow_solver_.stableStepIMPES();
if (ignore_impes_stability_) {
std::cout << "Timestep was " << stepsize << " and max stepsize was not computed." << std::endl;
} else {
std::cout << "Timestep was " << stepsize << " and max stepsize was " << max_dt << std::endl;
}
if (stepsize < max_dt || stepsize <= minimum_stepsize_) {
flow_solver_.doStepIMPES(state_.cell_z_, stepsize);
voldisc_ok = flow_solver_.volumeDiscrepancyAcceptable(state_.cell_pressure_, state_.face_pressure_,
state_.well_perf_pressure_, state_.cell_z_, stepsize);
} else {
// Restarting step.
stepsize = max_dt/1.5;
std::cout << "Restarting pressure step with new timestep " << stepsize << std::endl;
state_ = start_state;
wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
continue;
}
}
// If discrepancy too large, redo entire pressure step.
if (!voldisc_ok) {
std::cout << "********* Too large volume discrepancy: Shortening (pressure) stepsize, redoing step number " << step <<" **********" << std::endl;
stepsize *= 0.5;
state_ = start_state;
wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
continue;
}
// Adjust time.
current_time += stepsize;
if (voldisc_ok && increase_stepsize_ && stepsize < maximum_stepsize_) {
stepsize *= stepsize_increase_factor_;
stepsize = std::min(maximum_stepsize_, stepsize);
}
// If using given timesteps, set stepsize to match.
if (!report_times_.empty()) {
if (current_time >= report_times_[step]) {
bool output_now = ((step + 1) % output_interval_ == 0);
if (output_now) {
output(grid_, fluid_, state_, face_flux, step, output_name);
}
++step;
if (step == int(report_times_.size())) {
break;
}
}
stepsize = report_times_[step] - current_time;
} else {
bool output_now = ((step + 1) % output_interval_ == 0);
if (output_now) {
output(grid_, fluid_, state_, face_flux, step, output_name);
}
++step;
}
}
if (step % output_interval_ != 0) {
// Output was not written at last step, write final output.
output(grid_, fluid_, state_, face_flux, step - 1, output_name);
}
}
template<class Grid, class Rock, class Fluid, class Wells, class FlowSolver, class TransportSolver>
void
BlackoilSimulator<Grid, Rock, Fluid, Wells, FlowSolver, TransportSolver>::
output(const Grid& grid,
const Fluid& fluid,
const State& simstate,
const std::vector<double>& face_flux,
const int step,
const std::string& filebase)
{
// Compute saturations, total fluid volume density and mass fractions.
int num_cells = grid.numCells();
std::vector<typename Fluid::PhaseVec> sat(num_cells);
std::vector<typename Fluid::PhaseVec> mass_frac(num_cells);
std::vector<double> totflvol_dens(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
typename Fluid::FluidState fstate = fluid.computeState(simstate.cell_pressure_[cell], simstate.cell_z_[cell]);
sat[cell] = fstate.saturation_;
totflvol_dens[cell] = fstate.total_phase_volume_density_;
double totMass_dens = simstate.cell_z_[cell]*fluid.surfaceDensities();
mass_frac[cell][Fluid::Water] = simstate.cell_z_[cell][Fluid::Water]*fluid.surfaceDensities()[Fluid::Water]/totMass_dens;
mass_frac[cell][Fluid::Oil] = simstate.cell_z_[cell][Fluid::Oil]*fluid.surfaceDensities()[Fluid::Oil]/totMass_dens;
mass_frac[cell][Fluid::Gas] = simstate.cell_z_[cell][Fluid::Gas]*fluid.surfaceDensities()[Fluid::Gas]/totMass_dens;
}
// Ensure directory exists.
boost::filesystem::path fpath(filebase);
if (fpath.has_branch_path()) {
create_directories(fpath.branch_path());
}
// Output to VTK.
std::vector<typename Grid::Vector> cell_velocity;
Opm::estimateCellVelocitySimpleInterface(cell_velocity, grid, face_flux);
// Dune's vtk writer wants multi-component data to be flattened.
std::vector<double> cell_pressure_flat(&*simstate.cell_pressure_.front().begin(),
&*simstate.cell_pressure_.back().end());
std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
&*cell_velocity.back().end());
std::vector<double> z_flat(&*simstate.cell_z_.front().begin(),
&*simstate.cell_z_.back().end());
std::vector<double> sat_flat(&*sat.front().begin(),
&*sat.back().end());
std::vector<double> mass_frac_flat(&*mass_frac.front().begin(),
&*mass_frac.back().end());
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafGridView());
#else
Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafView());
#endif
vtkwriter.addCellData(cell_pressure_flat, "pressure", Fluid::numPhases);
vtkwriter.addCellData(cell_velocity_flat, "velocity", Grid::dimension);
vtkwriter.addCellData(z_flat, "z", Fluid::numComponents);
vtkwriter.addCellData(sat_flat, "sat", Fluid::numPhases);
vtkwriter.addCellData(mass_frac_flat, "massFrac", Fluid::numComponents);
vtkwriter.addCellData(totflvol_dens, "total fl. vol.");
vtkwriter.write(filebase + '-' + boost::lexical_cast<std::string>(step),
Dune::VTK::ascii);
// Dump data for Matlab.
std::vector<double> zv[Fluid::numComponents];
for (int comp = 0; comp < Fluid::numComponents; ++comp) {
zv[comp].resize(grid.numCells());
for (int cell = 0; cell < grid.numCells(); ++cell) {
zv[comp][cell] = simstate.cell_z_[cell][comp];
}
}
std::vector<double> sv[Fluid::numPhases];
for (int phase = 0; phase < Fluid::numPhases; ++phase) {
sv[phase].resize(grid.numCells());
for (int cell = 0; cell < grid.numCells(); ++cell) {
sv[phase][cell] = sat[cell][phase];
}
}
std::string matlabdumpname(filebase + "-");
matlabdumpname += boost::lexical_cast<std::string>(step);
matlabdumpname += ".dat";
std::ofstream dump(matlabdumpname.c_str());
dump.precision(15);
std::vector<double> liq_press(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
liq_press[cell] = simstate.cell_pressure_[cell][Fluid::Liquid];
}
// Liquid phase pressure.
std::copy(liq_press.begin(), liq_press.end(),
std::ostream_iterator<double>(dump, " "));
dump << '\n';
// z (3 components)
for (int comp = 0; comp < Fluid::numComponents; ++comp) {
std::copy(zv[comp].begin(), zv[comp].end(),
std::ostream_iterator<double>(dump, " "));
dump << '\n';
}
// s (3 components)
for (int phase = 0; phase < Fluid::numPhases; ++phase) {
std::copy(sv[phase].begin(), sv[phase].end(),
std::ostream_iterator<double>(dump, " "));
dump << '\n';
}
// Total fluid volume
std::copy(totflvol_dens.begin(), totflvol_dens.end(),
std::ostream_iterator<double>(dump, " "));
dump << '\n';
// Well report ...
const double seconds_pr_day = 3600.*24.;
for (unsigned int perf=0; perf<Wells::WellReport::report()->perfPressure.size(); ++perf) {
dump << std::setw(8) << Wells::WellReport::report()->cellId[perf] << " "
<< std::setw(22) << Wells::WellReport::report()->perfPressure[perf] << " "
<< std::setw(22) << Wells::WellReport::report()->cellPressure[perf] << " "
<< std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Water] << " "
<< std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Oil] << " "
<< std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Gas] << '\n';
}
dump << '\n';
}
} // namespace Opm
#endif // OPM_BLACKOILSIMULATOR_HEADER_INCLUDED

619
opm/porsol/blackoil/BlackoilWells.hpp
View File

@ -1,619 +0,0 @@
/*
Copyright 2011 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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILWELLS_HEADER_INCLUDED
#define OPM_BLACKOILWELLS_HEADER_INCLUDED
#include <opm/porsol/blackoil/fluid/BlackoilDefs.hpp>
#include <dune/grid/CpGrid.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/SparseTable.hpp>
#include <opm/porsol/common/Rock.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <dune/common/fvector.hh>
#include <vector>
#include <iostream>
namespace Opm
{
/// A class designed to encapsulate a set of rate- or
/// pressure-controlled wells in the black-oil setting.
class BlackoilWells : public BlackoilDefs
{
public:
void init(const Opm::Deck& deck,
const Dune::CpGrid& grid,
const Opm::Rock<3>& rock);
// Well-centric interface. Mostly used by pressure solver.
int numWells() const;
enum WellType { Injector, Producer };
WellType type(int wellnum) const;
enum WellControl { Rate, Pressure };
WellControl control(int wellnum) const;
double target(int wellnum) const;
double referenceDepth(int wellnum) const;
int numPerforations(int wellnum) const;
int wellCell(int wellnum, int perfnum) const;
double wellIndex(int wellnum, int perfnum) const;
// Updating rates and pressures after pressure solve.
void update(int num_cells,
const std::vector<double>& well_perf_pressures,
const std::vector<double>& well_perf_fluxes);
// Cell-centric interface. Mostly used by transport solver.
double perforationPressure(int cell) const;
double wellToReservoirFlux(int cell) const;
CompVec injectionMixture(int cell) const;
// Simple well report ...
class WellReport
{
public:
static WellReport* report()
{
if (well_report_ == 0x0)
well_report_ = new WellReport;
return well_report_;
}
void clearAll()
{
perfPressure.clear();
cellPressure.clear();
massRate.clear();
cellId.clear();
}
std::vector<double> perfPressure;
std::vector<double> cellPressure;
std::vector<BlackoilDefs::CompVec> massRate; // (surface volumes)
std::vector<int> cellId;
protected:
WellReport() {}
private:
static WellReport* well_report_;
};
private:
// Use the Peaceman well model to compute well indices
double computeWellIndex(double radius, const Dune::FieldVector<double, 3>& cubical,
const Opm::Rock<3>::PermTensor& permeability, double skin_factor) const;
struct WellData { WellType type; WellControl control; double target; double reference_bhp_depth; };
std::vector<WellData> well_data_;
struct PerfData { int cell; double well_index; };
Opm::SparseTable<PerfData> perf_data_;
std::vector<double> well_cell_flux_;
std::vector<double> well_cell_pressure_;
Dune::FieldVector<double, 3> injection_mixture_;
std::vector<std::string> well_names_;
};
BlackoilWells::WellReport* BlackoilWells::WellReport::well_report_ = 0x0;
// ------------ Method implementations --------------
// Forward declaration of some helper functions.
namespace
{
int prod_control_mode(const std::string& control);
int inje_control_mode(const std::string& control);
template<class grid_t>
const Dune::FieldVector<double,3> getCubeDim(const grid_t& grid, int cell);
} // anon namespace
inline void BlackoilWells::init(const Opm::Deck& deck,
const Dune::CpGrid& grid,
const Opm::Rock<3>& rock)
{
if (!deck.hasKeyword("WELSPECS")) {
OPM_MESSAGE("Missing keyword \"WELSPECS\" in deck, initializing no wells.");
return;
}
if (!deck.hasKeyword("COMPDAT")) {
OPM_MESSAGE("Missing keyword \"COMPDAT\" in deck, initializing no wells.");
return;
}
if (!(deck.hasKeyword("WCONINJE") || deck.hasKeyword("WCONPROD")) ) {
OPM_THROW(std::runtime_error, "Needed field is missing in file");
}
using namespace Opm;
// Get WELLSPECS data
const auto& welspecsKeyword = deck.getKeyword("WELSPECS");
const int num_welspecs = welspecsKeyword.size();
well_names_.reserve(num_welspecs);
well_data_.reserve(num_welspecs);
for (int i=0; i<num_welspecs; ++i) {
well_names_.push_back(welspecsKeyword.getRecord(i).getItem("WELL").get< std::string >(0));
WellData wd;
well_data_.push_back(wd);
well_data_.back().reference_bhp_depth =
welspecsKeyword.getRecord(i).getItem("REF_DEPTH").getSIDouble(0);
}
// Get COMPDAT data
const auto& compdatKeyword = deck.getKeyword("COMPDAT");
const int num_compdats = compdatKeyword.size();
std::vector<std::vector<PerfData> > wellperf_data(num_welspecs);
for (int kw=0; kw<num_compdats; ++kw) {
const auto& compdatRecord = compdatKeyword.getRecord(kw);
std::string name = compdatRecord.getItem("WELL").get< std::string >(0);
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
// global_cell is a map from compressed cells to Cartesian grid cells
const std::vector<int>& global_cell = grid.globalCell();
const std::array<int, 3>& cpgdim = grid.logicalCartesianSize();
std::map<int,int> cartesian_to_compressed;
for (int i=0; i<int(global_cell.size()); ++i) {
cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
}
bool found = false;
for (int wix=0; wix<num_welspecs; ++wix) {
if (well_names_[wix].compare(0,len, name) == 0) { //equal
int ix = compdatRecord.getItem("I").get< int >(0) - 1;
int jy = compdatRecord.getItem("J").get< int >(0) - 1;
int kz1 = compdatRecord.getItem("K1").get< int >(0) - 1;
int kz2 = compdatRecord.getItem("K2").get< int >(0) - 1;
for (int kz = kz1; kz <= kz2; ++kz) {
int cart_grid_indx = ix + cpgdim[0]*(jy + cpgdim[1]*kz);
std::map<int, int>::const_iterator cgit =
cartesian_to_compressed.find(cart_grid_indx);
if (cgit == cartesian_to_compressed.end()) {
OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << ix << ' ' << jy << ' '
<< kz << " not found!");
}
int cell = cgit->second;
PerfData pd;
pd.cell = cell;
if (compdatRecord.getItem("CF").getSIDouble(0) > 0.0) {
pd.well_index = compdatRecord.getItem("CF").getSIDouble(0);
} else {
double radius = 0.5*compdatRecord.getItem("DIAMETER").getSIDouble(0);
if (radius <= 0.0) {
radius = 0.5*unit::feet;
OPM_MESSAGE("Warning: Well bore internal radius set to " << radius);
}
Dune::FieldVector<double, 3> cubical = getCubeDim(grid, cell);
const Rock<3>::PermTensor permeability = rock.permeability(cell);
pd.well_index = computeWellIndex(radius, cubical, permeability,
compdatRecord.getItem("SKIN").getSIDouble(0));
}
wellperf_data[wix].push_back(pd);
}
found = true;
break;
}
}
if (!found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << compdatRecord.getItem("WELL").get< std::string >(0)
<< " in COMPDAT");
}
}
for (int w = 0; w < num_welspecs; ++w) {
perf_data_.appendRow(wellperf_data[w].begin(), wellperf_data[w].end());
if (well_data_[w].reference_bhp_depth == -1e100) {
// It was defaulted. Set reference depth to minimum perforation depth.
double min_depth = 1e100;
int num_wperfs = wellperf_data[w].size();
for (int perf = 0; perf < num_wperfs; ++perf) {
double depth = grid.cellCentroid(wellperf_data[w][perf].cell)[2];
min_depth = std::min(min_depth, depth);
}
well_data_[w].reference_bhp_depth = min_depth;
}
}
// Get WCONINJE data
injection_mixture_ = 0.0;
if (deck.hasKeyword("WCONINJE")) {
const auto& wconinjeKeyword = deck.getKeyword("WCONINJE");
const int num_wconinjes = wconinjeKeyword.size();
int injector_component = -1;
for (int kw=0; kw<num_wconinjes; ++kw) {
const auto& wconinjeRecord = wconinjeKeyword.getRecord(kw);
std::string name = wconinjeRecord.getItem("WELL").get< std::string >(0);
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix=0; wix<num_welspecs; ++wix) {
if (well_names_[wix].compare(0,len, name) == 0) { //equal
well_found = true;
well_data_[wix].type = Injector;
int m = inje_control_mode(wconinjeRecord.getItem("CMODE").get< std::string >(0));
switch(m) {
case 0: // RATE
well_data_[wix].control = Rate;
// TODO: convert rate to SI!
well_data_[wix].target = wconinjeRecord.getItem("RATE").get< double >(0);
break;
case 1: // RESV
well_data_[wix].control = Rate;
// TODO: convert rate to SI!
well_data_[wix].target = wconinjeRecord.getItem("RESV").get< double >(0);
break;
case 2: // BHP
well_data_[wix].control = Pressure;
well_data_[wix].target = wconinjeRecord.getItem("BHP").getSIDouble(0);
break;
case 3: // THP
well_data_[wix].control = Pressure;
well_data_[wix].target = wconinjeRecord.getItem("THP").getSIDouble(0);
break;
default:
OPM_THROW(std::runtime_error, "Unknown well control mode; WCONIJE = "
<< wconinjeRecord.getItem("CMODE").get< std::string >(0)
<< " in input file");
}
int itp = -1;
if (wconinjeRecord.getItem("TYPE").get< std::string >(0) == "WATER") {
itp = Water;
} else if (wconinjeRecord.getItem("TYPE").get< std::string >(0) == "OIL") {
itp = Oil;
} else if (wconinjeRecord.getItem("TYPE").get< std::string >(0) == "GAS") {
itp = Gas;
}
if (itp == -1 || (injector_component != -1 && itp != injector_component)) {
if (itp == -1) {
OPM_THROW(std::runtime_error, "Error in injector specification, found no known fluid type.");
} else {
OPM_THROW(std::runtime_error, "Error in injector specification, we can only handle a single injection fluid.");
}
} else {
injector_component = itp;
}
}
}
if (!well_found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << wconinjeRecord.getItem("WELL").get< std::string >(0)
<< " in WCONINJE");
}
}
if (injector_component != -1) {
injection_mixture_[injector_component] = 1.0;
}
} else {
// No WCONINJE.
// This default is only invoked if production wells
// start injecting, (and only if there are no injection wells).
// In other words, only if we have a primary production type scenario.
// In this case we expect the flow to be very small, and not affect
// the simulation in any significant way.
injection_mixture_[Oil] = 1.0;
}
// Get WCONPROD data
if (deck.hasKeyword("WCONPROD")) {
const auto& wconprodKeyword = deck.getKeyword("WCONPROD");
const int num_wconprods = wconprodKeyword.size();
for (int kw=0; kw<num_wconprods; ++kw) {
const auto& wconprodRecord = wconprodKeyword.getRecord(kw);
std::string name = wconprodRecord.getItem("WELL").get< std::string >(0);
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix=0; wix<num_welspecs; ++wix) {
if (well_names_[wix].compare(0,len, name) == 0) { //equal
well_found = true;
well_data_[wix].type = Producer;
int m = prod_control_mode(wconprodRecord.getItem("CMODE").get< std::string >(0));
switch(m) {
case 0: // ORAT
well_data_[wix].control = Rate;
well_data_[wix].target = wconprodRecord.getItem("ORAT").getSIDouble(0);
break;
case 1: // WRAT
well_data_[wix].control = Rate;
well_data_[wix].target = wconprodRecord.getItem("WRAT").getSIDouble(0);
break;
case 2: // GRAT
well_data_[wix].control = Rate;
well_data_[wix].target = wconprodRecord.getItem("GRAT").getSIDouble(0);
break;
case 3: // LRAT
well_data_[wix].control = Rate;
well_data_[wix].target = wconprodRecord.getItem("LRAT").getSIDouble(0);
break;
case 4: // RESV
well_data_[wix].control = Rate;
well_data_[wix].target = wconprodRecord.getItem("RESV").getSIDouble(0);
break;
case 5: // BHP
well_data_[wix].control = Pressure;
well_data_[wix].target = wconprodRecord.getItem("BHP").getSIDouble(0);
break;
case 6: // THP
well_data_[wix].control = Pressure;
well_data_[wix].target = wconprodRecord.getItem("THP").getSIDouble(0);
break;
default:
OPM_THROW(std::runtime_error, "Unknown well control mode; WCONPROD = "
<< wconprodRecord.getItem("CMODE").get< std::string >(0)
<< " in input file");
}
}
}
if (!well_found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << wconprodRecord.getItem("WELL").get< std::string >(0)
<< " in WCONPROD");
}
}
}
// Get WELTARG data
if (deck.hasKeyword("WELTARG")) {
const auto& weltargKeyword = deck.getKeyword("WELTARG");
const int num_weltargs = weltargKeyword.size();
for (int kw=0; kw<num_weltargs; ++kw) {
const auto& weltargRecord = weltargKeyword.getRecord(kw);
std::string name = weltargRecord.getItem("WELL").get< std::string >(0);
std::string::size_type len = name.find('*');
if (len != std::string::npos) {
name = name.substr(0, len);
}
bool well_found = false;
for (int wix=0; wix<num_welspecs; ++wix) {
if (well_names_[wix].compare(0,len, name) == 0) { //equal
well_found = true;
// TODO: convert to SI!
well_data_[wix].target = weltargRecord.getItem("NEW_VALUE").get< double >(0);
break;
}
}
if (!well_found) {
OPM_THROW(std::runtime_error, "Undefined well name: " << weltargRecord.getItem("WELL").get< std::string >(0)
<< " in WELTARG");
}
}
}
// Debug output.
std::cout << "\t WELL DATA" << std::endl;
for(int i=0; i< int(well_data_.size()); ++i) {
std::cout << i << ": " << well_data_[i].type << " "
<< well_data_[i].control << " " << well_data_[i].target
<< std::endl;
}
std::cout << "\n\t PERF DATA" << std::endl;
for(int i=0; i< int(perf_data_.size()); ++i) {
for(int j=0; j< int(perf_data_[i].size()); ++j) {
std::cout << i << ": " << perf_data_[i][j].cell << " "
<< perf_data_[i][j].well_index << std::endl;
}
}
// Ensuring that they have the right size.
well_cell_pressure_.resize(grid.numCells(), -1e100);
well_cell_flux_.resize(grid.numCells(), 0.0);
}
inline int BlackoilWells::numWells() const
{
return well_data_.size();
}
inline BlackoilWells::WellType BlackoilWells::type(int wellnum) const
{
return well_data_[wellnum].type;
}
inline BlackoilWells::WellControl BlackoilWells::control(int wellnum) const
{
return well_data_[wellnum].control;
}
inline double BlackoilWells::target(int wellnum) const
{
return well_data_[wellnum].target;
}
inline double BlackoilWells::referenceDepth(int wellnum) const
{
return well_data_[wellnum].reference_bhp_depth;
}
inline int BlackoilWells::numPerforations(int wellnum) const
{
return perf_data_[wellnum].size();
}
inline int BlackoilWells::wellCell(int wellnum, int perfnum) const
{
return perf_data_[wellnum][perfnum].cell;
}
inline double BlackoilWells::wellIndex(int wellnum, int perfnum) const
{
return perf_data_[wellnum][perfnum].well_index;
}
inline void BlackoilWells::update(int num_cells,
const std::vector<double>& well_perf_pressures,
const std::vector<double>& well_perf_fluxes)
{
// Input is per perforation, data members store for all cells.
assert(perf_data_.dataSize() == int(well_perf_pressures.size()));
well_cell_pressure_.resize(num_cells, -1e100);
well_cell_flux_.resize(num_cells, 0.0);
int pcount = 0;
for (int w = 0; w < numWells(); ++w) {
for (int perf = 0; perf < numPerforations(w); ++perf) {
int cell = wellCell(w, perf);
well_cell_pressure_[cell] = well_perf_pressures[pcount];
well_cell_flux_[cell] = well_perf_fluxes[pcount];
++pcount;
}
}
assert(pcount == perf_data_.dataSize());
}
inline double BlackoilWells::wellToReservoirFlux(int cell) const
{
return well_cell_flux_[cell];
}
inline double BlackoilWells::perforationPressure(int cell) const
{
return well_cell_pressure_[cell];
}
inline Dune::FieldVector<double, 3> BlackoilWells::injectionMixture(int /* cell */) const
{
return injection_mixture_;
}
inline double BlackoilWells::computeWellIndex(double radius,
const Dune::FieldVector<double, 3>& cubical,
const Opm::Rock<3>::PermTensor& permeability,
double skin_factor) const
{
// Use the Peaceman well model to compute well indices.
// radius is the radius of the well.
// cubical contains [dx, dy, dz] of the cell.
// (Note that the well model asumes that each cell is a cuboid).
// permeability is the permeability of the given cell.
// returns the well index of the cell.
// sse: Using the Peaceman modell.
// NOTE: The formula is valid for cartesian grids, so the result can be a bit
// (in worst case: there is no upper bound for the error) off the mark.
double effective_perm = sqrt(permeability(0,0) * permeability(1,1));
// sse: The formula for r_0 can be found on page 39 of
// "Well Models for Mimetic Finite Differerence Methods and Improved Representation
// of Wells in Multiscale Methods" by Ingeborg Skjelkvåle Ligaarden.
assert(permeability(0,0) > 0.0);
assert(permeability(1,1) > 0.0);
double kxoy = permeability(0,0) / permeability(1,1);
double kyox = permeability(1,1) / permeability(0,0);
double r0_denominator = pow(kyox, 0.25) + pow(kxoy, 0.25);
double r0_numerator = sqrt((sqrt(kyox)*cubical[0]*cubical[0]) +
(sqrt(kxoy)*cubical[1]*cubical[1]));
assert(r0_denominator > 0.0);
double r0 = 0.28 * r0_numerator / r0_denominator;
assert(radius > 0.0);
assert(r0 > 0.0);
if (r0 < radius) {
std::cout << "ERROR: Too big well radius detected.";
std::cout << "Specified well radius is " << radius
<< " while r0 is " << r0 << ".\n";
}
const double two_pi = 6.2831853071795864769252867665590057683943387987502116419498;
double wi_denominator = log(r0 / radius) + skin_factor;
double wi_numerator = two_pi * cubical[2];
assert(wi_denominator > 0.0);
double wi = effective_perm * wi_numerator / wi_denominator;
assert(wi > 0.0);
return wi;
}
// ------------- Helper functions for init() --------------
namespace
{
int prod_control_mode(const std::string& control){
const int num_prod_control_modes = 8;
static std::string prod_control_modes[num_prod_control_modes] =
{std::string("ORAT"), std::string("WRAT"), std::string("GRAT"),
std::string("LRAT"), std::string("RESV"), std::string("BHP"),
std::string("THP"), std::string("GRUP") };
int m = -1;
for (int i=0; i<num_prod_control_modes; ++i) {
if (control == prod_control_modes[i]) {
m = i;
break;
}
}
if (m >= 0) {
return m;
} else {
OPM_THROW(std::runtime_error, "Unknown well control mode = " << control << " in input file");
}
}
int inje_control_mode(const std::string& control)
{
const int num_inje_control_modes = 5;
static std::string inje_control_modes[num_inje_control_modes] =
{std::string("RATE"), std::string("RESV"), std::string("BHP"),
std::string("THP"), std::string("GRUP") };
int m = -1;
for (int i=0; i<num_inje_control_modes; ++i) {
if (control == inje_control_modes[i]) {
m = i;
break;
}
}
if (m >= 0) {
return m;
} else {
OPM_THROW(std::runtime_error, "Unknown well control mode = " << control << " in input file");
}
}
template<class grid_t>
const Dune::FieldVector<double,3> getCubeDim(const grid_t& grid, int cell)
{
Dune::FieldVector<double, 3> cube;
int num_local_faces = grid.numCellFaces(cell);
std::vector<double> x(num_local_faces);
std::vector<double> y(num_local_faces);
std::vector<double> z(num_local_faces);
for (int lf=0; lf<num_local_faces; ++ lf) {
int face = grid.cellFace(cell,lf);
const Dune::FieldVector<double,3>& centroid =
grid.faceCentroid(face);
x[lf] = centroid[0];
y[lf] = centroid[1];
z[lf] = centroid[2];
}
cube[0] = *max_element(x.begin(), x.end()) -
*min_element(x.begin(), x.end());
cube[1] = *max_element(y.begin(), y.end()) -
*min_element(y.begin(), y.end());
cube[2] = *max_element(z.begin(), z.end()) -
*min_element(z.begin(), z.end());
return cube;
}
} // anon namespace
} // namespace Opm
#endif // OPM_BLACKOILWELLS_HEADER_INCLUDED

582
opm/porsol/blackoil/ComponentTransport.hpp
View File

@ -1,582 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_COMPONENTTRANSPORT_HEADER_INCLUDED
#define OPM_COMPONENTTRANSPORT_HEADER_INCLUDED
#include <opm/porsol/blackoil/fluid/BlackoilDefs.hpp>
#include <opm/porsol/blackoil/BlackoilFluid.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <array>
#include <vector>
#include <iostream>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <limits.h>
#include <assert.h>
#include <stdio.h>
namespace Opm
{
template <class Grid, class Rock, class Fluid, class Wells>
class ExplicitCompositionalTransport : public BlackoilDefs
{
public:
/// @brief
/// Default constructor. Does nothing.
ExplicitCompositionalTransport()
: pgrid_(0), prock_(0), pfluid_(0), pwells_(0), ptrans_(0),
min_surfvol_threshold_(0.0),
single_step_only_(false),
min_vtime_(0.0)
{
}
void init(const Opm::ParameterGroup& param)
{
min_surfvol_threshold_ = param.getDefault("min_surfvol_threshold", min_surfvol_threshold_);
single_step_only_ = param.getDefault("single_step_only", single_step_only_);
min_vtime_ = param.getDefault("min_vtime", min_vtime_);
}
void setup(const Grid& grid,
const Rock& rock,
const Fluid& fluid,
const Wells& wells,
const std::vector<double>& face_trans,
const typename Grid::Vector& gravity)
{
pgrid_ = &grid;
prock_ = &rock;
pfluid_ = &fluid;
pwells_ = &wells;
ptrans_ = &face_trans;
gravity_ = gravity;
}
/// Return value is the time actually used, it may be smaller than dt if
/// we stop due to unacceptable volume discrepancy.
double transport(const PhaseVec& external_pressure,
const CompVec& external_composition,
const std::vector<double>& face_flux,
const std::vector<PhaseVec>& cell_pressure,
const std::vector<PhaseVec>& face_pressure,
const double dt,
const double voldisclimit,
std::vector<CompVec>& cell_z)
{
int num_cells = pgrid_->numCells();
std::vector<CompVec> comp_change;
std::vector<double> cell_outflux;
std::vector<double> cell_max_ff_deriv;
std::vector<double> cell_gravflux;
std::vector<double> cell_max_mob_deriv;
double cur_time = 0.0;
updateFluidProperties(cell_pressure, face_pressure, cell_z,
external_pressure, external_composition);
// if (!volumeDiscrepancyAcceptable(voldisclimit)) {
// return 0.0;
// }
std::vector<CompVec> cell_z_start;
std::cout << "Transport solver target time: " << dt << std::endl;
std::cout << " Step Stepsize Remaining time\n";
int count = 0;
while (cur_time < dt) {
cell_z_start = cell_z;
computeChange(face_flux, comp_change, cell_outflux, cell_max_ff_deriv);
double min_time = 1e100;
for (int cell = 0; cell < num_cells; ++cell) {
double pvol = prock_->porosity(cell)*pgrid_->cellVolume(cell);
double vtime = pvol/(cell_outflux[cell]*cell_max_ff_deriv[cell]);
double gtime = 1e100; // No working CFL for gravity yet.
double max_nonzero_time = 1e100;
for (int comp = 0; comp < numComponents; ++comp) {
if (comp_change[cell][comp] < 0.0) {
if (cell_z[cell][comp] > min_surfvol_threshold_) {
max_nonzero_time = std::min(max_nonzero_time,
-cell_z[cell][comp]*pvol/comp_change[cell][comp]);
} else {
comp_change[cell][comp] = 0.0;
cell_z[cell][comp] = 0.0;
}
}
}
vtime = std::max(vtime,min_vtime_);
double time = std::min(std::min(vtime, gtime), max_nonzero_time);
min_time = std::min(time, min_time);
}
min_time *= 0.49; // Semi-random CFL factor... \TODO rigorize
double step_time = dt - cur_time;
if (min_time < step_time) {
step_time = min_time;
cur_time += min_time;
} else {
cur_time = dt;
}
// Check if remaining number of steps is excessive.
if ((dt - cur_time)/step_time > 10000) {
std::cout << "Collapsing transport stepsize detected." << std::endl;
return cur_time;
}
std::cout.precision(10);
std::cout << std::setw(6) << count++
<< std::setw(24) << step_time
<< std::setw(24) << dt - cur_time << std::endl;
std::cout.precision(16);
for (int cell = 0; cell < num_cells; ++cell) {
double pv = pgrid_->cellVolume(cell)*prock_->porosity(cell);
comp_change[cell] *= (step_time/pv);
cell_z[cell] += comp_change[cell];
}
// After changing z, we recompute fluid properties.
updateFluidProperties(cell_pressure, face_pressure, cell_z,
external_pressure, external_composition);
bool ok = volumeDiscrepancyAcceptable(voldisclimit);
if (!ok) {
// Roll back to last ok step.
cell_z = cell_z_start;
cur_time -= step_time;
return cur_time;
}
if (single_step_only_) {
return cur_time;
}
}
return dt;
}
private: // Data
const Grid* pgrid_;
const Rock* prock_;
const Fluid* pfluid_;
const Wells* pwells_;
const std::vector<double>* ptrans_;
typename Grid::Vector gravity_;
struct AllTransportFluidData : public Opm::AllFluidData
{
std::vector<double> total_mobility;
std::vector<PhaseVec> fractional_flow;
template <class G, class R>
void computeNew(const G& grid,
const R& rock,
const BlackoilFluid& fluid,
const typename Grid::Vector gravity,
const std::vector<PhaseVec>& cell_pressure,
const std::vector<PhaseVec>& face_pressure,
const std::vector<CompVec>& cell_z,
const CompVec& bdy_z,
const double dt)
{
Opm::AllFluidData::computeNew(grid, rock, fluid, gravity,
cell_pressure, face_pressure,
cell_z, bdy_z, dt);
int num = grid.numCells();
total_mobility.resize(num);
fractional_flow.resize(num);
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
total_mobility[i] = 0.0;
for (int phase = 0; phase < numPhases; ++phase) {
total_mobility[i] += cell_data.mobility[i][phase];
}
fractional_flow[i] = cell_data.mobility[i];
fractional_flow[i] *= (1.0/total_mobility[i]);
}
}
};
AllTransportFluidData fluid_data_;
struct TransportFluidData
{
PhaseVec saturation;
PhaseVec mobility;
PhaseVec fractional_flow;
PhaseToCompMatrix phase_to_comp;
PhaseVec relperm;
PhaseVec viscosity;
};
TransportFluidData bdy_;
std::vector<int> perf_cells_;
std::vector<double> perf_flow_;
std::vector<TransportFluidData> perf_props_;
double min_surfvol_threshold_;
bool single_step_only_;
double min_vtime_;
private: // Methods
TransportFluidData computeProps(const PhaseVec& pressure,
const CompVec& composition)
{
BlackoilFluid::FluidState state = pfluid_->computeState(pressure, composition);
TransportFluidData data;
data.saturation = state.saturation_;
data.mobility = state.mobility_;
double total_mobility = 0.0;
for (int phase = 0; phase < numPhases; ++phase) {
total_mobility += state.mobility_[phase];
}
data.fractional_flow = state.mobility_;
data.fractional_flow /= total_mobility;
data.phase_to_comp = state.phase_to_comp_;
data.relperm = state.relperm_;
data.viscosity = state.viscosity_;
return data;
}
void updateFluidProperties(const std::vector<PhaseVec>& cell_pressure,
const std::vector<PhaseVec>& face_pressure,
const std::vector<CompVec>& cell_z,
const PhaseVec& external_pressure,
const CompVec& external_composition)
{
// Properties in reservoir.
const double dummy_dt = 1.0;
fluid_data_.computeNew(*pgrid_, *prock_, *pfluid_, gravity_,
cell_pressure, face_pressure, cell_z, external_composition, dummy_dt);
// Properties on boundary. \TODO no need to ever recompute this.
bdy_ = computeProps(external_pressure, external_composition);
// Properties at well perforations.
// \TODO only need to recompute this once per pressure update.
// No, that is false, at production perforations the cell z is
// used, which may change every step.
perf_cells_.clear();
perf_flow_.clear();
perf_props_.clear();
Wells::WellReport::report()->clearAll();
int num_cells = pgrid_->numCells();
for (int cell = 0; cell < num_cells; ++cell) {
double flow = pwells_->wellToReservoirFlux(cell);
if (flow != 0.0) {
perf_cells_.push_back(cell);
perf_flow_.push_back(flow);
// \TODO handle capillary in perforation pressure below?
PhaseVec well_pressure = flow > 0.0 ? PhaseVec(pwells_->perforationPressure(cell)) : cell_pressure[cell];
CompVec well_mixture = flow > 0.0 ? pwells_->injectionMixture(cell) : cell_z[cell];
perf_props_.push_back(computeProps(well_pressure, well_mixture));
Wells::WellReport::report()->perfPressure.push_back(pwells_->perforationPressure(cell));
Wells::WellReport::report()->cellPressure.push_back(cell_pressure[cell][0]);
Wells::WellReport::report()->cellId.push_back(cell);
}
}
}
void cellData(int cell, TransportFluidData& tfd) const
{
tfd.saturation = fluid_data_.cell_data.saturation[cell];
tfd.mobility = fluid_data_.cell_data.mobility[cell];
tfd.fractional_flow = fluid_data_.fractional_flow[cell];
tfd.phase_to_comp = fluid_data_.cell_data.state_matrix[cell];
tfd.relperm = fluid_data_.cell_data.relperm[cell];
tfd.viscosity = fluid_data_.cell_data.viscosity[cell];
}
bool volumeDiscrepancyAcceptable(const double voldisclimit) const
{
double rel_voldiscr = *std::max_element(fluid_data_.relvoldiscr.begin(), fluid_data_.relvoldiscr.end());
if (rel_voldiscr > voldisclimit) {
std::cout << " Relative volume discrepancy too large: " << rel_voldiscr << std::endl;
return false;
} else {
return true;
}
}
void computeChange(const std::vector<double>& face_flux,
std::vector<CompVec>& comp_change,
std::vector<double>& cell_outflux,
std::vector<double>& cell_max_ff_deriv)
{
const int num_cells = pgrid_->numCells();
comp_change.clear();
CompVec zero(0.0);
comp_change.resize(num_cells, zero);
cell_outflux.clear();
cell_outflux.resize(num_cells, 0.0);
cell_max_ff_deriv.clear();
cell_max_ff_deriv.resize(num_cells, 0.0);
CompVec surf_dens = pfluid_->surfaceDensities();
const int num_faces = pgrid_->numFaces();
std::vector<CompVec> face_change(num_faces);
std::vector<double> faces_max_ff_deriv(num_faces);
#pragma omp parallel for
for (int face = 0; face < num_faces; ++face) {
// Compute phase densities on face.
PhaseVec phase_dens(0.0);
for (int phase = 0; phase < numPhases; ++phase) {
const double* At = &fluid_data_.face_data.state_matrix[face][0][0]; // Already transposed since in Fortran order...
for (int comp = 0; comp < numPhases; ++comp) {
phase_dens[phase] += At[numPhases*phase + comp]*surf_dens[comp];
}
}
// Collect data from adjacent cells (or boundary).
int c[2];
TransportFluidData d[2];
for (int ix = 0; ix < 2; ++ix) {
c[ix] = pgrid_->faceCell(face, ix);
if (c[ix] >= 0) {
cellData(c[ix], d[ix]);
} else {
d[ix] = bdy_;
}
}
// Compute upwind directions.
int upwind_dir[numPhases] = { 0, 0, 0 };
PhaseVec vstar(face_flux[face]);
// double gravity_flux = gravity_*pgrid_->faceNormal(face)*pgrid_->faceArea(face);
typename Grid::Vector centroid_diff = c[0] >= 0 ? pgrid_->cellCentroid(c[0]) : pgrid_->faceCentroid(face);
centroid_diff -= c[1] >= 0 ? pgrid_->cellCentroid(c[1]) : pgrid_->faceCentroid(face);
double gravity_flux = gravity_*centroid_diff*ptrans_->operator[](face);
PhaseVec rho_star(phase_dens);
process_face(&(d[0].mobility[0]), &(d[1].mobility[0]),
&vstar[0], gravity_flux, numPhases, &rho_star[0], upwind_dir);
// Compute phase fluxes.
PhaseVec phase_mob;
double tot_mob = 0.0;
for (int phase = 0; phase < numPhases; ++phase) {
phase_mob[phase] = d[upwind_dir[phase] - 1].mobility[phase];
tot_mob += phase_mob[phase];
}
PhaseVec ff = phase_mob;
ff /= tot_mob;
PhaseVec phase_flux = ff;
phase_flux *= face_flux[face];
// Until we have proper bcs for transport, assume no gravity flow across bdys.
if (gravity_flux != 0.0 && c[0] >= 0 && c[1] >= 0) {
// Gravity contribution.
double omega = ff*phase_dens;
for (int phase = 0; phase < numPhases; ++phase) {
double gf = (phase_dens[phase] - omega)*gravity_flux;
phase_flux[phase] -= phase_mob[phase]*gf;
}
}
// Estimate max derivative of ff.
double face_max_ff_deriv = 0.0;
// Only using total flux upwinding for this purpose.
// Aim is to reproduce old results first. \TODO fix, include gravity.
int downwind_cell = c[upwind_dir[0]%2]; // Keep in mind that upwind_dir[] \in {1, 2}
if (downwind_cell >= 0) { // Only contribution on inflow and internal faces.
// Evaluating all functions at upwind viscosity.
// Added for this version.
PhaseVec upwind_viscosity = d[upwind_dir[0] - 1].viscosity;
PhaseVec upwind_sat = d[upwind_dir[0] - 1].saturation;
PhaseVec upwind_relperm = d[upwind_dir[0] - 1].relperm;
PhaseVec downwind_mob(0.0);
PhaseVec upwind_mob(0.0);
double downwind_totmob = 0.0;
double upwind_totmob = 0.0;
for (int phase = 0; phase < numPhases; ++phase) {
downwind_mob[phase] = fluid_data_.cell_data.relperm[downwind_cell][phase]/upwind_viscosity[phase];
downwind_totmob += downwind_mob[phase];
upwind_mob[phase] = upwind_relperm[phase]/upwind_viscosity[phase];
upwind_totmob += upwind_mob[phase];
}
PhaseVec downwind_ff = downwind_mob;
downwind_ff /= downwind_totmob;
PhaseVec upwind_ff = upwind_mob;
upwind_ff /= upwind_totmob;
PhaseVec ff_diff = upwind_ff;
ff_diff -= downwind_ff;
for (int phase = 0; phase < numPhases; ++phase) {
if (std::fabs(ff_diff[phase]) > 1e-10) {
if (face_flux[face] != 0.0) {
double ff_deriv = ff_diff[phase]/(upwind_sat[phase] - fluid_data_.cell_data.saturation[downwind_cell][phase]);
// assert(ff_deriv >= 0.0);
face_max_ff_deriv = std::max(face_max_ff_deriv, std::fabs(ff_deriv));
}
}
}
}
faces_max_ff_deriv[face] = face_max_ff_deriv;
// Compute z change.
CompVec change(0.0);
for (int phase = 0; phase < numPhases; ++phase) {
int upwind_ix = upwind_dir[phase] - 1; // Since process_face returns 1 or 2.
CompVec z_in_phase = d[upwind_ix].phase_to_comp[phase];
z_in_phase *= phase_flux[phase];
change += z_in_phase;
}
face_change[face] = change;
}
// Update output variables
#pragma omp parallel for
for (int cell = 0; cell < num_cells; ++cell) {
const int num_local_faces = pgrid_->numCellFaces(cell);
for (int local = 0; local < num_local_faces; ++local) {
int face = pgrid_->cellFace(cell,local);
if (cell == pgrid_->faceCell(face, 0)) {
comp_change[cell] -= face_change[face];
if (face_flux[face] >= 0.0) {
cell_outflux[cell] += std::fabs(face_flux[face]);
}
} else if (cell == pgrid_->faceCell(face, 1)) {
comp_change[cell] += face_change[face];
if (face_flux[face] < 0.0) {
cell_outflux[cell] += std::fabs(face_flux[face]);
}
}
cell_max_ff_deriv[cell] = std::max(cell_max_ff_deriv[cell],
faces_max_ff_deriv[face]);
}
}
// Done with all faces, now deal with well perforations.
int num_perf = perf_cells_.size();
for (int perf = 0; perf < num_perf; ++perf) {
int cell = perf_cells_[perf];
double flow = perf_flow_[perf];
assert(flow != 0.0);
const TransportFluidData& fl = perf_props_[perf];
// For injection, phase volumes depend on fractional
// flow of injection perforation, for production we
// use the fractional flow of the producing cell.
PhaseVec phase_flux = flow > 0.0 ? fl.fractional_flow : fluid_data_.fractional_flow[cell];
phase_flux *= flow;
// Conversion to mass flux is given at perforation state
// if injector, cell state if producer (this is ensured by
// updateFluidProperties()).
CompVec change(0.0);
for (int phase = 0; phase < numPhases; ++phase) {
CompVec z_in_phase = fl.phase_to_comp[phase];
z_in_phase *= phase_flux[phase];
change += z_in_phase;
}
comp_change[cell] += change;
Wells::WellReport::report()->massRate.push_back(change);
}
}
// Arguments are:
// [in] cmob1, cmob2: cell mobilities for adjacent cells
// [in, destroyed] vstar: total velocity (equal for all phases) on input, modified on output
// [in] gf: gravity flux
// [in] np: number of phases
// [in, destroyed] rho: density on face
// [out] ix: upwind cells for each phase (1 or 2)
static void
process_face(double *cmob1, double *cmob2, double *vstar, double gf,
int np, double *rho, int *ix)
{
int i,j,k,a,b,c;
double v, r, m;
for (i=0; i<np; ++i) {
/* ============= */
/* same sign */
/* r = max(rho)*/
j = 0; r = DBL_MIN;
for(k=0; k<np; ++k) { if (rho[k] > r) { r = rho[k]; j = k; } }
a = !(vstar[j]<0) && !(gf<0);
b = !(vstar[j]>0) && !(gf>0);
c = a || b;
if ( !c ) {
rho[j] = NAN;
v = vstar[j] - gf;
if (v < 0) { ix[j] = 2; m = cmob2[j]; }
else { ix[j] = 1; m = cmob1[j]; }
for (k=0; k<np; ++k) { vstar[k] -= m *(rho[k]-r)*gf;}
continue;
}
/* ============= */
/* opposite sign */
/* r = min(rho)*/
j = 0; r = DBL_MAX;
for(k=0; k<np; ++k) { if (rho[k] < r) { r = rho[k]; j = k; } }
if ( c ) {
rho[j] = NAN;
v = vstar[j] + gf;
if (v < 0) { ix[j] = 2; m = cmob2[j]; }
else { ix[j] = 1; m = cmob1[j]; }
for (k=0; k<np; ++k) { vstar[k] -= m *(rho[k]-r)*gf;}
continue;
}
}
}
static void
phase_upwind_directions(int number_of_faces, int *face_cells, double *dflux, double *gflux, int np,
double *cmobility, double *fdensity, int *ix)
{
int k,f;
int c1, c2;
double *cmob1, *cmob2;
double *fden = (double*)malloc(np * sizeof *fden);
double *vstar = (double*)malloc(np * sizeof *vstar);
double gf;
for (f=0; f<number_of_faces; ++f) {
gf = gflux[f];
for (k=0; k<np; ++k) {
vstar[k] = dflux[f];
fden [k] = fdensity[np*f+k];
}
c1 = face_cells[2*f + 0];
c2 = face_cells[2*f + 1];
if (c1 != -1) {cmob1 = cmobility+np*c1;}
else {cmob1 = cmobility+np*c2;}
if (c2 != -1) {cmob2 = cmobility+np*c2;}
else {cmob2 = cmobility+np*c1;}
process_face(cmob1, cmob2, vstar, gf, np, fden, ix+np*f);
}
free(fden);
free(vstar);
}
};
} // namespace Opm
#endif // OPM_COMPONENTTRANSPORT_HEADER_INCLUDED

500
opm/porsol/blackoil/co2fluid/BlackoilCo2PVT.hpp
View File

@ -1,500 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILCO2PVT_HEADER_INCLUDED
#define OPM_BLACKOILCO2PVT_HEADER_INCLUDED
#define OPM_BLACKOILPVT_HEADER_INCLUDED
#define OPM_DEPRECATED __attribute__((deprecated))
#define OPM_DEPRECATED_MSG(msg) __attribute__((deprecated))
#include <opm/material/fluidsystems/BrineCO2FluidSystem.hpp>
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
#include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp>
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/porsol/blackoil/co2fluid/benchmark3co2tables.hh>
#include <opm/porsol/blackoil/fluid/MiscibilityProps.hpp>
#include <opm/porsol/blackoil/fluid/BlackoilDefs.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
namespace Opm
{
class BlackoilCo2PVT : public BlackoilDefs
{
public:
typedef Opm::FluidSystems::BrineCO2</*Scalar=*/double, Opm::Benchmark3::CO2Tables> FluidSystem;
typedef Opm::CompositionalFluidState<double, FluidSystem> CompositionalFluidState;
void init(const Opm::Deck& deck);
void generateBlackOilTables(double temperature);
double getViscosity(double press,
const CompVec& surfvol,
PhaseIndex phase) const;
CompVec surfaceDensities() const;
double getSaturation(double press,
const CompVec& surfvol,
PhaseIndex phase) 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const;
void B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const;
void dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_B,
std::vector<PhaseVec>& output_dBdp) const;
void R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const;
void dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_R,
std::vector<PhaseVec>& output_dRdp) const;
private:
CompVec surfaceDensities_;
FluidSystem brineCo2_;
double temperature_;
struct SubState
{
Dune::FieldVector<double,2> density;
Dune::FieldMatrix<double,2,2> massfrac;
Dune::FieldVector<double,2> phaseVolume;
Dune::FieldVector<double,2> phaseViscosity;
double saturation;
};
void computeState(SubState& ss, double zBrine, double zCO2, double pressure) const;
enum {
wPhase = FluidSystem::liquidPhaseIdx,
nPhase = FluidSystem::gasPhaseIdx,
wComp = FluidSystem::BrineIdx,
nComp = FluidSystem::CO2Idx
};
}; // class BlackoilCo2PVT
// ------------ Method implementations --------------
void BlackoilCo2PVT::init(const Opm::Deck& /* deck */)
{
surfaceDensities_[Water] = 1000.;
surfaceDensities_[Gas] = 2.0;
surfaceDensities_[Oil] = 1000.;
temperature_ = 300.;
brineCo2_.init();
}
void BlackoilCo2PVT::generateBlackOilTables(double temperature)
{
std::cout << "\n Generating pvt tables for the eclipse black oil formulation\n using the oil component as brine and the gas component as co_2." << std::endl;
if (std::fabs(temperature-400.) > 100.0) {
std::cout << "T=" << temperature << " is outside the allowed range [300,500] Kelvin" << std::endl;
exit(-1);
}
temperature_ = temperature;
CompVec z;
z[Water] = 0.0;
z[Oil] = 1.0;
z[Gas] = 1.0e6;
std::ofstream black("blackoil_pvt");
black.precision(8);
black << std::fixed << std::showpoint;
const double pMin=150.e5;
const double pMax=400.e5;
const unsigned int np=11;
std::vector<double> pValue(np+1,0.0);
std::vector<double> rs(np+1,0.0);
pValue[0] = 101330.0;
rs[0] = 0.0;
// Buble points
z[Gas] = 1.0e4;
for (unsigned int i=0; i<np; ++i) {
pValue[i+1] = pMin + i*(pMax-pMin)/(np-1);
rs[i+1] = R(pValue[i+1], z, Liquid);
}
const unsigned int np_fill_in = 10;
const double dr = (rs[1] - rs[0])/(np_fill_in+1);
const double dp = (pValue[1] - pValue[0])/(np_fill_in+1);
rs.insert(rs.begin()+1,np_fill_in,0.0);
pValue.insert(pValue.begin()+1,np_fill_in,0.0);
for (unsigned int i=1; i<=np_fill_in; ++i) {
rs[i] = rs[i-1] + dr;
pValue[i] = pValue[i-1] + dp;
}
// Brine with dissolved co_2 ("live oil")
black << "PVTO\n";
black << "-- Rs Pbub Bo Vo\n";
black << "-- sm3/sm3 barsa rm3/sm3 cP\n";
// Undersaturated
for (unsigned int i=0; i<np+np_fill_in; ++i) {
z[Gas] = rs[i];
black << std::setw(14) << rs[i];
for (unsigned int j=i; j<np+1+np_fill_in; ++j) {
if (j<=np_fill_in) {
if (j==i) black << std::setw(14) << pValue[j]*1.e-5 << std::setw(14) << 1.0-j*0.001 << std::setw(14) << 1.06499;
continue;
}
if (j>i) black << std::endl << std::setw(14) << ' ';
black << std::setw(14) << pValue[j]*1.e-5
<< std::setw(14) << B(pValue[j], z, Liquid)
<< std::setw(14) << getViscosity(pValue[j], z, Liquid)*1.e3;
}
black << " /" << std::endl;
}
black << "/ " << std::endl;
// We provide tables for co_2 both with and without dissolved water:
// Co_2 neglecting dissolved water ("dry gas")
black << "\nPVDG\n";
black << "-- Pg Bg Vg\n";
black << "-- barsa rm3/sm3 cP\n";
for (unsigned int i=0; i<np; ++i) {
z[Oil] = 0.0;
z[Gas] = 1.0;
black << std::setw(14) << pValue[i+np_fill_in+1]*1.e-5
<< std::setw(14) << B(pValue[i+np_fill_in+1], z, Vapour)
<< std::setw(14) << getViscosity(pValue[i+np_fill_in+1], z, Vapour)*1.e3
<< std::endl;
}
black << "/ " << std::endl;
// Co_2 with dissolved water ("wet gas")
black << "\nPVTG\n";
black << "-- Pg Rv Bg Vg\n";
black << "-- barsa sm3/sm3 rm3/sm3 cP\n";
for (unsigned int i=0; i<np; ++i) {
z[Oil] = 1000.0;
z[Gas] = 1.0;
black << std::setw(14) << pValue[i+np_fill_in+1]*1.e-5
<< std::setw(14) << R(pValue[i+np_fill_in+1], z, Vapour)
<< std::setw(14) << B(pValue[i+np_fill_in+1], z, Vapour)
<< std::setw(14) << getViscosity(pValue[i+np_fill_in+1], z, Vapour)*1.e3
<< std::endl;
z[Oil] = 0.0;
black << std::setw(14) << ' '
<< std::setw(14) << R(pValue[i+np_fill_in+1], z, Vapour)
<< std::setw(14) << B(pValue[i+np_fill_in+1], z, Vapour)
<< std::setw(14) << getViscosity(pValue[i+np_fill_in+1], z, Vapour)*1.e3
<< " /" << std::endl;
}
black << "/ " << std::endl;
black << std::endl;
std::cout << " Pvt tables for temperature=" << temperature << " Kelvin is written to file blackoil_pvt. " << std::endl;
std::cout << " NOTE that the file contains tables for both PVDG (dry gas) and PVTG (wet gas)." << std::endl;
}
BlackoilCo2PVT::CompVec BlackoilCo2PVT::surfaceDensities() const
{
return surfaceDensities_;
}
double BlackoilCo2PVT::getViscosity(double press, const CompVec& surfvol, PhaseIndex phase) const
{
SubState ss;
computeState(ss, surfvol[Oil], surfvol[Gas], press);
switch(phase) {
case Aqua: return 1.0e-10;
case Liquid: return ss.phaseViscosity[wPhase];
case Vapour: return ss.phaseViscosity[nPhase];
};
return 1.0e-10;
}
double BlackoilCo2PVT::getSaturation(double press, const CompVec& surfvol, PhaseIndex phase) const
{
SubState ss;
computeState(ss, surfvol[Oil], surfvol[Gas], press);
switch(phase) {
case Aqua: return 0.0;
case Liquid: return ss.saturation;
case Vapour: return 1.0 - ss.saturation;
};
return 0.0;
}
double BlackoilCo2PVT::B(double press, const CompVec& surfvol, PhaseIndex phase) const
{
SubState ss;
computeState(ss, surfvol[Oil], surfvol[Gas], press);
switch(phase) {
case Aqua: return 1.0;
case Liquid: return surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
case Vapour: return surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
};
return 1.0;
}
double BlackoilCo2PVT::dBdp(double press, const CompVec& surfvol, PhaseIndex phase) const
{
const double dp = 100.;
return (B(press+dp,surfvol,phase)-B(press,surfvol,phase))/dp;
}
double BlackoilCo2PVT::R(double press, const CompVec& surfvol, PhaseIndex phase) const
{
SubState ss;
computeState(ss, surfvol[Oil], surfvol[Gas], press);
switch(phase) {
case Aqua: return 0.0;
case Liquid: return (ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10);
case Vapour: return (ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10);
};
return 0.0;
}
double BlackoilCo2PVT::dRdp(double press, const CompVec& surfvol, PhaseIndex phase) const
{
const double dp = 100.;
return (R(press+dp,surfvol,phase)-R(press,surfvol,phase))/dp;
}
void BlackoilCo2PVT::getViscosity(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const
{
int num = pressures.size();
output.resize(num);
SubState ss;
for (int i = 0; i < num; ++i) {
computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
output[i][Aqua] = 1.0e-10;
output[i][Liquid] = ss.phaseViscosity[wPhase];
output[i][Vapour] = ss.phaseViscosity[nPhase];
}
}
void BlackoilCo2PVT::B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const
{
int num = pressures.size();
output.resize(num);
SubState ss;
for (int i = 0; i < num; ++i) {
computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
output[i][Aqua] = 1.0;
output[i][Liquid] = surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
output[i][Vapour] = surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
}
}
void BlackoilCo2PVT::dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_B,
std::vector<PhaseVec>& output_dBdp) const
{
int num = pressures.size();
output_B.resize(num);
output_dBdp.resize(num);
SubState ss;
const double dp = 100.;
for (int i = 0; i < num; ++i) {
computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
output_B[i][Aqua] = 1.0;
output_B[i][Liquid] = surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10);
output_B[i][Vapour] = surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10);
computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]+dp);
output_dBdp[i][Aqua] = 0.0;
output_dBdp[i][Liquid] = (surfaceDensities_[Oil]/(ss.massfrac[wPhase][wComp]*ss.density[wPhase]+1.0e-10) - output_B[i][Liquid])/dp;
output_dBdp[i][Vapour] = (surfaceDensities_[Gas]/(ss.massfrac[nPhase][nComp]*ss.density[nPhase]+1.0e-10) - output_B[i][Vapour])/dp;
}
}
void BlackoilCo2PVT::R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const
{
int num = pressures.size();
output.resize(num);
SubState ss;
for (int i = 0; i < num; ++i) {
computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
output[i][Aqua] = 0.0;
output[i][Liquid] = (ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10);
output[i][Vapour] = (ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10);
}
}
void BlackoilCo2PVT::dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_R,
std::vector<PhaseVec>& output_dRdp) const
{
int num = pressures.size();
output_R.resize(num);
output_dRdp.resize(num);
SubState ss;
const double dp = 100.;
for (int i = 0; i < num; ++i) {
computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]);
output_R[i][Aqua] = 0.0;
output_R[i][Liquid] = (ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10);
output_R[i][Vapour] = (ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10);
computeState(ss, surfvol[i][Oil], surfvol[i][Gas], pressures[i][Liquid]+dp);
output_dRdp[i][Aqua] = 0.0;
output_dRdp[i][Liquid] = ((ss.massfrac[wPhase][nComp]*surfaceDensities_[Oil])/(ss.massfrac[wPhase][wComp]*surfaceDensities_[Gas]+1.0e-10) - output_R[i][Liquid])/dp;
output_dRdp[i][Vapour] = ((ss.massfrac[nPhase][wComp]*surfaceDensities_[Gas])/(ss.massfrac[nPhase][nComp]*surfaceDensities_[Oil]+1.0e-10) - output_R[i][Vapour])/dp;
}
}
void BlackoilCo2PVT::computeState(BlackoilCo2PVT::SubState& ss, double zBrine, double zCO2, double pressure) const
{
CompositionalFluidState fluidState;
fluidState.setTemperature(temperature_);
fluidState.setPressure(wPhase, pressure);
fluidState.setPressure(nPhase, pressure);
double massH20 = surfaceDensities_[Oil]*zBrine;
double massCO2 = surfaceDensities_[Gas]*zCO2;
// A priori, assume presence of both phases
FluidSystem::ParameterCache<double> paramCache;
typedef Opm::MiscibleMultiPhaseComposition</*Scalar=*/double, FluidSystem> MMPC;
MMPC::solve(fluidState, paramCache, /*setViscosity=*/false, /*setEnthalpy=*/false);
ss.density[wPhase] = fluidState.density(wPhase);
ss.density[nPhase] = fluidState.density(nPhase);
ss.massfrac[wPhase][nComp] = fluidState.massFraction(wPhase, nComp);
ss.massfrac[nPhase][wComp] = fluidState.massFraction(nPhase, wComp);
ss.massfrac[wPhase][wComp] = 1.0 - ss.massfrac[wPhase][nComp];
ss.massfrac[nPhase][nComp] = 1.0 - ss.massfrac[nPhase][wComp];
double detX = ss.massfrac[wPhase][wComp]*ss.massfrac[nPhase][nComp]-ss.massfrac[wPhase][nComp]*ss.massfrac[nPhase][wComp];
ss.phaseVolume[wPhase] = (massH20*ss.massfrac[nPhase][nComp] - massCO2*ss.massfrac[nPhase][wComp])/(ss.density[wPhase]*detX);
ss.phaseVolume[nPhase] = (massCO2*ss.massfrac[wPhase][wComp] - massH20*ss.massfrac[wPhase][nComp])/(ss.density[nPhase]*detX);
// Determine number of phase
if (ss.phaseVolume[wPhase] > 0.0 && ss.phaseVolume[nPhase] > 0.0) { // Both phases
ss.saturation = ss.phaseVolume[wPhase]/(ss.phaseVolume[wPhase]+ss.phaseVolume[nPhase]);
fluidState.setSaturation(wPhase, ss.saturation);
fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
}
else if (ss.phaseVolume[wPhase] <= 0.0) { // Wetting phase only
ss.saturation = 0.0;
// Gas phase:
ss.massfrac[nPhase][nComp] = massCO2/(massCO2+massH20);
ss.massfrac[nPhase][wComp] = 1.0 - ss.massfrac[nPhase][nComp];
double M1 = FluidSystem::molarMass(wComp);
double M2 = FluidSystem::molarMass(nComp);
double avgMolarMass = M1*M2/(M2 + ss.massfrac[nPhase][nComp]*(M1 - M2));
fluidState.setMoleFraction(nPhase, nComp, ss.massfrac[nPhase][nComp]*avgMolarMass/M2);
fluidState.setMoleFraction(nPhase, wComp, ss.massfrac[nPhase][wComp]*avgMolarMass/M1);
ss.density[nPhase] = brineCo2_.density(fluidState, paramCache, nPhase);
fluidState.setDensity(nPhase, ss.density[nPhase]);
ss.phaseVolume[nPhase] = (massH20+massCO2)/ss.density[nPhase];
fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
// Virtual properties of non-existing liquid phase:
paramCache.updatePhase(fluidState, /*phaseIdx=*/nPhase);
typedef Opm::ComputeFromReferencePhase</*Scalar=*/double, FluidSystem> CFRP;
CFRP::solve(fluidState,
paramCache,
/*refPhaseIdx=*/nPhase,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
ss.massfrac[wPhase][wComp] = fluidState.massFraction(wPhase, wComp);
ss.massfrac[wPhase][nComp] = fluidState.massFraction(wPhase, nComp);
ss.density[wPhase] = fluidState.density(wPhase);
ss.phaseVolume[wPhase] = 0.0;
fluidState.setSaturation(wPhase, ss.saturation);
}
else if (ss.phaseVolume[nPhase] <= 0.0) { // Non-wetting phase only
ss.saturation = 1.0;
// Liquid phase:
ss.massfrac[wPhase][wComp] = massH20/(massCO2+massH20);
ss.massfrac[wPhase][nComp] = 1.0 - ss.massfrac[wPhase][wComp];
double M1 = FluidSystem::molarMass(wComp);
double M2 = FluidSystem::molarMass(nComp);
double avgMolarMass = M1*M2/(M2 + ss.massfrac[wPhase][nComp]*(M1 - M2));
fluidState.setMoleFraction(wPhase, nComp, ss.massfrac[wPhase][nComp]*avgMolarMass/M2);
fluidState.setMoleFraction(wPhase, wComp, ss.massfrac[wPhase][wComp]*avgMolarMass/M1);
ss.density[wPhase] = brineCo2_.density(fluidState, paramCache, wPhase);
fluidState.setDensity(wPhase, ss.density[wPhase]);
ss.phaseVolume[wPhase] = (massH20+massCO2)/ss.density[wPhase];
fluidState.setSaturation(wPhase, ss.saturation);
// Virtual properties of non-existing gas phase:
paramCache.updatePhase(fluidState, /*phaseIdx=*/nPhase);
typedef ComputeFromReferencePhase</*Scalar=*/double, FluidSystem> CFRP;
CFRP::solve(fluidState,
paramCache,
/*refPhaseIdx=*/wPhase,
/*setViscosity=*/false,
/*setEnthalpy=*/false);
ss.massfrac[nPhase][nComp] = fluidState.massFraction(nPhase, nComp);
ss.massfrac[nPhase][wComp] = fluidState.massFraction(nPhase, wComp);
ss.density[nPhase] = fluidState.density(nPhase);
ss.phaseVolume[nPhase] = 0.0;
fluidState.setSaturation(nPhase, 1.0 - ss.saturation);
}
ss.phaseViscosity[wPhase] = brineCo2_.viscosity(fluidState, paramCache, wPhase);
ss.phaseViscosity[nPhase] = brineCo2_.viscosity(fluidState, paramCache, nPhase);
}
typedef BlackoilCo2PVT BlackoilPVT;
} // Opm
#endif // OPM_BLACKOILCO2PVT_HEADER_INCLUDED

44
opm/porsol/blackoil/co2fluid/benchmark3co2tables.hh
View File

@ -1,44 +0,0 @@
/*****************************************************************************
* Copyright (C) 2009-2010 by Melanie Darcis *
* Copyright (C) 2010 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* This program 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 2 of the License, or *
* (at your option) any later version. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/**
* \file
*
* \brief Provides the class with the tabulated values of CO2 for the
* benchmark3 problem
*/
#ifndef OPM_BENCHMARK3_CO2TABLES_HH
#define OPM_BENCHMARK3_CO2TABLES_HH
#include <opm/material/common/UniformTabulated2DFunction.hpp>
#include <assert.h>
namespace Opm
{
namespace Benchmark3
{
// the real work is done by some external program which provides
// ready-to-use tables.
#include "benchmark3co2values.inc"
}
}
#endif

25298
opm/porsol/blackoil/co2fluid/benchmark3co2values.inc
View File

File diff suppressed because it is too large Load Diff

165
opm/porsol/blackoil/fluid/BlackoilComponent.hpp
View File

@ -1,165 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
#define OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
#include <dune/common/stdstreams.hh>
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 Scalar>
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

52
opm/porsol/blackoil/fluid/BlackoilDefs.hpp
View File

@ -1,52 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILDEFS_HEADER_INCLUDED
#define OPM_BLACKOILDEFS_HEADER_INCLUDED
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
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<Scalar, numComponents> CompVec;
typedef Dune::FieldVector<Scalar, numPhases> PhaseVec;
static_assert(int(numComponents) == int(numPhases), "");
typedef Dune::FieldMatrix<Scalar, numComponents, numPhases> PhaseToCompMatrix;
typedef Dune::FieldMatrix<Scalar, numPhases, numPhases> PhaseJacobian;
};
} // namespace Opm
#endif // OPM_BLACKOILDEFS_HEADER_INCLUDED

220
opm/porsol/blackoil/fluid/BlackoilPVT.cpp
View File

@ -1,220 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include "BlackoilPVT.hpp"
#include <opm/parser/eclipse/Units/Units.hpp>
#include "MiscibilityDead.hpp"
#include "MiscibilityLiveOil.hpp"
#include "MiscibilityLiveGas.hpp"
#include "MiscibilityWater.hpp"
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdoTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdgTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
using namespace Opm;
namespace Opm
{
void BlackoilPVT::init(const Deck& deck)
{
Opm::ParseContext parseContext;
EclipseState eclipseState( deck , parseContext);
region_number_ = 0;
// Surface densities. Accounting for different orders in eclipse and our code.
if (deck.hasKeyword("DENSITY")) {
const auto& densityRecord =
deck.getKeyword("DENSITY").getRecord(/*regionIdx=*/0);
densities_[Aqua] = densityRecord.getItem("WATER").getSIDouble(0);
densities_[Vapour] = densityRecord.getItem("GAS").getSIDouble(0);
densities_[Liquid] = densityRecord.getItem("OIL").getSIDouble(0);
} else {
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
}
// Water PVT
if (deck.hasKeyword("PVTW")) {
water_props_.reset(new MiscibilityWater(deck.getKeyword("PVTW")));
} else {
water_props_.reset(new MiscibilityWater(0.5*Opm::prefix::centi*Opm::unit::Poise)); // Eclipse 100 default
}
// Oil PVT
const auto& tables = eclipseState.getTableManager();
const auto& pvdoTables = tables.getPvdoTables();
const auto& pvtoTables = tables.getPvtoTables();
if (!pvdoTables.empty()) {
const auto& pvdoTable = pvdoTables.getTable<PvdoTable>(0);
oil_props_.reset(new MiscibilityDead(pvdoTable));
} else if (pvtoTables.empty()) {
// PVTOTables is a std::vector<>
const auto& pvtoTable = pvtoTables[0];
oil_props_.reset(new MiscibilityLiveOil(pvtoTable));
} else if (deck.hasKeyword("PVCDO")) {
auto *misc_water = new MiscibilityWater(0);
misc_water->initFromPvcdo(deck.getKeyword("PVCDO"));
oil_props_.reset(misc_water);
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDO and PVTO\n");
}
// Gas PVT
const auto& pvdgTables = tables.getPvdgTables();
const auto& pvtgTables = tables.getPvtgTables();
if (!pvdgTables.empty()) {
const auto& pvdgTable = pvdgTables.getTable<PvdgTable>(0);
gas_props_.reset(new MiscibilityDead(pvdgTable));
} else if (pvtgTables.empty()) {
gas_props_.reset(new MiscibilityLiveGas(pvtgTables[0]));
} else {
OPM_THROW(std::runtime_error, "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:
OPM_THROW(std::runtime_error, "Unknown phase accessed: " << phase);
}
}
void BlackoilPVT::getViscosity(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_B,
std::vector<PhaseVec>& 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_R,
std::vector<PhaseVec>& 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

88
opm/porsol/blackoil/fluid/BlackoilPVT.hpp
View File

@ -1,88 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILPVT_HEADER_INCLUDED
#define OPM_BLACKOILPVT_HEADER_INCLUDED
#include "MiscibilityProps.hpp"
#include "BlackoilDefs.hpp"
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <boost/scoped_ptr.hpp>
#include <string>
namespace Opm
{
class BlackoilPVT : public BlackoilDefs
{
public:
void init(const Deck& deck);
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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const;
void B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const;
void dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_B,
std::vector<PhaseVec>& output_dBdp) const;
void R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output) const;
void dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
std::vector<PhaseVec>& output_R,
std::vector<PhaseVec>& output_dRdp) const;
private:
int region_number_;
const MiscibilityProps& propsForPhase(PhaseIndex phase) const;
boost::scoped_ptr<MiscibilityProps> water_props_;
boost::scoped_ptr<MiscibilityProps> oil_props_;
boost::scoped_ptr<MiscibilityProps> gas_props_;
CompVec densities_;
mutable std::vector<double> data1_;
mutable std::vector<double> data2_;
};
}
#endif // OPM_BLACKOILPVT_HEADER_INCLUDED

192
opm/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp
View File

@ -1,192 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
#define OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SgofTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include "BlackoilDefs.hpp"
#include <iostream>
#include <stdexcept>
namespace Opm
{
// Forward declaration needed by associated parameters class.
template <class ScalarT, class ParamsT>
class FluidMatrixInteractionBlackoil;
template <class ScalarT>
class FluidMatrixInteractionBlackoilParams
{
public:
typedef ScalarT Scalar;
void init(const Opm::Deck& deck)
{
ParseContext parseContext;
EclipseState eclipseState(deck , parseContext);
// Extract input data.
const auto& tables = eclipseState.getTableManager();
const auto& swofTable = tables.getSwofTables().getTable<SwofTable>(0);
const auto& sgofTable = tables.getSgofTables().getTable<SgofTable>(0);
std::vector<double> sw = swofTable.getSwColumn().vectorCopy();
std::vector<double> krw = swofTable.getKrwColumn().vectorCopy();
std::vector<double> krow = swofTable.getKrowColumn().vectorCopy();
std::vector<double> pcow = swofTable.getPcowColumn().vectorCopy();
std::vector<double> sg = sgofTable.getSgColumn().vectorCopy();
std::vector<double> krg = sgofTable.getKrgColumn().vectorCopy();
std::vector<double> krog = sgofTable.getKrogColumn().vectorCopy();
std::vector<double> pcog = sgofTable.getPcogColumn().vectorCopy();
// 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 <class S, class P>
friend class FluidMatrixInteractionBlackoil;
Opm::utils::UniformTableLinear<Scalar> krw_;
Opm::utils::UniformTableLinear<Scalar> krow_;
Opm::utils::UniformTableLinear<Scalar> pcow_;
Opm::utils::UniformTableLinear<Scalar> krg_;
Opm::utils::UniformTableLinear<Scalar> krog_;
Opm::utils::UniformTableLinear<Scalar> pcog_;
Scalar krocw_; // = krow_(s_wc)
};
/*!
* \ingroup material
*
* \brief Capillary pressures and relative permeabilities for a black oil system.
*/
template <class ScalarT, class ParamsT = FluidMatrixInteractionBlackoilParams<ScalarT> >
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 <class pcContainerT, class SatContainerT>
static void pC(pcContainerT &pc,
const Params &params,
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 relative permeability of all phases.
*/
template <class krContainerT, class SatContainerT>
static void kr(krContainerT &kr,
const Params &params,
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 <class krContainerT, class SatContainerT>
static void dkr(krContainerT &dkr,
const Params &params,
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

91
opm/porsol/blackoil/fluid/FluidStateBlackoil.hpp
View File

@ -1,91 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#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<CompVec> surface_volume_density; // z
std::vector<PhaseVec> phase_pressure; // p
// Variables from PVT functions.
std::vector<PhaseVec> formation_volume_factor; // B
std::vector<PhaseVec> formation_volume_factor_deriv; // dB/dp
std::vector<PhaseVec> solution_factor; // R
std::vector<PhaseVec> solution_factor_deriv; // dR/dp
std::vector<PhaseVec> viscosity; // mu
// Variables computed from PVT data.
// The A matrices are all in Fortran order (or, equivalently,
// we store the transposes).
std::vector<PhaseToCompMatrix> state_matrix; // A' = (RB^{-1})'
std::vector<PhaseVec> phase_volume_density; // u
std::vector<Scalar> total_phase_volume_density; // sum(u)
std::vector<PhaseVec> saturation; // s = u/sum(u)
std::vector<PhaseVec> phase_compressibility; // c
std::vector<Scalar> total_compressibility; // cT
std::vector<Scalar> experimental_term; // ex = sum(Ai*dA*Ai*z)
// Variables computed from saturation.
std::vector<PhaseVec> relperm; // kr
std::vector<PhaseJacobian> relperm_deriv; // dkr/ds
std::vector<PhaseVec> mobility; // lambda
std::vector<PhaseJacobian> mobility_deriv; // dlambda/ds
};
} // end namespace Opm
#endif // OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED

199
opm/porsol/blackoil/fluid/MiscibilityDead.cpp
View File

@ -1,199 +0,0 @@
//===========================================================================
//
// File: MiscibilityDead.cpp
//
// Created: Wed Feb 10 09:06:04 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <algorithm>
#include "MiscibilityDead.hpp"
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdgTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdoTable.hpp>
#include <boost/lexical_cast.hpp>
#include <string>
#include <fstream>
using namespace std;
using namespace Opm;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityDead::MiscibilityDead(const Opm::PvdgTable& pvdgTable)
{
// Copy data
std::vector<double> press = pvdgTable.getPressureColumn().vectorCopy( );
std::vector<double> B_inv = pvdgTable.getFormationFactorColumn().vectorCopy( );
std::vector<double> visc = pvdgTable.getViscosityColumn().vectorCopy( );
const int sz = B_inv.size();
for (int i = 0; i < sz; ++i) {
B_inv[i] = 1.0 / B_inv[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<std::string>(count++)).c_str());
// os.precision(15);
// os << "1/B\n\n" << one_over_B_
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
}
/// Constructor
MiscibilityDead::MiscibilityDead(const Opm::PvdoTable& pvdoTable)
{
// Copy data
std::vector<double> press = pvdoTable.getPressureColumn().vectorCopy( );
std::vector<double> B_inv = pvdoTable.getFormationFactorColumn().vectorCopy( );
std::vector<double> visc = pvdoTable.getViscosityColumn().vectorCopy( );
const int sz = B_inv.size();
for (int i = 0; i < sz; ++i) {
B_inv[i] = 1.0 / B_inv[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<std::string>(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<PhaseVec>& pressures,
const std::vector<CompVec>&,
int phase,
std::vector<double>& 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<PhaseVec>& pressures,
const std::vector<CompVec>&,
int phase,
std::vector<double>& 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvols,
int phase,
std::vector<double>& output_B,
std::vector<double>& 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<PhaseVec>& pressures,
const std::vector<CompVec>&,
int,
std::vector<double>& 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<PhaseVec>& pressures,
const std::vector<CompVec>&,
int,
std::vector<double>& output_R,
std::vector<double>& 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);
}
}

94
opm/porsol/blackoil/fluid/MiscibilityDead.hpp
View File

@ -1,94 +0,0 @@
//===========================================================================
//
// File: MiscibilityDead.hpp
//
// Created: Wed Feb 10 09:05:47 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#ifndef SINTEF_MISCIBILITYDEAD_HEADER
#define SINTEF_MISCIBILITYDEAD_HEADER
/** Class for immiscible dead oil and dry gas.
* Detailed description.
*/
#include "MiscibilityProps.hpp"
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
namespace Opm
{
class PvdoTable;
class PvdgTable;
class MiscibilityDead : public MiscibilityProps
{
public:
MiscibilityDead(const PvdoTable& pvdoTable);
MiscibilityDead(const PvdgTable& pvdgTable);
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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_B,
std::vector<double>& output_dBdp) const;
virtual void R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_R,
std::vector<double>& output_dRdp) const;
private:
// PVT properties of dry gas or dead oil
Opm::utils::UniformTableLinear<double> one_over_B_;
Opm::utils::UniformTableLinear<double> viscosity_;
};
}
#endif // SINTEF_MISCIBILITYDEAD_HEADER

275
opm/porsol/blackoil/fluid/MiscibilityLiveGas.cpp
View File

@ -1,275 +0,0 @@
//===========================================================================
//
// File: MiscibilityLiveGas.cpp
//
// Created: Wed Feb 10 09:21:53 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <algorithm>
#include "MiscibilityLiveGas.hpp"
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvtgTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
using namespace std;
using namespace Opm;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityLiveGas::MiscibilityLiveGas(const Opm::PvtgTable& pvtgTable)
{
// GAS, PVTG
const auto &saturatedPvtgTable = pvtgTable.getSaturatedTable();
saturated_gas_table_.resize(4);
saturated_gas_table_[0] = saturatedPvtgTable.getColumn("PG").vectorCopy();
saturated_gas_table_[1] = saturatedPvtgTable.getColumn("BG").vectorCopy();
saturated_gas_table_[2] = saturatedPvtgTable.getColumn("MUG").vectorCopy();
saturated_gas_table_[3] = saturatedPvtgTable.getColumn("RV").vectorCopy();
int sz = saturated_gas_table_[0].size();
undersat_gas_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
const auto &undersaturatedPvtgTable = pvtgTable.getUnderSaturatedTable(i);
undersat_gas_tables_[i][0] = undersaturatedPvtgTable.getColumn("RV").vectorCopy();
undersat_gas_tables_[i][1] = undersaturatedPvtgTable.getColumn("BG").vectorCopy();
undersat_gas_tables_[i][2] = undersaturatedPvtgTable.getColumn("MUG").vectorCopy();
}
}
// Destructor
MiscibilityLiveGas::~MiscibilityLiveGas()
{
}
double MiscibilityLiveGas::getViscosity(int /*region*/, double press, const surfvol_t& surfvol) const
{
return miscible_gas(press, surfvol, 2, false);
}
void MiscibilityLiveGas::getViscosity(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& 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 = linearInterpolation(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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& 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 = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (R < maxR ) { // Saturated case
return linearInterpolationDerivative(saturated_gas_table_[0],
saturated_gas_table_[3],
press);
} else {
return 0.0; // Undersaturated case
}
}
void MiscibilityLiveGas::dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_R,
std::vector<double>& 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& 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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_B,
std::vector<double>& 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 = linearInterpolation(saturated_gas_table_[0],
saturated_gas_table_[3], press,
section);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (deriv) {
if (R < maxR ) { // Saturated case
return linearInterpolationDerivative(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 =
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolation(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 linearInterpolation(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 linearInterpolation(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 linearInterpolation(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 =
linearInterpolation(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolation(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

94
opm/porsol/blackoil/fluid/MiscibilityLiveGas.hpp
View File

@ -1,94 +0,0 @@
//===========================================================================
//
// File: MiscibilityLiveGas.hpp
//
// Created: Wed Feb 10 09:21:26 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#ifndef SINTEF_MISCIBILITYLIVEGAS_HEADER
#define SINTEF_MISCIBILITYLIVEGAS_HEADER
/** Class for miscible wet gas.
* Detailed description.
*/
#include "MiscibilityProps.hpp"
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
namespace Opm
{
class PvtgTable;
class MiscibilityLiveGas : public MiscibilityProps
{
public:
MiscibilityLiveGas(const Opm::PvtgTable& pvtgTable);
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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_B,
std::vector<double>& output_dBdp) const;
virtual void R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_R,
std::vector<double>& 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<std::vector<double> > saturated_gas_table_;
std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;
};
}
#endif // SINTEF_MISCIBILITYLIVEGAS_HEADER

396
opm/porsol/blackoil/fluid/MiscibilityLiveOil.cpp
View File

@ -1,396 +0,0 @@
//===========================================================================
//
// File: MiscibiltyLiveOil.cpp
//
// Created: Wed Feb 10 09:08:25 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <algorithm>
#include "MiscibilityLiveOil.hpp"
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvtoTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
using namespace std;
using namespace Opm;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityLiveOil::MiscibilityLiveOil(const PvtoTable& pvtoTable)
{
const auto& saturatedPvtoTable = pvtoTable.getSaturatedTable();
saturated_oil_table_.resize(4);
saturated_oil_table_[0] = saturatedPvtoTable.getColumn("P").vectorCopy();
saturated_oil_table_[1] = saturatedPvtoTable.getColumn("BO").vectorCopy();
saturated_oil_table_[2] = saturatedPvtoTable.getColumn("MU").vectorCopy();
saturated_oil_table_[3] = saturatedPvtoTable.getColumn("RS").vectorCopy();
// store the inverse of Bo...
int sz = saturated_oil_table_[1].size();
for (size_t i=0; i < saturated_oil_table_[1].size(); ++i) {
saturated_oil_table_[1][i] = 1.0/saturated_oil_table_[1][i];
}
undersat_oil_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
const auto &undersaturatedPvtoTable = pvtoTable.getUnderSaturatedTable(i);
undersat_oil_tables_[i][0] = undersaturatedPvtoTable.getColumn("P").vectorCopy();
undersat_oil_tables_[i][1] = undersaturatedPvtoTable.getColumn("BO").vectorCopy();
undersat_oil_tables_[i][2] = undersaturatedPvtoTable.getColumn("MU").vectorCopy();
// store the inverted oil formation factor
for (size_t j=0; j<undersat_oil_tables_[i][1].size(); ++j) {
undersat_oil_tables_[i][1][j] = 1.0/undersat_oil_tables_[i][1][j];
}
}
// Fill in additional entries in undersaturated tables by interpolating/extrapolating 1/Bo and mu_o ...
int iPrev = -1;
int iNext = 1;
while (undersat_oil_tables_[iNext][0].size() < 2) {
++iNext;
}
assert(iNext < sz);
for (int i=0; i<sz; ++i) {
if (undersat_oil_tables_[i][0].size() > 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<double>::iterator itPrev = upper_bound(undersat_oil_tables_[iPrev][0].begin(),
undersat_oil_tables_[iPrev][0].end(),pressure0+1.);
if (itPrev == undersat_oil_tables_[iPrev][0].end()) {
--itPrev; // Extrapolation ...
} else if (itPrev == undersat_oil_tables_[iPrev][0].begin()) {
++itPrev;
}
if (itPrev == undersat_oil_tables_[iPrev][0].end()-1) {
flagPrev = false; // Last data set for "prev" ...
}
double dPPrev = *itPrev - *(itPrev-1);
pressure = *itPrev;
int index = int(itPrev - undersat_oil_tables_[iPrev][0].begin());
slopePrevBinv = (undersat_oil_tables_[iPrev][1][index] - undersat_oil_tables_[iPrev][1][index-1])/dPPrev;
slopePrevVisc = (undersat_oil_tables_[iPrev][2][index] - undersat_oil_tables_[iPrev][2][index-1])/dPPrev;
}
if (flagNext) {
std::vector<double>::iterator itNext = upper_bound(undersat_oil_tables_[iNext][0].begin(),
undersat_oil_tables_[iNext][0].end(),pressure0+1.);
if (itNext == undersat_oil_tables_[iNext][0].end()) {
--itNext; // Extrapolation ...
} else if (itNext == undersat_oil_tables_[iNext][0].begin()) {
++itNext;
}
if (itNext == undersat_oil_tables_[iNext][0].end()-1) {
flagNext = false; // Last data set for "next" ...
}
double dPNext = *itNext - *(itNext-1);
if (flagPrev) {
pressure = std::min(pressure,*itNext);
} else {
pressure = *itNext;
}
int index = int(itNext - undersat_oil_tables_[iNext][0].begin());
slopeNextBinv = (undersat_oil_tables_[iNext][1][index] - undersat_oil_tables_[iNext][1][index-1])/dPNext;
slopeNextVisc = (undersat_oil_tables_[iNext][2][index] - undersat_oil_tables_[iNext][2][index-1])/dPNext;
}
double dP = pressure - pressure0;
if (iPrev >= 0) {
double w = (saturated_oil_table_[3][i] - saturated_oil_table_[3][iPrev]) /
(saturated_oil_table_[3][iNext] - saturated_oil_table_[3][iPrev]);
undersat_oil_tables_[i][0].push_back(pressure0+dP);
undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back() +
dP*(slopePrevBinv+w*(slopeNextBinv-slopePrevBinv)));
undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back() +
dP*(slopePrevVisc+w*(slopeNextVisc-slopePrevVisc)));
} else {
undersat_oil_tables_[i][0].push_back(pressure0+dP);
undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back()+dP*slopeNextBinv);
undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back()+dP*slopeNextVisc);
}
}
}
}
// Destructor
MiscibilityLiveOil::~MiscibilityLiveOil()
{
}
double MiscibilityLiveOil::getViscosity(int /*region*/, double press, const surfvol_t& surfvol) const
{
return miscible_oil(press, surfvol, 2, false);
}
void MiscibilityLiveOil::getViscosity(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const
{
assert(pressures.size() == surfvol.size());
int num = pressures.size();
output.resize(num);
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
output[i] = miscible_oil(pressures[i][phase], surfvol[i], 2, false);
}
}
// Dissolved gas-oil ratio
double MiscibilityLiveOil::R(int /*region*/, double press, const surfvol_t& surfvol) const
{
return evalR(press, surfvol);
}
void MiscibilityLiveOil::R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const
{
assert(pressures.size() == surfvol.size());
int num = pressures.size();
output.resize(num);
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
output[i] = evalR(pressures[i][phase], surfvol[i]);
}
}
// Dissolved gas-oil ratio derivative
double MiscibilityLiveOil::dRdp(int /*region*/, double press, const surfvol_t& surfvol) const
{
double R = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (R < maxR ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[3],
press);
} else {
return 0.0; // Undersaturated case
}
}
void MiscibilityLiveOil::dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_R,
std::vector<double>& output_dRdp) const
{
assert(pressures.size() == surfvol.size());
int num = pressures.size();
output_R.resize(num);
output_dRdp.resize(num);
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
evalRDeriv(pressures[i][phase], surfvol[i], output_R[i], output_dRdp[i]);
}
}
double MiscibilityLiveOil::B(int /*region*/, double press, const surfvol_t& surfvol) const
{
return evalB(press, surfvol);
}
void MiscibilityLiveOil::B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const
{
assert(pressures.size() == surfvol.size());
int num = pressures.size();
output.resize(num);
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
output[i] = evalB(pressures[i][phase], surfvol[i]);
}
}
double MiscibilityLiveOil::dBdp(int /*region*/, double press, const surfvol_t& surfvol) const
{
// if (surfvol[Liquid] == 0.0) return 0.0; // To handle no-oil case.
double Bo = evalB(press, surfvol); // \TODO check if we incur virtual call overhead here.
return -Bo*Bo*miscible_oil(press, surfvol, 1, true);
}
void MiscibilityLiveOil::dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_B,
std::vector<double>& output_dBdp) const
{
assert(pressures.size() == surfvol.size());
B(pressures, surfvol, phase, output_B);
int num = pressures.size();
output_dBdp.resize(num);
#pragma omp parallel for
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 MiscibilityLiveOil::evalR(double press, const surfvol_t& surfvol) const
{
if (surfvol[Vapour] == 0.0) {
return 0.0;
}
double R = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (R < maxR ) { // Saturated case
return R;
} else {
return maxR; // Undersaturated case
}
}
void MiscibilityLiveOil::evalRDeriv(const double press, const surfvol_t& surfvol,
double& R, double& dRdp) const
{
if (surfvol[Vapour] == 0.0) {
R = 0.0;
dRdp = 0.0;
return;
}
R = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (R < maxR ) {
// Saturated case
dRdp = linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[3],
press);
} else {
// Undersaturated case
R = maxR;
dRdp = 0.0;
}
}
double MiscibilityLiveOil::evalB(double press, const surfvol_t& surfvol) const
{
// if (surfvol[Liquid] == 0.0) return 1.0; // To handle no-oil case.
return 1.0/miscible_oil(press, surfvol, 1, false);
}
void MiscibilityLiveOil::evalBDeriv(const double press, const surfvol_t& surfvol,
double& B, double& dBdp) const
{
B = evalB(press, surfvol);
dBdp = -B*B*miscible_oil(press, surfvol, 1, true);
}
double MiscibilityLiveOil::miscible_oil(double press, const surfvol_t& surfvol,
int item, bool deriv) const
{
int section;
double R = linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[3],
press, section);
double maxR = (surfvol[Liquid] == 0.0) ? 0.0 : surfvol[Vapour]/surfvol[Liquid];
if (deriv) {
if (R < maxR ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[3], maxR);
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
} else {
if (R < maxR ) { // Saturated case
return linearInterpolation(saturated_oil_table_[0],
saturated_oil_table_[item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[3], maxR);
double w = (maxR - saturated_oil_table_[3][is]) /
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
assert(undersat_oil_tables_[is][0].size() >= 2);
assert(undersat_oil_tables_[is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

100
opm/porsol/blackoil/fluid/MiscibilityLiveOil.hpp
View File

@ -1,100 +0,0 @@
//===========================================================================
//
// File: MiscibilityLiveOil.hpp
//
// Created: Wed Feb 10 09:08:09 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#ifndef SINTEF_MISCIBILITYLIVEOIL_HEADER
#define SINTEF_MISCIBILITYLIVEOIL_HEADER
/** Class for miscible live oil.
* Detailed description.
*/
#include "MiscibilityProps.hpp"
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
namespace Opm
{
class PvtoTable;
class MiscibilityLiveOil : public MiscibilityProps
{
public:
MiscibilityLiveOil(const PvtoTable& pvtoTable);
virtual ~MiscibilityLiveOil();
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<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_B,
std::vector<double>& output_dBdp) const;
virtual void R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const;
virtual void dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_R,
std::vector<double>& output_dRdp) const;
protected:
double evalR(double press, const surfvol_t& surfvol) const;
void evalRDeriv(double press, const surfvol_t& surfvol, double& R, double& dRdp) const;
double evalB(double press, const surfvol_t& surfvol) const;
void evalBDeriv(double press, const surfvol_t& surfvol, double& B, double& dBdp) const;
// item: 1=B 2=mu;
double miscible_oil(double press, const surfvol_t& surfvol, int item,
bool deriv = false) const;
// PVT properties of live oil (with dissolved gas)
std::vector<std::vector<double> > saturated_oil_table_;
std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;
};
}
#endif // SINTEF_MISCIBILITYLIVEOIL_HEADER

53
opm/porsol/blackoil/fluid/MiscibilityProps.cpp
View File

@ -1,53 +0,0 @@
//===========================================================================
//
// File: MiscibilityProps.cpp
//
// Created: Wed Feb 10 09:05:05 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include "MiscibilityProps.hpp"
using namespace std;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityProps::MiscibilityProps()
{
}
MiscibilityProps::~MiscibilityProps()
{
}
} // namespace Opm

87
opm/porsol/blackoil/fluid/MiscibilityProps.hpp
View File

@ -1,87 +0,0 @@
//===========================================================================
//
// File: MiscibilityProps.hpp
//
// Created: Wed Feb 10 09:04:35 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// 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 <http://www.gnu.org/licenses/>.
*/
#ifndef SINTEF_MISCIBILITYPROPS_HEADER
#define SINTEF_MISCIBILITYPROPS_HEADER
/** Base class for properties of fluids and rocks.
* Detailed description.
*/
#include "BlackoilDefs.hpp"
#include <vector>
#include <array>
namespace Opm
{
class MiscibilityProps : public BlackoilDefs
{
public:
typedef CompVec surfvol_t;
MiscibilityProps();
virtual ~MiscibilityProps();
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const = 0;
virtual double B (int region, double press, const surfvol_t& surfvol) const = 0;
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const = 0;
virtual double R (int region, double press, const surfvol_t& surfvol) const = 0;
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const = 0;
virtual void getViscosity(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const = 0;
virtual void B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const = 0;
virtual void dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_B,
std::vector<double>& output_dBdp) const = 0;
virtual void R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output) const = 0;
virtual void dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvol,
int phase,
std::vector<double>& output_R,
std::vector<double>& output_dRdp) const = 0;
};
} // namespace Opm
#endif // SINTEF_MISCIBILITYPROPS_HEADER

210
opm/porsol/blackoil/fluid/MiscibilityWater.hpp
View File

@ -1,210 +0,0 @@
//===========================================================================
//
// File: MiscibilityWater.hpp
//
// Created: Tue May 18 10:26:13 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Bjørn Spjelkavik <bsp@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPENRS_MISCIBILITYWATER_HEADER
#define OPENRS_MISCIBILITYWATER_HEADER
#include "MiscibilityProps.hpp"
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
#include <opm/parser/eclipse/Deck/DeckItem.hpp>
// Forward declaration.
class PVTW;
namespace Opm
{
class MiscibilityWater : public MiscibilityProps
{
public:
typedef std::vector<std::vector<double> > table_t;
MiscibilityWater(const DeckKeyword& pvtwKeyword)
{
const auto& pvtwRecord = pvtwKeyword.getRecord(0);
ref_press_ = pvtwRecord.getItem("P_REF").getSIDouble(0);
ref_B_ = pvtwRecord.getItem("WATER_VOL_FACTOR").getSIDouble(0);
comp_ = pvtwRecord.getItem("WATER_COMPRESSIBILITY").getSIDouble(0);
viscosity_ = pvtwRecord.getItem("WATER_VISCOSITY").getSIDouble(0);
if (pvtwRecord.getItem("WATER_VISCOSIBILITY").getSIDouble(0) != 0.0) {
OPM_THROW(std::runtime_error, "MiscibilityWater does not support 'viscosibility'.");
}
}
MiscibilityWater(double visc)
: ref_press_(0.0),
ref_B_(1.0),
comp_(0.0),
viscosity_(visc)
{
}
// WTF?? we initialize a class for water from a keyword for oil?
void initFromPvcdo(const DeckKeyword& pvcdoKeyword)
{
const auto& pvcdoRecord = pvcdoKeyword.getRecord(0);
ref_press_ = pvcdoRecord.getItem("P_REF").getSIDouble(0);
ref_B_ = pvcdoRecord.getItem("OIL_VOL_FACTOR").getSIDouble(0);
comp_ = pvcdoRecord.getItem("OIL_COMPRESSIBILITY").getSIDouble(0);
viscosity_ = pvcdoRecord.getItem("OIL_VISCOSITY").getSIDouble(0);
if (pvcdoRecord.getItem("OIL_VISCOSIBILITY").getSIDouble(0) != 0.0) {
OPM_THROW(std::runtime_error, "MiscibilityWater does not support 'viscosibility'.");
}
}
virtual ~MiscibilityWater()
{
}
virtual double getViscosity(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return viscosity_;
}
virtual void getViscosity(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>&,
int,
std::vector<double>& output) const
{
int num = pressures.size();
output.clear();
output.resize(num, viscosity_);
}
virtual double B(int /*region*/, double press, const surfvol_t& /*surfvol*/) const
{
if (comp_) {
// Computing a polynomial approximation to the exponential.
double x = comp_*(press - ref_press_);
return ref_B_/(1.0 + x + 0.5*x*x);
} else {
return ref_B_;
}
}
virtual void B(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>&,
int phase,
std::vector<double>& output) const
{
int num = pressures.size();
if (comp_) {
output.resize(num);
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
// Computing a polynomial approximation to the exponential.
double x = comp_*(pressures[i][phase] - ref_press_);
output[i] = ref_B_/(1.0 + x + 0.5*x*x);
}
} else {
output.clear();
output.resize(num, ref_B_);
}
}
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const
{
if (comp_) {
return -comp_*B(region, press, surfvol);
} else {
return 0.0;
}
}
virtual void dBdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>& surfvols,
int phase,
std::vector<double>& output_B,
std::vector<double>& output_dBdp) const
{
B(pressures, surfvols, phase, output_B);
int num = pressures.size();
if (comp_) {
output_dBdp.resize(num);
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
output_dBdp[i] = -comp_*output_B[i];
}
} else {
output_dBdp.clear();
output_dBdp.resize(num, 0.0);
}
}
virtual double R(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 0.0;
}
virtual void R(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>&,
int,
std::vector<double>& output) const
{
int num = pressures.size();
output.clear();
output.resize(num, 0.0);
}
virtual double dRdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
{
return 0.0;
}
virtual void dRdp(const std::vector<PhaseVec>& pressures,
const std::vector<CompVec>&,
int,
std::vector<double>& output_R,
std::vector<double>& 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);
}
private:
double ref_press_;
double ref_B_;
double comp_;
double viscosity_;
};
}
#endif // OPENRS_MISCIBILITYWATER_HEADER

963
opm/porsol/mimetic/TpfaCompressible.hpp
View File

@ -1,963 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED
#define OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED
#include <opm/porsol/mimetic/TpfaCompressibleAssembler.hpp>
#include <opm/porsol/blackoil/BlackoilFluid.hpp>
#include <opm/core/linalg/LinearSolverIstl.hpp>
#include <opm/porsol/common/BoundaryConditions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/SparseTable.hpp>
#include <array>
#include <iostream>
#include <memory>
namespace Opm
{
template <class GridInterface,
class RockInterface,
class FluidInterface,
class WellsInterface,
class BCInterface>
class TpfaCompressible
{
public:
typedef TpfaCompressibleAssembler PressureAssembler;
/// @brief
/// Default constructor. Does nothing.
TpfaCompressible()
: pgrid_(0), prock_(0), pfluid_(0)
{
}
/// @brief
/// Initializes run-time parameters of the solver.
void init(const Opm::ParameterGroup& param)
{
// Initialize inflow mixture to a fixed, user-provided mix.
typename FluidInterface::CompVec mix(0.0);
const int nc = FluidInterface::numComponents;
double inflow_mixture_gas = param.getDefault("inflow_mixture_gas", 1.0);
double inflow_mixture_oil = param.getDefault("inflow_mixture_oil", 0.0);
switch (nc) {
case 2:
mix[FluidInterface::Gas] = inflow_mixture_gas;
mix[FluidInterface::Oil] = inflow_mixture_oil;
break;
case 3: {
double inflow_mixture_water = param.getDefault("inflow_mixture_water", 0.0);
mix[FluidInterface::Water] = inflow_mixture_water;
mix[FluidInterface::Gas] = inflow_mixture_gas;
mix[FluidInterface::Oil] = inflow_mixture_oil;
break;
}
default:
OPM_THROW(std::runtime_error, "Unhandled number of components: " << nc);
}
inflow_mixture_ = mix;
linsolver_.reset(new LinearSolverIstl(param));
flux_rel_tol_ = param.getDefault("flux_rel_tol", 1e-5);
press_rel_tol_ = param.getDefault("press_rel_tol", 1e-5);
max_num_iter_ = param.getDefault("max_num_iter", 15);
max_relative_voldiscr_ = param.getDefault("max_relative_voldiscr", 0.15);
relax_time_voldiscr_ = param.getDefault("relax_time_voldiscr", 0.0);
relax_weight_pressure_iteration_ = param.getDefault("relax_weight_pressure_iteration", 1.0);
experimental_jacobian_ = param.getDefault("experimental_jacobian", false);
nonlinear_residual_tolerance_ = param.getDefault("nonlinear_residual_tolerance", 0.0);
output_residual_ = param.getDefault("output_residual", false);
voldisc_factor_ = param.getDefault("voldisc_factor", 1.0);
}
/// @brief
/// Accessor for the inflow mixture.
typename FluidInterface::CompVec inflowMixture() const
{
return inflow_mixture_;
}
/// @brief
/// Setup routine, does grid/rock-dependent initialization.
///
/// @param [in] grid
/// The grid.
///
/// @param [in] rock
/// The cell-wise permeabilities and porosities.
///
/// @param [in] fluid
/// Fluid properties.
///
/// @param [in] wells
/// Well specifications.
///
/// @param [in] grav
/// Gravity vector. Its Euclidian two-norm value
/// represents the strength of the gravity field (in units
/// of m/s^2) while its direction is the direction of
/// gravity in the current model.
///
/// @param [in] bc
/// Boundary conditions.
///
/// @param [in] face_pressure
/// Face pressures. Allow optional micromanagement of
/// pressure boundary conditions.
///
void setup(const GridInterface& grid,
const RockInterface& rock,
const FluidInterface& fluid,
const WellsInterface& wells,
const typename GridInterface::Vector& grav,
const BCInterface& bc,
const std::vector<typename FluidInterface::PhaseVec>* face_pressure = 0)
{
pgrid_ = &grid;
prock_ = &rock;
pfluid_ = &fluid;
pwells_ = &wells;
gravity_ = grav;
// Extract perm tensors.
const double* perm = &(rock.permeability(0)(0,0));
poro_.clear();
poro_.resize(grid.numCells(), 1.0);
for (int i = 0; i < grid.numCells(); ++i) {
poro_[i] = rock.porosity(i);
}
// Initialize
psolver_.init(grid, wells, perm, &poro_[0], grav);
// Build bctypes_ and bcvalues_.
int num_faces = grid.numFaces();
bctypes_.clear();
bctypes_.resize(num_faces, PressureAssembler::FBC_UNSET);
bcvalues_.clear();
bcvalues_.resize(num_faces, 0.0);
for (int face = 0; face < num_faces; ++face) {
int bid = pgrid_->boundaryId(face);
if (bid == 0) {
bctypes_[face] = PressureAssembler::FBC_UNSET;
continue;
}
FlowBC face_bc = bc.flowCond(bid);
if (face_bc.isDirichlet()) {
bctypes_[face] = PressureAssembler::FBC_PRESSURE;
bcvalues_[face] = face_bc.pressure();
if (face_bc.pressure() < 0.0) {
bcvalues_[face] = (*face_pressure)[face][FluidInterface::Liquid];
}
} else if (face_bc.isNeumann()) {
bctypes_[face] = PressureAssembler::FBC_FLUX;
bcvalues_[face] = face_bc.outflux(); // TODO: may have to switch sign here depending on orientation.
if (bcvalues_[face] != 0.0) {
OPM_THROW(std::runtime_error, "Nonzero Neumann conditions not yet properly implemented "
"(signs must be fixed, also face pressures are not correctly computed for this case)");
}
} else {
OPM_THROW(std::runtime_error, "Unhandled boundary condition type.");
}
}
// Setup unchanging well data structures.
perf_wells_.clear();
perf_cells_.clear();
perf_A_.clear();
perf_mob_.clear();
perf_sat_.clear();
int num_wells = pwells_->numWells();
for (int well = 0; well < num_wells; ++well) {
int num_perf = pwells_->numPerforations(well);
for (int perf = 0; perf < num_perf; ++perf) {
int cell = pwells_->wellCell(well, perf);
perf_wells_.push_back(well);
perf_cells_.push_back(cell);
}
}
int num_perf = perf_wells_.size();
perf_A_.resize(num_perf*numPhases*numComponents);
perf_mob_.resize(num_perf*numPhases);
perf_sat_.resize(num_perf);
}
double volumeDiscrepancyLimit() const
{
return max_relative_voldiscr_;
}
const std::vector<double>& faceTransmissibilities()
{
return psolver_.faceTransmissibilities();
}
bool volumeDiscrepancyAcceptable(const std::vector<typename FluidInterface::PhaseVec>& cell_pressure,
const std::vector<typename FluidInterface::PhaseVec>& face_pressure,
const std::vector<double>& well_perf_pressure,
const std::vector<typename FluidInterface::CompVec>& cell_z,
const double dt)
{
computeFluidProps(cell_pressure, face_pressure, well_perf_pressure, cell_z, dt);
double rel_voldiscr = *std::max_element(fp_.relvoldiscr.begin(), fp_.relvoldiscr.end());
if (rel_voldiscr > max_relative_voldiscr_) {
std::cout << " Relative volume discrepancy too large: " << rel_voldiscr << std::endl;
return false;
} else {
std::cout << " Relative volume discrepancy ok: " << rel_voldiscr << std::endl;
return true;
}
}
enum ReturnCode { SolveOk, VolumeDiscrepancyTooLarge, FailedToConverge };
/// @brief
/// Construct and solve system of linear equations for the
/// phase pressure values on cells and faces, also compute
/// total face fluxes.
///
/// @param [inout] cell_pressure
/// Phase pressures per cell.
///
/// @param [inout] face_pressure
/// Phase pressures per face.
///
/// @param [inout] cell_z
/// Surface volume per cell. Only changed if the @code
/// transport @endcode argument is true.
///
/// @param [out] face_flux
/// Total (summed over all phases) volume flux (signed)
/// across each face.
///
/// @param [out] well_perf_pressures
/// Pressure in each well perforation.
///
/// @param [out] well_perf_fluxes
/// Total (summed over all phases) volume flux (signed,
/// positive meaning injection) from each well perforation.
///
/// @param [in] src
/// Explicit source terms. One scalar value for each grid
/// cell representing the rate (in units of m^3/s) of fluid
/// being injected into (>0) or extracted from (<0) a given
/// grid cell.
///
/// @param [in] dt
/// Timestep for pressure solver.
///
/// @param [in] transport
/// If true, modify @code cell_z @endcode by IMPES scheme.
///
ReturnCode solve(std::vector<typename FluidInterface::PhaseVec>& cell_pressure,
std::vector<typename FluidInterface::PhaseVec>& face_pressure,
const std::vector<typename FluidInterface::CompVec>& cell_z,
std::vector<double>& face_flux,
std::vector<double>& well_bhp_pressures,
std::vector<double>& well_perf_pressures,
std::vector<double>& well_perf_fluxes,
const std::vector<double>& src,
const double dt)
{
// Set up initial state.
// \TODO We are currently using Liquid phase pressure,
// what is correct with capillary pressure?
SolverState state;
int num_cells = cell_pressure.size();
int num_faces = face_pressure.size();
state.cell_pressure.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
state.cell_pressure[cell] = cell_pressure[cell][FluidInterface::Liquid];
}
state.face_pressure.resize(num_faces);
for (int face = 0; face < num_faces; ++face) {
state.face_pressure[face] = face_pressure[face][FluidInterface::Liquid];
}
state.face_flux.clear();
state.face_flux.resize(num_faces);
state.well_bhp_pressure = well_bhp_pressures;
state.well_perf_pressure = well_perf_pressures;
state.well_perf_flux = well_perf_fluxes;
// Run solver.
ReturnCode retcode;
if (nonlinear_residual_tolerance_ == 0.0) {
// Old version.
retcode = solveImpl(cell_z, src, dt, state);
} else {
// New version using residual reduction as
// convergence criterion.
retcode = solveImplNew(cell_z, src, dt, state);
}
// Copy results to output variables.
if (retcode == SolveOk) {
// \TODO Currently all phase pressures will be equal,
// what is correct with capillary pressure?
for (int cell = 0; cell < num_cells; ++cell) {
cell_pressure[cell] = state.cell_pressure[cell];
}
for (int face = 0; face < num_faces; ++face) {
face_pressure[face] = state.face_pressure[face];
}
face_flux = state.face_flux;
well_bhp_pressures = state.well_bhp_pressure;
well_perf_pressures = state.well_perf_pressure;
well_perf_fluxes = state.well_perf_flux;
}
return retcode;
}
/// Call this function after solve().
double stableStepIMPES()
{
return psolver_.explicitTimestepLimit(&fp_.face_data.state_matrix[0][0][0],
&fp_.face_data.mobility[0][0],
&fp_.face_data.mobility_deriv[0][0][0],
&(pfluid_->surfaceDensities()[0]));
}
/// Do an IMPES step using the facilities of the pressure solver.
void doStepIMPES(std::vector<typename FluidInterface::CompVec>& cell_z,
const double dt)
{
psolver_.explicitTransport(dt, &(cell_z[0][0]));
}
private:
const GridInterface* pgrid_;
const RockInterface* prock_;
const FluidInterface* pfluid_;
const WellsInterface* pwells_;
typename GridInterface::Vector gravity_;
// typename FluidInterface::FluidData fp_;
Opm::AllFluidData fp_;
std::vector<double> poro_;
PressureAssembler psolver_;
std::unique_ptr<LinearSolverIstl> linsolver_;
std::vector<PressureAssembler::FlowBCTypes> bctypes_;
std::vector<double> bcvalues_;
typename FluidInterface::CompVec inflow_mixture_;
double flux_rel_tol_;
double press_rel_tol_;
int max_num_iter_;
double max_relative_voldiscr_;
double relax_time_voldiscr_;
double relax_weight_pressure_iteration_;
bool experimental_jacobian_;
bool output_residual_;
double nonlinear_residual_tolerance_;
double voldisc_factor_;
typedef typename FluidInterface::PhaseVec PhaseVec;
typedef typename FluidInterface::CompVec CompVec;
enum { numPhases = FluidInterface::numPhases,
numComponents = FluidInterface::numComponents };
std::vector<int> perf_wells_;
std::vector<int> perf_cells_;
std::vector<double> perf_A_; // Flat storage.
std::vector<double> perf_mob_; // Flat storage.
std::vector<PhaseVec> perf_sat_;
std::vector<double> perf_gpot_; // Flat storage.
// Only intended to avoid extra allocations.
mutable std::vector<PhaseVec> helper_cell_pressure_;
mutable std::vector<PhaseVec> helper_face_pressure_;
struct SolverState
{
std::vector<double> cell_pressure;
std::vector<double> face_pressure;
std::vector<double> face_flux;
std::vector<double> well_bhp_pressure;
std::vector<double> well_perf_pressure;
std::vector<double> well_perf_flux;
};
void computeFluidProps(const std::vector<typename FluidInterface::PhaseVec>& phase_pressure,
const std::vector<typename FluidInterface::PhaseVec>& phase_pressure_face,
const std::vector<double>& well_perf_pressure,
const std::vector<typename FluidInterface::CompVec>& cell_z,
const double dt)
{
fp_.computeNew(*pgrid_, *prock_, *pfluid_, gravity_, phase_pressure, phase_pressure_face, cell_z, inflow_mixture_, dt);
// Properties at well perforations.
// \TODO only need to recompute this once per pressure update.
// No, that is false, at production perforations the cell z is
// used, which may change every step.
unsigned int perfcount = 0;
int num_wells = pwells_->numWells();
for (int well = 0; well < num_wells; ++well) {
bool inj = pwells_->type(well) == WellsInterface::Injector;
int num_perf = pwells_->numPerforations(well);
for (int perf = 0; perf < num_perf; ++perf) {
int cell = pwells_->wellCell(well, perf);
// \TODO handle capillary in perforation pressure below?
PhaseVec well_pressure = inj ? PhaseVec(well_perf_pressure[perfcount]) : phase_pressure[cell];
CompVec well_mixture = inj ? pwells_->injectionMixture(cell) : cell_z[cell];
typename FluidInterface::FluidState state = pfluid_->computeState(well_pressure, well_mixture);
std::copy(&state.phase_to_comp_[0][0], &state.phase_to_comp_[0][0] + numComponents*numPhases,
&perf_A_[perfcount*numPhases*numComponents]);
std::copy(state.mobility_.begin(), state.mobility_.end(),
&perf_mob_[perfcount*numPhases]);
perf_sat_[perfcount] = state.saturation_;
++perfcount;
}
}
assert(perfcount == perf_wells_.size());
}
// Convert scalar pressures to phase pressures, then call computeFluidProps().
void computeFluidPropsScalarPress(const std::vector<double>& cell_pressure_scalar,
const std::vector<double>& face_pressure_scalar,
const std::vector<double>& well_perf_pressure,
const std::vector<typename FluidInterface::CompVec>& cell_z,
const double dt)
{
// Copy to phase pressures. \TODO handle capillary pressure.
int num_cells = pgrid_->numCells();
int num_faces = pgrid_->numFaces();
helper_cell_pressure_.resize(num_cells);
helper_face_pressure_.resize(num_faces);
for (int cell = 0; cell < num_cells; ++cell) {
helper_cell_pressure_[cell] = cell_pressure_scalar[cell];
}
for (int face = 0; face < num_faces; ++face) {
helper_face_pressure_[face] = face_pressure_scalar[face];
}
computeFluidProps(helper_cell_pressure_, helper_face_pressure_, well_perf_pressure, cell_z, dt);
}
// Compute res = Ax - b.
void computeLinearResidual(const PressureAssembler::LinearSystem& s, std::vector<double>& res)
{
res.resize(s.n);
for (int row = 0; row < s.n; ++row) {
res[row] = -s.b[row];
for (int i = s.ia[row]; i < s.ia[row + 1]; ++i) {
res[row] += s.sa[i]*s.x[s.ja[i]];
}
}
}
// Compute residual and Jacobian of the new formulation.
void computeResidualJacobian(const std::vector<double>& initial_voldiscr,
const std::vector<double>& cell_pressure_scalar,
const std::vector<double>& cell_pressure_scalar_initial,
const std::vector<double>& well_bhp,
const std::vector<double>& src,
const double dt,
PressureAssembler::LinearSystem& linsys,
std::vector<double>& res)
{
std::vector<double> zero(initial_voldiscr.size(), 0.0);
// Assemble system matrix and rhs.
psolver_.assemble(&src[0],
&bctypes_[0],
&bcvalues_[0],
dt,
&fp_.cell_data.total_compressibility[0],
&zero[0],
&fp_.cell_data.state_matrix[0][0][0],
&fp_.face_data.state_matrix[0][0][0],
&perf_A_[0],
&fp_.face_data.mobility[0][0],
&perf_mob_[0],
&cell_pressure_scalar_initial[0],
&fp_.face_data.gravity_potential[0][0],
&perf_gpot_[0],
&(pfluid_->surfaceDensities()[0]));
psolver_.linearSystem(linsys);
// The linear system is for direct evaluation, we want a residual based approach.
// First we compute the residual for the original code.
int num_cells = pgrid_->numCells();
std::copy(cell_pressure_scalar.begin(), cell_pressure_scalar.end(), linsys.x);
std::copy(well_bhp.begin(), well_bhp.end(), linsys.x + num_cells);
computeLinearResidual(linsys, res);
// Then we compute the residual we actually want by subtracting terms that do not
// appear in the new formulation and adding the new terms.
for (int cell = 0; cell < num_cells; ++cell) {
double tc = fp_.cell_data.total_compressibility[cell];
double dres = tc*(cell_pressure_scalar[cell] - cell_pressure_scalar_initial[cell]);
dres -= 1.0 - fp_.cell_data.total_phase_volume_density[cell];
dres *= pgrid_->cellVolume(cell)*prock_->porosity(cell)/dt;
res[cell] -= dres;
}
// Change the jacobian by adding/subtracting the necessary terms.
for (int cell = 0; cell < num_cells; ++cell) {
for (int i = linsys.ia[cell]; i < linsys.ia[cell + 1]; ++i) {
if (linsys.ja[i] == cell) {
double tc = fp_.cell_data.total_compressibility[cell];
linsys.sa[i] -= tc*pgrid_->cellVolume(cell)*prock_->porosity(cell)/dt;
double ex = fp_.cell_data.experimental_term[cell];
linsys.sa[i] += ex*pgrid_->cellVolume(cell)*prock_->porosity(cell)/dt;
}
}
}
}
// Compute the well potentials. Assumes that the perforation variables
// have been set properly: perf_[wells_|cells_|A_].
void computeWellPotentials(std::vector<double>& wellperf_gpot) const
{
int num_perf = perf_cells_.size();
wellperf_gpot.resize(num_perf*numPhases);
for (int perf = 0; perf < num_perf; ++perf) {
int well = perf_wells_[perf];
int cell = perf_cells_[perf];
typename GridInterface::Vector pos = pgrid_->cellCentroid(cell);
// With wells, we assume that gravity is in the z-direction.
assert(gravity_[0] == 0.0 && gravity_[1] == 0.0);
double depth_delta = pos[2] - pwells_->referenceDepth(well);
double gh = gravity_[2]*depth_delta;
// At is already transposed since in Fortran order.
const double* At = &perf_A_[perf*numPhases*numComponents];
PhaseVec rho = pfluid_->phaseDensities(At);
for (int phase = 0; phase < numPhases; ++phase) {
// Gravity potential is (by phase) \rho_\alpha g h
wellperf_gpot[numPhases*perf + phase] = rho[phase]*gh;
}
}
}
// Compute the relative changes in fluxes and pressures.
static std::pair<double, double>
computeFluxPressChanges(const std::vector<double>& face_flux,
const std::vector<double>& well_perf_fluxes,
const std::vector<double>& cell_pressure_scalar,
const std::vector<double>& start_face_flux,
const std::vector<double>& start_perf_flux,
const std::vector<double>& start_cell_press)
{
int num_faces = face_flux.size();
int num_perf = well_perf_fluxes.size();
int num_cells = cell_pressure_scalar.size();
double max_flux_face = std::max(std::fabs(*std::min_element(face_flux.begin(), face_flux.end())),
std::fabs(*std::max_element(face_flux.begin(), face_flux.end())));
double max_flux_perf = num_perf == 0 ? 0.0
: std::max(std::fabs(*std::min_element(well_perf_fluxes.begin(), well_perf_fluxes.end())),
std::fabs(*std::max_element(well_perf_fluxes.begin(), well_perf_fluxes.end())));
double max_flux = std::max(max_flux_face, max_flux_perf);
double max_press = std::max(std::fabs(*std::min_element(cell_pressure_scalar.begin(),
cell_pressure_scalar.end())),
std::fabs(*std::max_element(cell_pressure_scalar.begin(),
cell_pressure_scalar.end())));
double flux_change_infnorm = 0.0;
double press_change_infnorm = 0.0;
for (int face = 0; face < num_faces; ++face) {
flux_change_infnorm = std::max(flux_change_infnorm,
std::fabs(face_flux[face] - start_face_flux[face]));
}
for (int perf = 0; perf < num_perf; ++perf) {
flux_change_infnorm = std::max(flux_change_infnorm,
std::fabs(well_perf_fluxes[perf] - start_perf_flux[perf]));
}
for (int cell = 0; cell < num_cells; ++cell) {
press_change_infnorm = std::max(press_change_infnorm,
std::fabs(cell_pressure_scalar[cell] - start_cell_press[cell]));
}
double flux_rel_difference = flux_change_infnorm/max_flux;
double press_rel_difference = press_change_infnorm/max_press;
return std::make_pair(flux_rel_difference, press_rel_difference);
}
// Compute well perforation pressures.
void computeWellPerfPressures(const std::vector<double>& well_perf_fluxes,
const std::vector<double>& well_bhp,
const std::vector<double>& wellperf_gpot,
std::vector<double>& well_perf_pressures) const
{
// Compute averaged saturations for each well. This code
// assumes that flow is either in or out of any single
// well, not both.
int num_perf = well_perf_fluxes.size();
std::vector<PhaseVec> well_sat(pwells_->numWells(), PhaseVec(0.0));
std::vector<double> well_flux(pwells_->numWells(), 0.0);
for (int perf = 0; perf < num_perf; ++perf) {
int well = perf_wells_[perf];
double flux = well_perf_fluxes[perf];
well_flux[well] += flux;
PhaseVec tmp = perf_sat_[perf];
tmp *= flux;
well_sat[well] += tmp;
}
for (int well = 0; well < pwells_->numWells(); ++well) {
well_sat[well] *= 1.0/well_flux[well];
}
// Compute well_perf_pressures
for (int perf = 0; perf < num_perf; ++perf) {
well_perf_pressures[perf] = well_bhp[perf_wells_[perf]];
PhaseVec sat = well_sat[perf_wells_[perf]];
for (int phase = 0; phase < numPhases; ++phase) {
well_perf_pressures[perf]
+= sat[phase]*wellperf_gpot[numPhases*perf + phase];
}
}
}
// Sets the initial volume discrepancy, applying relaxation if requested.
bool getInitialVolumeDiscrepancy(const double dt, std::vector<double>& voldiscr_initial)
{
voldiscr_initial = fp_.voldiscr;
double rel_voldiscr = *std::max_element(fp_.relvoldiscr.begin(), fp_.relvoldiscr.end());
if (rel_voldiscr > max_relative_voldiscr_) {
std::cout << "Initial relative volume discrepancy too large: " << rel_voldiscr << std::endl;
return false;
}
if (relax_time_voldiscr_ > 0.0) {
double relax = std::min(1.0,dt/relax_time_voldiscr_);
std::transform(voldiscr_initial.begin(), voldiscr_initial.end(), voldiscr_initial.begin(),
std::binder1st<std::multiplies<double> >(std::multiplies<double>() , relax));
}
if (voldisc_factor_ != 1.0) {
std::transform(voldiscr_initial.begin(), voldiscr_initial.end(), voldiscr_initial.begin(),
std::binder1st<std::multiplies<double> >(std::multiplies<double>(), voldisc_factor_));
}
return true;
}
// Modifies cell pressures, face pressures and face fluxes of
// state such that
// state <- weight*state + (1 - weight)*start_state
// So weight == 1.0 does nothing.
void relaxState(const double weight,
const SolverState& start_state,
SolverState& state)
{
int num_cells = pgrid_->numCells();
int num_faces = pgrid_->numFaces();
for (int cell = 0; cell < num_cells; ++cell) {
state.cell_pressure[cell] = weight*state.cell_pressure[cell] + (1.0-weight)*start_state.cell_pressure[cell];
}
for (int face = 0; face < num_faces; ++face) {
state.face_pressure[face] = weight*state.face_pressure[face] + (1.0-weight)*start_state.face_pressure[face];
state.face_flux[face] = weight*state.face_flux[face] + (1.0-weight)*start_state.face_flux[face];
}
}
// Implements the main nonlinear loop of the pressure solver.
ReturnCode solveImpl(const std::vector<typename FluidInterface::CompVec>& cell_z,
const std::vector<double>& src,
const double dt,
SolverState& state)
{
int num_cells = pgrid_->numCells();
std::vector<double> cell_pressure_initial = state.cell_pressure;
std::vector<double> voldiscr_initial;
SolverState start_state;
// ------------ Main iteration loop -------------
for (int iter = 0; iter < max_num_iter_; ++iter) {
start_state = state;
// (Re-)compute fluid properties.
computeFluidPropsScalarPress(state.cell_pressure, state.face_pressure, state.well_perf_pressure, cell_z, dt);
// Initialization for the first iteration only.
if (iter == 0) {
bool voldisc_ok = getInitialVolumeDiscrepancy(dt, voldiscr_initial);
if (!voldisc_ok) {
return VolumeDiscrepancyTooLarge;
}
// perf_gpot_ is computed once per pressure solve,
// while perf_A_, perf_mob_ are recoomputed
// for every iteration.
computeWellPotentials(perf_gpot_);
}
if (experimental_jacobian_) {
// Compute residual and jacobian.
PressureAssembler::LinearSystem s;
std::vector<double> residual;
computeResidualJacobian(voldiscr_initial, state.cell_pressure, cell_pressure_initial,
state.well_bhp_pressure, src, dt, s, residual);
if (output_residual_) {
// Temporary hack to get output of residual.
static int psolve_iter = -1;
if (iter == 0) {
++psolve_iter;
}
std::ostringstream oss;
oss << "residual-" << psolve_iter << '-' << iter << ".dat";
std::ofstream outres(oss.str().c_str());
std::copy(residual.begin(), residual.end(), std::ostream_iterator<double>(outres, "\n"));
}
// Solve system for dp, that is, we use residual as the rhs.
LinearSolverIstl::LinearSolverReport result
= linsolver_->solve(s.n, s.nnz, s.ia, s.ja, s.sa, &residual[0], s.x);
if (!result.converged) {
OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
<< "Residual reduction achieved is " << result.residual_reduction << '\n');
}
// Set x so that the call to computePressuresAndFluxes() will work.
// Recall that x now contains dp, and we want it to contain p - dp
for (int cell = 0; cell < num_cells; ++cell) {
s.x[cell] = state.cell_pressure[cell] - s.x[cell];
}
for (int well = 0; well < pwells_->numWells(); ++well) {
s.x[num_cells + well] = state.well_bhp_pressure[well] - s.x[num_cells + well];
}
} else {
// Assemble system matrix and rhs.
psolver_.assemble(&src[0],
&bctypes_[0],
&bcvalues_[0],
dt,
&fp_.cell_data.total_compressibility[0],
&voldiscr_initial[0],
&fp_.cell_data.state_matrix[0][0][0],
&fp_.face_data.state_matrix[0][0][0],
&perf_A_[0],
&fp_.face_data.mobility[0][0],
&perf_mob_[0],
&cell_pressure_initial[0],
&fp_.face_data.gravity_potential[0][0],
&perf_gpot_[0],
&(pfluid_->surfaceDensities()[0]));
PressureAssembler::LinearSystem s;
psolver_.linearSystem(s);
// Solve system.
LinearSolverIstl::LinearSolverReport res = linsolver_->solve(s.n, s.nnz, s.ia, s.ja, s.sa, s.b, s.x);
if (!res.converged) {
OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << res.iterations << " iterations.\n"
<< "Residual reduction achieved is " << res.residual_reduction << '\n');
}
}
// Get pressures and face fluxes.
psolver_.computePressuresAndFluxes(state.cell_pressure, state.face_pressure, state.face_flux,
state.well_bhp_pressure, state.well_perf_flux);
// Relaxation
if (relax_weight_pressure_iteration_ != 1.0) {
relaxState(relax_weight_pressure_iteration_, start_state, state);
}
// Compute state.well_perf_pressure.
computeWellPerfPressures(state.well_perf_flux, state.well_bhp_pressure,
perf_gpot_, state.well_perf_pressure);
// Compute relative changes for pressure and flux.
std::pair<double, double> rel_changes
= computeFluxPressChanges(state.face_flux, state.well_perf_flux, state.cell_pressure,
start_state.face_flux, start_state.well_perf_flux, start_state.cell_pressure);
double flux_rel_difference = rel_changes.first;
double press_rel_difference = rel_changes.second;
// Test for convergence.
if (iter == 0) {
std::cout << "Iteration Rel. flux change Rel. pressure change\n";
}
std::cout.precision(5);
std::cout << std::setw(6) << iter
<< std::setw(24) << flux_rel_difference
<< std::setw(24) << press_rel_difference << std::endl;
std::cout.precision(16);
if (flux_rel_difference < flux_rel_tol_ || press_rel_difference < press_rel_tol_) {
std::cout << "Pressure solver converged. Number of iterations: " << iter + 1 << '\n' << std::endl;
return SolveOk;
}
}
return FailedToConverge;
}
// Implements the main nonlinear loop of the pressure solver.
ReturnCode solveImplNew(const std::vector<typename FluidInterface::CompVec>& cell_z,
const std::vector<double>& src,
const double dt,
SolverState& state)
{
int num_cells = pgrid_->numCells();
int num_faces = pgrid_->numFaces();
std::vector<double> cell_pressure_initial = state.cell_pressure;
std::vector<double> voldiscr_initial;
SolverState start_state;
// ------------ Main iteration loop -------------
for (int iter = 0; iter < max_num_iter_; ++iter) {
start_state = state;
// (Re-)compute fluid properties.
computeFluidPropsScalarPress(state.cell_pressure, state.face_pressure, state.well_perf_pressure, cell_z, dt);
// Initialization for the first iteration only.
if (iter == 0) {
bool voldisc_ok = getInitialVolumeDiscrepancy(dt, voldiscr_initial);
if (!voldisc_ok) {
return VolumeDiscrepancyTooLarge;
}
// perf_gpot_ is computed once per pressure solve,
// while perf_A_, perf_mob_ are recoomputed
// for every iteration.
computeWellPotentials(perf_gpot_);
}
// Compute residual and jacobian.
PressureAssembler::LinearSystem s;
std::vector<double> residual;
computeResidualJacobian(voldiscr_initial, state.cell_pressure, cell_pressure_initial,
state.well_bhp_pressure, src, dt, s, residual);
if (output_residual_) {
// Temporary hack to get output of residual.
static int psolve_iter = -1;
if (iter == 0) {
++psolve_iter;
}
std::ostringstream oss;
oss << "residual-" << psolve_iter << '-' << iter << ".dat";
std::ofstream outres(oss.str().c_str());
std::copy(residual.begin(), residual.end(), std::ostream_iterator<double>(outres, "\n"));
}
// Find the maxnorm of the residual.
double maxres = std::max(std::fabs(*std::max_element(residual.begin(), residual.end())),
std::fabs(*std::min_element(residual.begin(), residual.end())));
// Solve system for dp, that is, we use residual as the rhs.
LinearSolverIstl::LinearSolverReport result
= linsolver_->solve(s.n, s.nnz, s.ia, s.ja, s.sa, &residual[0], s.x);
if (!result.converged) {
OPM_THROW(std::runtime_error, "Linear solver failed to converge in " << result.iterations << " iterations.\n"
<< "Residual reduction achieved is " << result.residual_reduction << '\n');
}
// Set x so that the call to computePressuresAndFluxes() will work.
// Recall that x now contains dp, and we want it to contain p - dp
for (int cell = 0; cell < num_cells; ++cell) {
s.x[cell] = state.cell_pressure[cell] - s.x[cell];
}
for (int well = 0; well < pwells_->numWells(); ++well) {
s.x[num_cells + well] = state.well_bhp_pressure[well] - s.x[num_cells + well];
}
// Get pressures and face fluxes.
psolver_.computePressuresAndFluxes(state.cell_pressure, state.face_pressure, state.face_flux,
state.well_bhp_pressure, state.well_perf_flux);
// Relaxation
if (relax_weight_pressure_iteration_ != 1.0) {
double ww = relax_weight_pressure_iteration_;
for (int cell = 0; cell < num_cells; ++cell) {
state.cell_pressure[cell] = ww*state.cell_pressure[cell] + (1.0-ww)*start_state.cell_pressure[cell];
}
if (iter > 0) {
for (int face = 0; face < num_faces; ++face) {
state.face_pressure[face] = ww*state.face_pressure[face] + (1.0-ww)*start_state.face_pressure[face];
state.face_flux[face] = ww*state.face_flux[face] + (1.0-ww)*start_state.face_flux[face];
}
}
}
// Compute state.well_perf_pressure.
computeWellPerfPressures(state.well_perf_flux, state.well_bhp_pressure,
perf_gpot_, state.well_perf_pressure);
// Compute relative changes for pressure and flux.
std::pair<double, double> rel_changes
= computeFluxPressChanges(state.face_flux, state.well_perf_flux, state.cell_pressure,
start_state.face_flux, start_state.well_perf_flux, start_state.cell_pressure);
double flux_rel_difference = rel_changes.first;
double press_rel_difference = rel_changes.second;
// Test for convergence.
if (iter == 0) {
std::cout << "Iteration Rel. flux change Rel. pressure change Residual\n";
}
std::cout.precision(5);
std::cout << std::setw(6) << iter
<< std::setw(24) << flux_rel_difference
<< std::setw(24) << press_rel_difference
<< std::setw(24) << maxres << std::endl;
std::cout.precision(16);
if (maxres < nonlinear_residual_tolerance_) {
std::cout << "Pressure solver converged. Number of iterations: " << iter + 1 << '\n' << std::endl;
return SolveOk;
}
}
return FailedToConverge;
}
};
} // namespace Opm
#endif // OPM_TPFACOMPRESSIBLE_HEADER_INCLUDED

527
opm/porsol/mimetic/TpfaCompressibleAssembler.hpp
View File

@ -1,527 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TPFACOMPRESSIBLEASSEMBLER_HEADER_INCLUDED
#define OPM_TPFACOMPRESSIBLEASSEMBLER_HEADER_INCLUDED
#include <opm/core/pressure/tpfa/cfs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/pressure/legacy_well.h>
#include <opm/core/pressure/tpfa/compr_quant.h>
#include <dune/grid/common/GridAdapter.hpp>
#include <stdexcept>
/// @brief
/// Encapsulates the cfs_tpfa (= compressible flow solver
/// two-point flux approximation) solver modules.
class TpfaCompressibleAssembler
{
public:
/// @brief
/// Default constructor, does nothing.
TpfaCompressibleAssembler()
: state_(Uninitialized), data_(0), bc_(0)
{
wells_.number_of_wells = 0;
}
/// @brief
/// Destructor.
~TpfaCompressibleAssembler()
{
flow_conditions_destroy(bc_);
cfs_tpfa_destroy(data_);
}
/// @brief
/// Initialize the solver's structures for a given grid, for well setup also call initWells().
/// @tparam Grid This must conform to the SimpleGrid concept.
/// @tparam Wells This must conform to the SimpleWells concept.
/// @param grid The grid object.
/// @param wells Well specifications.
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
/// @param perm Porosity by cell.
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
template <class Grid, class Wells>
void init(const Grid& grid, const Wells& wells, const double* perm, const double* porosity,
const typename Grid::Vector& gravity)
{
initWells(wells);
init(grid, perm, porosity, gravity);
}
/// @brief
/// Initialize the solver's structures for a given grid, for well setup also call initWells().
/// @tparam Grid This must conform to the SimpleGrid concept.
/// @param grid The grid object.
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
template <class Grid>
void init(const Grid& grid, const double* perm, const double* porosity, const typename Grid::Vector& gravity)
{
// Build C grid structure.
grid_.init(grid);
// Initialize data.
int num_phases = 3;
well_t* w = 0;
if (wells_.number_of_wells != 0) {
w = &wells_;
}
data_ = cfs_tpfa_construct(grid_.c_grid(), w, num_phases);
if (!data_) {
throw std::runtime_error("Failed to initialize cfs_tpfa solver.");
}
// Compute half-transmissibilities
int num_cells = grid.numCells();
int ngconn = grid_.c_grid()->cell_facepos[num_cells];
ncf_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
int num_local_faces = grid.numCellFaces(cell);
ncf_[cell] = num_local_faces;
}
htrans_.resize(ngconn);
tpfa_htrans_compute(grid_.c_grid(), perm, &htrans_[0]);
// Compute transmissibilities.
trans_.resize(grid_.numFaces());
tpfa_trans_compute(grid_.c_grid(), &htrans_[0], &trans_[0]);
// Compute pore volumes.
porevol_.resize(num_cells);
for (int i = 0; i < num_cells; ++i) {
porevol_[i] = porosity[i]*grid.cellVolume(i);
}
// Set gravity.
if (Grid::dimension != 3) {
throw std::logic_error("Only 3 dimensions supported currently.");
}
std::copy(gravity.begin(), gravity.end(), gravity_);
state_ = Initialized;
}
/// Boundary condition types.
enum FlowBCTypes { FBC_UNSET = BC_NOFLOW ,
FBC_PRESSURE = BC_PRESSURE ,
FBC_FLUX = BC_FLUX_TOTVOL };
/// @brief
/// Assemble the sparse system.
/// You must call init() prior to calling assemble().
/// @param sources Source terms, one per cell. Positive numbers
/// are sources, negative are sinks.
/// @param total_mobilities Scalar total mobilities, one per cell.
/// @param omegas Gravity term, one per cell. In a multi-phase
/// flow setting this is equal to
/// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f]
/// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a
/// phase density and \f$\lambda_t\f$ is the total mobility.
void assemble(const double* sources,
const FlowBCTypes* bctypes,
const double* bcvalues,
const double dt,
const double* totcompr,
const double* voldiscr,
const double* cellA, // num phases^2 * num cells, fortran ordering!
const double* faceA, // num phases^2 * num faces, fortran ordering!
const double* wellperfA,
const double* phasemobf,
const double* phasemobwellperf,
const double* cell_pressure,
const double* gravcapf,
const double* wellperf_gpot,
const double* /* surf_dens */)
{
if (state_ == Uninitialized) {
throw std::runtime_error("Error in TpfaCompressibleAssembler::assemble(): You must call init() prior to calling assemble().");
}
UnstructuredGrid* g = grid_.c_grid();
// Boundary conditions.
gather_boundary_conditions(bctypes, bcvalues);
// Source terms from user.
double* src = const_cast<double*>(sources); // Ugly? Yes. Safe? I think so.
// Wells.
well_t* wells = NULL;
well_control_t* wctrl = NULL;
struct completion_data* wcompl = NULL;
if (wells_.number_of_wells != 0) {
wells = &wells_;
wctrl = &wctrl_;
wcompl = &wcompl_;
// The next objects already have the correct sizes.
std::copy(wellperf_gpot, wellperf_gpot + well_gpot_storage_.size(), well_gpot_storage_.begin());
std::copy(wellperfA, wellperfA + well_A_storage_.size(), well_A_storage_.begin());
std::copy(phasemobwellperf, phasemobwellperf + well_phasemob_storage_.size(), well_phasemob_storage_.begin());
}
// Assemble the embedded linear system.
compr_quantities cq = { 3 ,
const_cast<double *>(totcompr ) ,
const_cast<double *>(voldiscr ) ,
const_cast<double *>(cellA ) ,
const_cast<double *>(faceA ) ,
const_cast<double *>(phasemobf) };
// Call the assembly routine. After this, linearSystem() may be called.
cfs_tpfa_assemble(g, dt, wells, bc_, src,
&cq, &trans_[0], gravcapf,
wctrl, wcompl,
cell_pressure, &porevol_[0],
data_);
phasemobf_.assign(phasemobf, phasemobf + grid_.numFaces()*3);
gravcapf_.assign(gravcapf, gravcapf + grid_.numFaces()*3);
state_ = Assembled;
}
/// Encapsulate a sparse linear system in CSR format.
struct LinearSystem
{
int n;
int nnz;
int* ia;
int* ja;
double* sa;
double* b;
double* x;
};
/// @brief
/// Access the linear system assembled.
/// You must call assemble() prior to calling linearSystem().
/// @param[out] s The linear system encapsulation to modify.
/// After this call, s will point to linear system structures
/// that are owned and allocated internally.
void linearSystem(LinearSystem& s)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TpfaCompressibleAssembler::linearSystem(): "
"You must call assemble() prior to calling linearSystem().");
}
s.n = data_->A->m;
s.nnz = data_->A->nnz;
s.ia = data_->A->ia;
s.ja = data_->A->ja;
s.sa = data_->A->sa;
s.b = data_->b;
s.x = data_->x;
}
/// @brief
/// Compute cell pressures and face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// computePressuresAndFluxes().
/// @param[out] cell_pressures Cell pressure values.
/// @param[out] face_areas Face flux values.
void computePressuresAndFluxes(std::vector<double>& cell_pressures,
std::vector<double>& face_pressures,
std::vector<double>& face_fluxes,
std::vector<double>& well_pressures,
std::vector<double>& well_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TpfaCompressibleAssembler::computePressuresAndFluxes(): "
"You must call assemble() (and solve the linear system) "
"prior to calling computePressuresAndFluxes().");
}
int num_cells = grid_.c_grid()->number_of_cells;
int num_faces = grid_.c_grid()->number_of_faces;
cell_pressures.clear();
cell_pressures.resize(num_cells, 0.0);
face_pressures.clear();
face_pressures.resize(num_faces, 0.0);
face_fluxes.clear();
face_fluxes.resize(num_faces, 0.0);
int np = 3; // Number of phases.
// Wells.
well_t* wells = NULL;
struct completion_data* wcompl = NULL;
double* wpress = 0;
double* wflux = 0;
if (wells_.number_of_wells != 0) {
wells = &wells_;
wcompl = &wcompl_;
well_pressures.resize(wells_.number_of_wells);
well_fluxes.resize(well_cells_storage_.size());
wpress = &well_pressures[0];
wflux = &well_fluxes[0];
}
cfs_tpfa_press_flux(grid_.c_grid(),
bc_, wells,
np, &trans_[0], &phasemobf_[0], &gravcapf_[0],
wcompl,
data_, &cell_pressures[0], &face_fluxes[0],
wpress, wflux);
cfs_tpfa_fpress(grid_.c_grid(), bc_, np, &htrans_[0],
&phasemobf_[0], &gravcapf_[0], data_, &cell_pressures[0],
&face_fluxes[0], &face_pressures[0]);
}
/// @brief
/// Explicit IMPES time step limit.
double explicitTimestepLimit(const double* faceA, // num phases^2 * num faces, fortran ordering!
const double* phasemobf,
const double* phasemobf_deriv,
const double* surf_dens)
{
compr_quantities cq = { 3, // nphases
0, // totcompr
0, // voldiscr
0, // Ac
const_cast<double *>(faceA) ,
const_cast<double *>(phasemobf) };
return cfs_tpfa_impes_maxtime(grid_.c_grid(), &cq, &trans_[0], &porevol_[0], data_,
phasemobf_deriv, surf_dens, gravity_);
}
/// @brief
/// Explicit IMPES transport.
void explicitTransport(const double dt,
double* cell_surfvols)
{
int np = 3; // Number of phases.
well_t* wells = NULL;
if (wells_.number_of_wells != 0) {
wells = &wells_;
}
cfs_tpfa_expl_mass_transport(grid_.c_grid(), wells, &wcompl_,
np, dt, &porevol_[0],
data_, cell_surfvols);
}
/// @brief
/// Compute cell fluxes from face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// faceFluxToCellFlux().
/// @param face_fluxes
/// @param face_areas Face flux values (usually output from computePressuresAndFluxes()).
/// @param[out] cell_fluxes Cell-wise flux values.
/// They are given in cell order, and for each cell there is
/// one value for each adjacent face (in the same order as the
/// cell-face topology of the grid). Positive values represent
/// fluxes out of the cell.
void faceFluxToCellFlux(const std::vector<double>& face_fluxes,
std::vector<double>& cell_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TpfaCompressibleAssembler::faceFluxToCellFlux(): "
"You must call assemble() (and solve the linear system) "
"prior to calling faceFluxToCellFlux().");
}
const UnstructuredGrid& g = *(grid_.c_grid());
int num_cells = g.number_of_cells;
cell_fluxes.resize(g.cell_facepos[num_cells]);
for (int cell = 0; cell < num_cells; ++cell) {
for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
int face = g.cell_faces[hface];
bool pos = (g.face_cells[2*face] == cell);
cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
}
}
}
/// @brief
/// Access the number of connections (faces) per cell. Deprecated, will be removed.
const std::vector<int>& numCellFaces() const
{
return ncf_;
}
const std::vector<double>& faceTransmissibilities() const
{
return trans_;
}
private:
// Disabling copy and assigment for now.
TpfaCompressibleAssembler(const TpfaCompressibleAssembler&);
TpfaCompressibleAssembler& operator=(const TpfaCompressibleAssembler&);
enum State { Uninitialized, Initialized, Assembled };
State state_;
// Solver data.
cfs_tpfa_data* data_;
// Grid.
GridAdapter grid_;
// Number of faces per cell.
std::vector<int> ncf_;
// Transmissibility storage.
std::vector<double> htrans_;
std::vector<double> trans_;
// Pore volumes.
std::vector<double> porevol_;
// Phase mobilities per face.
std::vector<double> phasemobf_;
// Gravity and capillary contributions (per face).
std::vector<double> gravcapf_;
// Gravity
double gravity_[3];
// Boundary conditions.
FlowBoundaryConditions *bc_;
// Well data
well_t wells_;
std::vector<int> well_connpos_storage_;
std::vector<int> well_cells_storage_;
well_control_t wctrl_;
std::vector<well_type> wctrl_type_storage_;
std::vector<well_control> wctrl_ctrl_storage_;
std::vector<double> wctrl_target_storage_;
struct completion_data wcompl_;
std::vector<double> well_prodind_storage_;
std::vector<double> well_gpot_storage_;
std::vector<double> well_A_storage_;
std::vector<double> well_phasemob_storage_;
/// @brief
/// Initialize wells in solver structure.
/// @tparam Wells
/// This must conform to the SimpleWells concept.
/// @param w
/// The well object.
template <class Wells>
void initWells(const Wells& w)
{
int num_wells = w.numWells();
if (num_wells == 0) {
wells_.number_of_wells = 0;
return;
}
wctrl_type_storage_.resize(num_wells);
wctrl_ctrl_storage_.resize(num_wells);
wctrl_target_storage_.resize(num_wells);
for (int i = 0; i < num_wells; ++i) {
wctrl_type_storage_[i] = (w.type(i) == Wells::Injector) ? INJECTOR : PRODUCER;
wctrl_ctrl_storage_[i] = (w.control(i) == Wells::Rate) ? RATE : BHP;
wctrl_target_storage_[i] = w.target(i);
int num_perf = w.numPerforations(i);
well_connpos_storage_.push_back(well_cells_storage_.size());
for (int j = 0; j < num_perf; ++j) {
well_cells_storage_.push_back(w.wellCell(i, j));
well_prodind_storage_.push_back(w.wellIndex(i, j));
}
}
well_connpos_storage_.push_back(well_cells_storage_.size());
int tot_num_perf = well_prodind_storage_.size();
well_gpot_storage_.resize(tot_num_perf*3);
well_A_storage_.resize(3*3*tot_num_perf);
well_phasemob_storage_.resize(3*tot_num_perf);
// Setup 'wells_'
wells_.number_of_wells = num_wells;
wells_.well_connpos = &well_connpos_storage_[0];
wells_.well_cells = &well_cells_storage_[0];
// Setup 'wctrl_'
wctrl_.type = &wctrl_type_storage_[0];
wctrl_.ctrl = &wctrl_ctrl_storage_[0];
wctrl_.target = &wctrl_target_storage_[0];
// Setup 'wcompl_'
wcompl_.WI = &well_prodind_storage_[0];
wcompl_.gpot = &well_gpot_storage_[0];
wcompl_.A = &well_A_storage_[0];
wcompl_.phasemob = &well_phasemob_storage_[0];
}
void
gather_boundary_conditions(const FlowBCTypes* bctypes ,
const double* bcvalues)
{
if (bc_ == 0) {
bc_ = flow_conditions_construct(0);
}
else {
flow_conditions_clear(bc_);
}
int ok = bc_ != 0;
for (std::size_t i = 0, nf = grid_.numFaces(); ok && (i < nf); ++i) {
if (bctypes[ i ] == FBC_PRESSURE) {
ok = flow_conditions_append(BC_PRESSURE,
static_cast<int>(i),
bcvalues[ i ],
bc_);
}
else if (bctypes[ i ] == FBC_FLUX) {
ok = flow_conditions_append(BC_FLUX_TOTVOL,
static_cast<int>(i),
bcvalues[ i ],
bc_);
}
}
if (! ok) {
flow_conditions_destroy(bc_);
bc_ = 0;
}
}
}; // class TpfaCompressibleAssembler
#endif // OPM_TPFACOMPRESSIBLEASSEMBLER_HEADER_INCLUDED