Copied and renamed basic building blocks for black oil pvt.

This commit is contained in:
Atgeirr Flø Rasmussen 2012-01-04 15:44:53 +01:00
parent dd4c0c2b1b
commit f8ad8b8bf0
10 changed files with 1475 additions and 0 deletions

View File

@ -0,0 +1,101 @@
/*
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 <tr1/array>
#include <iostream>
#include <boost/static_assert.hpp>
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 };
// We need types with operator= and constructor taking scalars
// for the small vectors and matrices, to save us from having to
// rewrite a large amount of code.
template <typename T, int N>
class SmallVec
{
public:
SmallVec()
{}
SmallVec(const T& elem)
{ assign(elem); }
SmallVec& operator=(const T& elem)
{ assign(elem); return *this; }
const T& operator[](int i) const
{ return data_[i]; }
T& operator[](int i)
{ return data_[i]; }
template <typename U>
void assign(const U& elem)
{
for (int i = 0; i < N; ++i) {
data_[i] = elem;
}
}
private:
T data_[N];
};
template <typename T, int Rows, int Cols>
class SmallMat
{
public:
SmallMat()
{}
SmallMat(const T& elem)
{ data_.assign(elem); }
SmallMat& operator=(const T& elem)
{ data_.assign(elem); return *this; }
typedef SmallVec<T, Cols> RowType;
const RowType& operator[](int i) const
{ return data_[i]; }
RowType& operator[](int i)
{ return data_[i]; }
private:
SmallVec<RowType, Rows> data_;
};
typedef double Scalar;
typedef SmallVec<Scalar, numComponents> CompVec;
typedef SmallVec<Scalar, numPhases> PhaseVec;
BOOST_STATIC_ASSERT(int(numComponents) == int(numPhases));
typedef SmallMat<Scalar, numComponents, numPhases> PhaseToCompMatrix;
typedef SmallMat<Scalar, numPhases, numPhases> PhaseJacobian;
// Attempting to guard against alignment issues.
BOOST_STATIC_ASSERT(sizeof(CompVec) == numComponents*sizeof(Scalar));
BOOST_STATIC_ASSERT(sizeof(PhaseVec) == numPhases*sizeof(Scalar));
BOOST_STATIC_ASSERT(sizeof(PhaseToCompMatrix) == numComponents*numPhases*sizeof(Scalar));
BOOST_STATIC_ASSERT(sizeof(PhaseJacobian) == numPhases*numPhases*sizeof(Scalar));
};
} // namespace Opm
#endif // OPM_BLACKOILDEFS_HEADER_INCLUDED

View File

@ -0,0 +1,97 @@
//===========================================================================
//
// 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"
namespace Opm
{
class MiscibilityLiveOil : public MiscibilityProps
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
MiscibilityLiveOil(const table_t& pvto);
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

View File

@ -0,0 +1,179 @@
//===========================================================================
//
// 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 <opm/core/fluid/blackoil/MiscibilityDead.hpp>
#include <algorithm>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <boost/lexical_cast.hpp>
#include <string>
#include <fstream>
using namespace std;
using namespace Dune;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityDead::MiscibilityDead(const table_t& pvd_table)
{
const int region_number = 0;
if (pvd_table.size() != 1) {
THROW("More than one PVT-region");
}
// Copy data
const int sz = pvd_table[region_number][0].size();
std::vector<double> press(sz);
std::vector<double> B_inv(sz);
std::vector<double> visc(sz);
for (int i = 0; i < sz; ++i) {
press[i] = pvd_table[region_number][0][i];
B_inv[i] = 1.0 / pvd_table[region_number][1][i];
visc[i] = pvd_table[region_number][2][i];
}
int samples = 1025;
buildUniformMonotoneTable(press, B_inv, samples, one_over_B_);
buildUniformMonotoneTable(press, visc, samples, viscosity_);
// Dumping the created tables.
// static int count = 0;
// std::ofstream os((std::string("dump-") + boost::lexical_cast<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);
}
}

View File

@ -0,0 +1,89 @@
//===========================================================================
//
// 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 <opm/core/fluid/blackoil/MiscibilityProps.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
namespace Opm
{
class MiscibilityDead : public MiscibilityProps
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
MiscibilityDead(const table_t& pvd_table);
virtual ~MiscibilityDead();
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
virtual double B(int region, double press, const surfvol_t& surfvol) const;
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
virtual double R(int region, double press, const surfvol_t& surfvol) const;
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
virtual void getViscosity(const std::vector<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
Dune::utils::UniformTableLinear<double> one_over_B_;
Dune::utils::UniformTableLinear<double> viscosity_;
};
}
#endif // SINTEF_MISCIBILITYDEAD_HEADER

View File

@ -0,0 +1,52 @@
//===========================================================================
//
// 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 "MiscibilityProps.hpp"
using namespace std;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityProps::MiscibilityProps()
{
}
MiscibilityProps::~MiscibilityProps()
{
}
} // namespace Opm

View File

@ -0,0 +1,87 @@
//===========================================================================
//
// 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 <tr1/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

View File

@ -0,0 +1,286 @@
//===========================================================================
//
// 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 <opm/core/fluid/blackoil/MiscibilityLiveGas.hpp>
#include <algorithm>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
using namespace std;
using namespace Dune;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityLiveGas::MiscibilityLiveGas(const table_t& pvtg)
{
// GAS, PVTG
const int region_number = 0;
if (pvtg.size() != 1) {
THROW("More than one PVD-region");
}
saturated_gas_table_.resize(4);
const int sz = pvtg[region_number].size();
for (int k=0; k<4; ++k) {
saturated_gas_table_[k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_gas_table_[0][i] = pvtg[region_number][i][0]; // p
saturated_gas_table_[1][i] = pvtg[region_number][i][2]; // Bg
saturated_gas_table_[2][i] = pvtg[region_number][i][3]; // mu_g
saturated_gas_table_[3][i] = pvtg[region_number][i][1]; // Rv
}
undersat_gas_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
undersat_gas_tables_[i].resize(3);
int tsize = (pvtg[region_number][i].size() - 1)/3;
undersat_gas_tables_[i][0].resize(tsize);
undersat_gas_tables_[i][1].resize(tsize);
undersat_gas_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) {
undersat_gas_tables_[i][0][j] = pvtg[region_number][i][++k]; // Rv
undersat_gas_tables_[i][1][j] = pvtg[region_number][i][++k]; // Bg
undersat_gas_tables_[i][2][j] = pvtg[region_number][i][++k]; // mu_g
}
}
}
// 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 = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (R < maxR ) { // Saturated case
return R;
} else {
return maxR; // Undersaturated case
}
}
void MiscibilityLiveGas::R(const std::vector<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 = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (R < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_gas_table_[0],
saturated_gas_table_[3],
press);
} else {
return 0.0; // Undersaturated case
}
}
void MiscibilityLiveGas::dRdp(const std::vector<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 = linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[3], press,
section);
double maxR = surfvol[Liquid]/surfvol[Vapour];
if (deriv) {
if (R < maxR ) { // Saturated case
return linearInterpolDerivative(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = section;
if (undersat_gas_tables_[is][0].size() < 2) {
double val = (saturated_gas_table_[item][is+1]
- saturated_gas_table_[item][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
return val;
}
double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = (val2 - val1)/
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
return val;
}
} else {
if (R < maxR ) { // Saturated case
return linearInterpolationExtrap(saturated_gas_table_[0],
saturated_gas_table_[item],
press);
} else { // Undersaturated case
int is = section;
// Extrapolate from first table section
if (is == 0 && press < saturated_gas_table_[0][0]) {
return linearInterpolationExtrap(undersat_gas_tables_[0][0],
undersat_gas_tables_[0][item],
maxR);
}
// Extrapolate from last table section
int ltp = saturated_gas_table_[0].size() - 1;
if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
return linearInterpolationExtrap(undersat_gas_tables_[ltp][0],
undersat_gas_tables_[ltp][item],
maxR);
}
// Interpolate between table sections
double w = (press - saturated_gas_table_[0][is]) /
(saturated_gas_table_[0][is+1] -
saturated_gas_table_[0][is]);
if (undersat_gas_tables_[is][0].size() < 2) {
double val = saturated_gas_table_[item][is] +
w*(saturated_gas_table_[item][is+1] -
saturated_gas_table_[item][is]);
return val;
}
double val1 =
linearInterpolationExtrap(undersat_gas_tables_[is][0],
undersat_gas_tables_[is][item],
maxR);
double val2 =
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
undersat_gas_tables_[is+1][item],
maxR);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

View File

@ -0,0 +1,92 @@
//===========================================================================
//
// 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"
namespace Opm
{
class MiscibilityLiveGas : public MiscibilityProps
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
MiscibilityLiveGas(const table_t& pvto);
virtual ~MiscibilityLiveGas();
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
virtual double R(int region, double press, const surfvol_t& surfvol) const;
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
virtual double B(int region, double press, const surfvol_t& surfvol) const;
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
virtual void getViscosity(const std::vector<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

View File

@ -0,0 +1,395 @@
//===========================================================================
//
// 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 <opm/core/fluid/blackoil/MiscibilityLiveOil.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <algorithm>
using namespace std;
using namespace Dune;
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
MiscibilityLiveOil::MiscibilityLiveOil(const table_t& pvto)
{
// OIL, PVTO
const int region_number = 0;
if (pvto.size() != 1) {
THROW("More than one PVD-region");
}
saturated_oil_table_.resize(4);
const int sz = pvto[region_number].size();
for (int k=0; k<4; ++k) {
saturated_oil_table_[k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_oil_table_[0][i] = pvto[region_number][i][1]; // p
saturated_oil_table_[1][i] = 1.0/pvto[region_number][i][2]; // 1/Bo
saturated_oil_table_[2][i] = pvto[region_number][i][3]; // mu_o
saturated_oil_table_[3][i] = pvto[region_number][i][0]; // Rs
}
undersat_oil_tables_.resize(sz);
for (int i=0; i<sz; ++i) {
undersat_oil_tables_[i].resize(3);
int tsize = (pvto[region_number][i].size() - 1)/3;
undersat_oil_tables_[i][0].resize(tsize);
undersat_oil_tables_[i][1].resize(tsize);
undersat_oil_tables_[i][2].resize(tsize);
for (int j=0, k=0; j<tsize; ++j) {
undersat_oil_tables_[i][0][j] = pvto[region_number][i][++k]; // p
undersat_oil_tables_[i][1][j] = 1.0/pvto[region_number][i][++k]; // 1/Bo
undersat_oil_tables_[i][2][j] = pvto[region_number][i][++k]; // mu_o
}
}
// 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 = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (R < maxR ) { // Saturated case
return linearInterpolDerivative(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 = linearInterpolationExtrap(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 = linearInterpolationExtrap(saturated_oil_table_[0],
saturated_oil_table_[3], press);
double maxR = surfvol[Vapour]/surfvol[Liquid];
if (R < maxR ) {
// Saturated case
dRdp = linearInterpolDerivative(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 = linearInterpolationExtrap(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 linearInterpolDerivative(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 =
linearInterpolDerivative(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolDerivative(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 linearInterpolationExtrap(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 =
linearInterpolationExtrap(undersat_oil_tables_[is][0],
undersat_oil_tables_[is][item],
press);
double val2 =
linearInterpolationExtrap(undersat_oil_tables_[is+1][0],
undersat_oil_tables_[is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

View File

@ -0,0 +1,97 @@
//===========================================================================
//
// 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"
namespace Opm
{
class MiscibilityLiveOil : public MiscibilityProps
{
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
MiscibilityLiveOil(const table_t& pvto);
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