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

Added gravity transport solver functionality. Not tested.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-08-29 13:49:02 +02:00
parent becf485795
commit eb52ee58c1
2 changed files with 76 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
opm/polymer/TransportModelCompressiblePolymer.cpp
View File

@ -201,11 +201,11 @@ namespace Opm
fractionalflow_(grid.number_of_cells, -1.0),
mc_(grid.number_of_cells, -1.0)
{
const int np = props.numPhases();
const int num_cells = grid.number_of_cells;
if (props.numPhases() != 2) {
THROW("Property object must have 2 phases");
}
int np = props.numPhases();
int num_cells = props.numCells();
visc_.resize(np*num_cells);
A_.resize(np*np*num_cells);
A0_.resize(np*np*num_cells);
@ -215,7 +215,7 @@ namespace Opm
for (int i = 0; i < num_cells; ++i) {
allcells_[i] = i;
}
props.satRange(props.numCells(), &allcells_[0], &smin_[0], &smax_[0]);
props.satRange(num_cells, &allcells_[0], &smin_[0], &smax_[0]);
}
@ -251,9 +251,9 @@ namespace Opm
res_counts.clear();
#endif
props_.viscosity(props_.numCells(), &pressure[0], NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(props_.numCells(), &pressure0[0], NULL, &allcells_[0], &A0_[0], NULL);
props_.matrix(props_.numCells(), &pressure[0], NULL, &allcells_[0], &A_[0], NULL);
props_.viscosity(grid_.number_of_cells, &pressure[0], NULL, &allcells_[0], &visc_[0], NULL);
props_.matrix(grid_.number_of_cells, &pressure0[0], NULL, &allcells_[0], &A0_[0], NULL);
props_.matrix(grid_.number_of_cells, &pressure[0], NULL, &allcells_[0], &A_[0], NULL);
computePorosity(grid_, porosity_standard_, rock_comp_, pressure0, porosity0_);
computePorosity(grid_, porosity_standard_, rock_comp_, pressure, porosity_);
@ -1124,7 +1124,8 @@ namespace Opm
double dmob_ds[4];
double dmob_dc[2];
double dmobwat_dc;
polyprops_.effectiveMobilitiesBoth(c, cmax, &visc_[cell], relperm, drelperm_ds,
const int np = props_.numPhases();
polyprops_.effectiveMobilitiesBoth(c, cmax, &visc_[np*cell], relperm, drelperm_ds,
mob, dmob_ds, dmobwat_dc, if_with_der);
ff = mob[0]/(mob[0] + mob[1]);
@ -1192,7 +1193,7 @@ namespace Opm
c0(tm.concentration_[cell]),
cmax0(tm.cmax0_[cell]),
porosity(tm.porosity_[cell]),
dtpv(tm.dt_/tm.porevolume_[cell]),
dtpv(tm.dt_/(tm.grid_.cell_volumes[cell]*porosity)),
dps(tm.polyprops_.deadPoreVol()),
rhor(tm.polyprops_.rockDensity())
@ -1283,30 +1284,52 @@ namespace Opm
{
double sat[2] = { s, 1.0 - s };
double relperm[2];
const int np = props_.numPhases();
props_.relperm(1, sat, &cell, relperm, 0);
polyprops_.effectiveMobilities(c, cmax0_[cell], visc_, relperm, mob);
polyprops_.effectiveMobilities(c, cmax0_[cell], &visc_[np*cell], relperm, mob);
}
void TransportModelCompressiblePolymer::initGravity(const double* grav)
{
// Set up gravflux_ = T_ij g (rho_w - rho_o) (z_i - z_j)
// Set up transmissibilities.
std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
const int nf = grid_.number_of_faces;
const int dim = grid_.dimensions;
trans_.resize(nf);
gravflux_.resize(nf);
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &gravflux_[0]);
const double delta_rho = props_.density()[0] - props_.density()[1];
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &trans_[0]);
// Remember gravity vector.
gravity_ = grav;
}
void TransportModelCompressiblePolymer::initGravityDynamic()
{
// Set up gravflux_ = T_ij g [ (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f)
// + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ]
// But b_w,i * rho_w,S = rho_w,i, which we compute with a call to props_.density().
// We assume that we already have stored T_ij in trans_.
// We also assume that the A_ matrices are updated from an earlier call to solve().
const int nc = grid_.number_of_cells;
const int nf = grid_.number_of_faces;
const int np = props_.numPhases();
ASSERT(np == 2);
const int dim = grid_.dimensions;
density_.resize(nc*np);
props_.density(grid_.number_of_cells, &A_[0], &density_[0]);
std::fill(gravflux_.begin(), gravflux_.end(), 0.0);
for (int f = 0; f < nf; ++f) {
const int* c = &grid_.face_cells[2*f];
double gdz = 0.0;
const double signs[2] = { 1.0, -1.0 };
if (c[0] != -1 && c[1] != -1) {
for (int d = 0; d < dim; ++d) {
gdz += grav[d]*(grid_.cell_centroids[dim*c[0] + d] - grid_.cell_centroids[dim*c[1] + d]);
for (int ci = 0; ci < 2; ++ci) {
double gdz = 0.0;
for (int d = 0; d < dim; ++d) {
gdz += gravity_[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]);
}
gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]);
}
}
gravflux_[f] *= delta_rho*gdz;
}
}
@ -1400,14 +1423,17 @@ namespace Opm
void TransportModelCompressiblePolymer::solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation,
std::vector<double>& concentration,
std::vector<double>& cmax)
{
// Assume that solve() has already been called, so that A_ and
// porosity_ are current.
initGravityDynamic();
// initialize variables.
porevolume_ = porevolume;
dt_ = dt;
toWaterSat(saturation, saturation_);
concentration_ = &concentration[0];

68
opm/polymer/TransportModelCompressiblePolymer.hpp
View File

@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTMODELPOLYMER_HEADER_INCLUDED
#define OPM_TRANSPORTMODELPOLYMER_HEADER_INCLUDED
#ifndef OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#define OPM_TRANSPORTMODELCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/polymer/PolymerProperties.hpp>
@ -87,6 +87,10 @@ namespace Opm
std::vector<double>& concentration,
std::vector<double>& cmax);
/// Initialise quantities needed by gravity solver.
/// \param[in] grav Gravity vector
void initGravity(const double* grav);
/// Solve for gravity segregation.
/// This uses a column-wise nonlinear Gauss-Seidel approach.
/// It assumes that the input columns contain cells in a single
@ -99,49 +103,20 @@ namespace Opm
/// \param[in, out] concentration Polymer concentration.
/// \param[in, out] cmax Highest concentration that has occured in a given cell.
void solveGravity(const std::vector<std::vector<int> >& columns,
const double* porevolume,
const double dt,
std::vector<double>& saturation,
std::vector<double>& concentration,
std::vector<double>& cmax);
public: // But should be made private...
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
void solveSingleCellBracketing(int cell);
void solveSingleCellNewton(int cell);
void solveSingleCellGradient(int cell);
void solveSingleCellNewtonSimple(int cell,bool use_sc);
class ResidualEquation;
void initGravity(const double* grav);
void solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux);
int solveGravityColumn(const std::vector<int>& cells);
// This should be private:
class ResidualEquation;
friend class TransportModelCompressiblePolymer::ResidualEquation;
void scToc(const double* x, double* x_c) const;
//
#ifdef PROFILING
class Newton_Iter {
public:
bool res_s;
int cell;
double s;
double c;
private:
Newton_Iter(bool res_s_val, int cell_val, double s_val, double c_val) {
res_s = res_s_val;
cell = cell_val;
s = s_val;
c = c_val;
}
};
std::list<Newton_Iter> res_counts;
#endif
private:
const UnstructuredGrid& grid_;
const double* porosity_standard_;
const BlackoilPropertiesInterface& props_;
@ -162,7 +137,7 @@ namespace Opm
double* cmax_;
std::vector<double> fractionalflow_; // one per cell
std::vector<double> mc_; // one per cell
std::vector<double> visc_;
std::vector<double> visc_; // viscosity (without polymer, for given pressure)
std::vector<double> A_;
std::vector<double> A0_;
std::vector<double> porosity0_;
@ -171,6 +146,9 @@ namespace Opm
std::vector<double> smax_;
// For gravity segregation.
const double* gravity_;
std::vector<double> trans_;
std::vector<double> density_;
std::vector<double> gravflux_;
std::vector<double> mob_;
std::vector<double> cmax0_;
@ -183,6 +161,7 @@ namespace Opm
std::vector<int> ja_upw_;
std::vector<int> ia_downw_;
std::vector<int> ja_downw_;
struct ResidualC;
struct ResidualS;
@ -190,6 +169,19 @@ namespace Opm
class ResidualCGrav;
class ResidualSGrav;
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
void solveSingleCellBracketing(int cell);
void solveSingleCellNewton(int cell);
void solveSingleCellGradient(int cell);
void solveSingleCellNewtonSimple(int cell,bool use_sc);
void solveSingleCellGravity(const std::vector<int>& cells,
const int pos,
const double* gravflux);
int solveGravityColumn(const std::vector<int>& cells);
void initGravityDynamic();
void fracFlow(double s, double c, double cmax, int cell, double& ff) const;
void fracFlowWithDer(double s, double c, double cmax, int cell, double& ff,