OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-26 03:00:17 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added adsorbtion in residual computation.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-03-27 11:20:20 +02:00
parent 6d9f0ea3ec
commit acedd4b596
4 changed files with 110 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
examples/polymer_reorder.cpp
View File

@ -197,6 +197,28 @@ public:
return props_.density()[phase];
}
template <class PolyC,
class CAds,
class DCAdsDc>
void adsorbtion(const PolyC& c, const PolyC& cmax, CAds& cads, DCAdsDc& dcadsdc) {
cads = polyprops_.adsorbtionWithDer(c, cmax, &dcadsdc);
}
const std::vector<double> porosity() const {
const int num_cells = props_.numCells();
std::vector<double> porosity(num_cells, 0.);
const double* poro = props_.porosity();
for (std::vector<double>::iterator it = porosity.begin(); it != porosity.end(); ++it, ++poro) {
*it = poro[0];
}
return porosity;
}
double rockdensity() const {
return polyprops_.rockDensity();
}
template <class Sat,
class PolyC,
class Mob,
@ -364,14 +386,14 @@ public:
std::vector<double>& faceflux () { return flux_ ; }
std::vector<double>& saturation () { return sat_ ; }
std::vector<double>& concentration() { return concentration_; }
std::vector<double>& cmax() { return cmax_; }
std::vector<double>& maxconcentration() { return cmax_; }
const std::vector<double>& pressure () const { return press_ ; }
const std::vector<double>& facepressure() const { return fpress_; }
const std::vector<double>& faceflux () const { return flux_ ; }
const std::vector<double>& saturation () const { return sat_ ; }
const std::vector<double>& concentration() const { return concentration_; }
const std::vector<double>& cmax() const { return cmax_; }
const std::vector<double>& maxconcentration() const { return cmax_; }
private:
std::vector<double> press_ ;
@ -401,7 +423,7 @@ static void outputState(const UnstructuredGrid& grid,
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
dm["cmax"] = &state.cmax();
dm["cmax"] = &state.maxconcentration();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
@ -905,7 +927,7 @@ main(int argc, char** argv)
if (use_reorder) {
Opm::toWaterSat(state.saturation(), reorder_sat);
reorder_model.solve(&state.faceflux()[0], &reorder_src[0], stepsize, inflow_c,
&reorder_sat[0], &state.concentration()[0], &state.cmax()[0]);
&reorder_sat[0], &state.concentration()[0], &state.maxconcentration()[0]);
Opm::toBothSat(reorder_sat, state.saturation());
Opm::computeInjectedProduced(*props, state.saturation(), src, simtimer.currentStepLength(), injected, produced);
if (use_segregation_split) {
@ -915,7 +937,8 @@ main(int argc, char** argv)
// reorder_model.solveGravity(columns, simtimer.currentStepLength(), reorder_sat);
// Opm::toBothSat(reorder_sat, state.saturation());
} else {
colsolver.solve(columns, simtimer.currentStepLength(), state.saturation(), state.concentration());
colsolver.solve(columns, simtimer.currentStepLength(), state.saturation(), state.concentration(),
state.maxconcentration());
}
} else {
// // Not implemented for polymer
@ -941,7 +964,7 @@ main(int argc, char** argv)
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), polyprop.deadPoreVol());
polymass_adsorbed = Opm::computePolymerAdsorbed(*props, polyprop, porevol, state.cmax());
polymass_adsorbed = Opm::computePolymerAdsorbed(*props, polyprop, porevol, state.maxconcentration());
Opm::computeInjectedProduced(*props, polyprop, state.saturation(), state.concentration(),
src, simtimer.currentStepLength(), inflow_c,
injected, produced, polyinj, polyprod);

5
opm/polymer/GravityColumnSolverPolymer.hpp
View File

@ -48,12 +48,15 @@ namespace Opm
void solve(const std::map<int, std::vector<int> >& columns,
const double dt,
std::vector<double>& s,
std::vector<double>& c);
std::vector<double>& c,
std::vector<double>& cmax);
private:
void solveSingleColumn(const std::vector<int>& column_cells,
const double dt,
std::vector<double>& s,
std::vector<double>& c,
std::vector<double>& cmax,
std::vector<double>& sol_vec
);
Model& model_;

36
opm/polymer/GravityColumnSolverPolymer_impl.hpp
View File

@ -59,16 +59,20 @@ namespace Opm
};
struct StateWithZeroFlux
{
StateWithZeroFlux(std::vector<double>& s, std::vector<double>& c) : sat(s), cpoly(c) {}
StateWithZeroFlux(std::vector<double>& s, std::vector<double>& c, std::vector<double>& cmax_arg) : sat(s), cpoly(c), cmax(cmax_arg) {}
const ZeroVec& faceflux() const { return zv; }
const std::vector<double>& saturation() const { return sat; }
std::vector<double>& saturation() { return sat; }
const std::vector<double>& concentration() const { return cpoly; }
std::vector<double>& concentration() { return cpoly; }
const std::vector<double>& maxconcentration() const { return cmax; }
std::vector<double>& maxconcentration() { return cmax; }
ZeroVec zv;
std::vector<double>& sat;
std::vector<double>& cpoly;
std::vector<double>& cmax;
};
struct Vecs
{
Vecs(int sz) : sol(sz, 0.0) {}
@ -84,13 +88,13 @@ namespace Opm
Vecs v;
typedef std::vector<double> vector_type;
};
struct BandMatrixCoeff
struct BandMatrixCoeff
{
BandMatrixCoeff(int N, int ku, int kl) : N_(N), ku_(ku), kl_(kl), nrow_(2*kl + ku + 1) {
}
// compute the position where to store the coefficient of a matrix A_{i,j} (i,j=0,...,N-1)
// in a array which is sent to the band matrix solver of LAPACK.
int operator ()(int i, int j) const {
@ -116,11 +120,12 @@ namespace Opm
void GravityColumnSolverPolymer<Model>::solve(const std::map<int, std::vector<int> >& columns,
const double dt,
std::vector<double>& s,
std::vector<double>& c
std::vector<double>& c,
std::vector<double>& cmax
)
{
// Initialize model. These things are done for the whole grid!
StateWithZeroFlux state(s, c); // This holds s and c by reference.
StateWithZeroFlux state(s, c, cmax); // This holds s, c and cmax by reference.
JacSys sys(2*grid_.number_of_cells);
std::vector<double> increment(2*grid_.number_of_cells, 0.0);
model_.initStep(state, grid_, sys);
@ -131,7 +136,7 @@ namespace Opm
model_.initIteration(state, grid_, sys);
std::map<int, std::vector<int> >::const_iterator it;
for (it = columns.begin(); it != columns.end(); ++it) {
solveSingleColumn(it->second, dt, s, c, increment);
solveSingleColumn(it->second, dt, s, c, cmax, increment);
}
for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
sys.vector().writableSolution()[2*cell + 0] += increment[2*cell + 0];
@ -150,9 +155,9 @@ namespace Opm
THROW("Failed to converge!");
}
// Finalize.
// model_.finishIteration(); // Doesn't do anything in th 2p model.
// model_.finishIteration(); //
// finishStep() writes to state, which holds s by reference.
// This will update the entire grid's state...
// This will update the entire grid's state... cmax is updated here.
model_.finishStep(grid_, sys.vector().solution(), state);
}
@ -167,12 +172,13 @@ namespace Opm
const double dt,
std::vector<double>& s,
std::vector<double>& c,
std::vector<double>& cmax,
std::vector<double>& sol_vec)
{
// This is written only to work with SinglePointUpwindTwoPhase,
// not with arbitrary problem models.
int col_size = column_cells.size();
StateWithZeroFlux state(s, c); // This holds s by reference.
StateWithZeroFlux state(s, c, cmax); // This holds s by reference.
// Assemble.
const int kl = 3;
@ -185,7 +191,7 @@ namespace Opm
for (int ci = 0; ci < col_size; ++ci) {
std::vector<double> F(2, 0.);
std::vector<double> F(2, 0.);
std::vector<double> dFd1(4, 0.);
std::vector<double> dFd2(4, 0.);
std::vector<double> dF(4, 0.);
@ -198,7 +204,7 @@ namespace Opm
const int c1 = grid_.face_cells[2*face + 0];
const int c2 = grid_.face_cells[2*face + 1];
if (c1 == prev_cell || c2 == prev_cell || c1 == next_cell || c2 == next_cell) {
F.assign(2, 0.);
F.assign(2, 0.);
dFd1.assign(4, 0.);
dFd2.assign(4, 0.);
model_.fluxConnection(state, grid_, dt, cell, face, &F[0], &dFd1[0], &dFd2[0]);
@ -218,7 +224,7 @@ namespace Opm
hm[bmc(2*ci + 0, 2*ci + 1)] += dFd1[1];
hm[bmc(2*ci + 1, 2*ci + 0)] += dFd1[2];
hm[bmc(2*ci + 1, 2*ci + 1)] += dFd1[3];
rhs[2*ci + 0] += F[0];
rhs[2*ci + 1] += F[1];
}
@ -235,7 +241,7 @@ namespace Opm
rhs[2*ci + 1] += F[1];
}
// model_.sourceTerms(); // Not needed
// model_.sourceTerms(); // Not needed
// Solve.
const int num_rhs = 1;
int info = 0;

76
opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp
View File

@ -48,9 +48,9 @@ namespace Opm {
class ModelParameterStorage {
public:
ModelParameterStorage(int nc, int totconn)
: drho_(0.0), mob_(0), dmobds_(0), dmobwatdc_(0), mc_(0), dmcdc_(0),
porevol_(0), deadporespace_(0), dg_(0), sw_(0), c_(0), ds_(0), dsc_(0), pc_(0), dpc_(0), trans_(0),
data_()
: drho_(0.0), deadporespace_(0.0), rockdensity_(0.0), mob_(0), dmobds_(0), dmobwatdc_(0), mc_(0),
dmcdc_(0), porevol_(0), porosity_(0), dg_(0), sw_(0), c_(0),
ds_(0), dsc_(0), dcads_(0), dcadsdc_(0), pc_(0), dpc_(0), trans_(0), data_()
{
size_t alloc_sz;
@ -60,12 +60,15 @@ namespace Opm {
alloc_sz += nc; // mc_
alloc_sz += nc; // dmcdc_
alloc_sz += 1 * nc; // porevol_
alloc_sz += 1 * nc; // deadporespave__
alloc_sz += 1 * nc; // porosity_
alloc_sz += 1 * nc; // deadporespace_
alloc_sz += 1 * totconn; // dg_
alloc_sz += 1 * nc; // sw_
alloc_sz += 1 * nc; // c_
alloc_sz += 1 * nc; // dc_
alloc_sz += 1 * nc; // dsc_
alloc_sz += 1 * nc; // dcads_
alloc_sz += 1 * nc; // dcadsdc_
alloc_sz += 1 * nc; // pc_
alloc_sz += 1 * nc; // dpc_
alloc_sz += 1 * totconn; // trans_
@ -77,13 +80,15 @@ namespace Opm {
mc_ = dmobwatdc_ + (1 * nc );
dmcdc_ = mc_ + (1 * nc );
porevol_ = dmcdc_ + (1 * nc );
deadporespace_ = porevol_ + (1 * nc );
dg_ = deadporespace_ + (1 * nc );
porosity_ = porevol_ + (1 * nc );
dg_ = porosity_ + (1 * nc );
sw_ = dg_ + (1 * totconn);
c_ = sw_ + (1 * nc );
ds_ = c_ + (1 * nc );
dsc_ = ds_ + (1 * nc );
pc_ = dsc_ + (1 * nc );
dcads_ = dsc_ + (1 * nc );
dcadsdc_ = dcads_ + (1 * nc );
pc_ = dcadsdc_ + (1 * nc );
dpc_ = pc_ + (1 * nc );
trans_ = dpc_ + (1 * nc );
}
@ -91,6 +96,12 @@ namespace Opm {
double& drho () { return drho_ ; }
double drho () const { return drho_ ; }
double& deadporespace() { return deadporespace_ ; }
double deadporespace() const { return deadporespace_ ; }
double& rockdensity() { return rockdensity_ ; }
double rockdensity() const { return rockdensity_ ; }
double* mob (int cell) { return mob_ + (2*cell + 0); }
const double* mob (int cell) const { return mob_ + (2*cell + 0); }
@ -109,8 +120,9 @@ namespace Opm {
double* porevol() { return porevol_ ; }
double porevol(int cell) const { return porevol_[cell] ; }
double* deadporespace() { return deadporespace_ ; }
double deadporespace(int cell) const { return deadporespace_[cell] ; }
double* porosity() { return porosity_ ; }
double porosity(int cell) const { return porosity_[cell] ; }
double& dg(int i) { return dg_[i] ; }
double dg(int i) const { return dg_[i] ; }
@ -127,6 +139,12 @@ namespace Opm {
double& dsc(int cell) { return dsc_[cell] ; }
double dsc(int cell) const { return dsc_[cell] ; }
double& dcads(int cell) { return dcads_[cell] ; }
double dcads(int cell) const { return dcads_[cell] ; }
double& dcadsdc(int cell) { return dcadsdc_[cell] ; }
double dcadsdc(int cell) const { return dcadsdc_[cell] ; }
double& pc(int cell) { return pc_[cell] ; }
double pc(int cell) const { return pc_[cell] ; }
@ -138,18 +156,22 @@ namespace Opm {
private:
double drho_ ;
double deadporespace_ ;
double rockdensity_ ;
double *mob_ ;
double *dmobds_ ;
double *dmobwatdc_ ;
double *mc_ ;
double *dmcdc_ ;
double *porevol_ ;
double *deadporespace_ ;
double *porosity_ ;
double *dg_ ;
double *sw_ ;
double *c_ ;
double *ds_ ;
double *dsc_ ;
double *dcads_ ; // difference of cads to compute residual
double *dcadsdc_ ; // derivative of cads
double *pc_ ;
double *dpc_ ;
double *trans_ ;
@ -190,7 +212,11 @@ namespace Opm {
}
std::copy(porevol.begin(), porevol.end(), store_.porevol());
const std::vector<double> porosity = fluid.porosity();
std::copy(porosity.begin(), porosity.end(), store_.porosity());
store_.rockdensity() = fluid.rockdensity();
}
void makefhfQPeriodic( const std::vector<int>& p_faces,const std::vector<int>& hf_faces,
const std::vector<int>& nb_faces)
{
@ -347,15 +373,17 @@ namespace Opm {
) const {
(void) g;
const double pv = store_.porevol(cell);
const double dps = store_.deadporespace(cell);
const double pv = store_.porevol(cell);
const double dps = store_.deadporespace();
const double rhor = fluid_.rockdensity();
const double poro = store_.porosity(cell);
F[0] += pv * store_.ds(cell);
F[1] += pv * (1 - dps) * store_.dsc(cell);
F[0] += pv * store_.ds(cell);
F[1] += pv * (1 - dps) * store_.dsc(cell) + rhor*(1 - poro)/poro*pv*store_.dcads(cell);
dF[0] += pv;
dF[1] += 0.;
dF[2] += pv * (1 - dps) * store_.c(cell);
dF[3] += pv * (1 - dps) * store_.sw(cell);
dF[3] += pv * (1 - dps) * store_.sw(cell) + rhor*(1 - poro)/poro*pv*store_.dcadsdc(cell);
}
template <class Grid ,
@ -475,6 +503,7 @@ namespace Opm {
sys.vector().solution();
const ::std::vector<double>& sat = state.saturation();
const ::std::vector<double>& cpoly = state.concentration();
const ::std::vector<double>& cmax = state.maxconcentration();
bool in_range = true;
for (int cell = 0; cell < g.number_of_cells; ++cell) {
@ -485,8 +514,13 @@ namespace Opm {
store_.sw(cell) = s[0];
store_.c(cell) = c;
store_.dsc(cell) = s[0]*c - sat[cell*2 + 0]*cpoly[cell];
double dcadsdc;
double cads;
fluid_.adsorbtion(cpoly[cell], cmax[cell], cads, dcadsdc);
store_.dcads(cell) = -cads;
fluid_.adsorbtion(c, cmax[cell], cads, dcadsdc);
store_.dcads(cell) += cads;
store_.dcadsdc(cell) = dcadsdc;
double s_min = fluid_.s_min(cell);
double s_max = fluid_.s_max(cell);
@ -543,12 +577,14 @@ namespace Opm {
const SolutionVector& x ,
ReservoirState& state) {
double *s = &state.saturation()[0*2 + 0];
double *c = &state.concentration()[0*1 + 0];
double *s = &state.saturation()[0*2 + 0];
double *c = &state.concentration()[0*1 + 0];
double *cmax = &state.maxconcentration()[0*1 + 0];
for (int cell = 0; cell < g.number_of_cells; ++cell, s += 2, c += 1) {
for (int cell = 0; cell < g.number_of_cells; ++cell, s += 2, c += 1, cmax +=1) {
s[0] += x[2*cell + 0];
c[0] += x[2*cell + 1];
cmax[0] = std::max(c[0], cmax[0]);
double s_min = fluid_.s_min(cell);
double s_max = fluid_.s_max(cell);
assert(s[0] >= s_min - sat_tol_);