Switch to design based on deriving from 'Model'.
This commit is contained in:
Bård Skaflestad
@@ -45,8 +45,10 @@ namespace Opm {
|
||||
namespace spu_2p {
|
||||
class ModelParameterStorage {
|
||||
public:
|
||||
ModelParameterStorage(int nc, int totconn, double drho = 0.0)
|
||||
: drho_(dhro), data_()
|
||||
ModelParameterStorage(int nc, int totconn)
|
||||
: drho_(0.0), mob_(0), dmob_(0),
|
||||
porevol_(0), dg_(0), ds_(0),
|
||||
data_()
|
||||
{
|
||||
size_t alloc_sz;
|
||||
|
||||
@@ -93,218 +95,6 @@ namespace Opm {
|
||||
|
||||
std::vector<double> data_;
|
||||
};
|
||||
|
||||
class ModelParameters {
|
||||
public:
|
||||
template <class Grid>
|
||||
ModelParameters(const Grid* g,
|
||||
const ModelParameterStorage* s,
|
||||
const bool gravity)
|
||||
: gravity_(gravity),
|
||||
f2hf_ (2 * g->number_of_faces, -1),
|
||||
data_ (s)
|
||||
{
|
||||
for (int c = 0, i = 0; c < g->number_of_cells; ++c) {
|
||||
for (; i < g->cell_facepos[c + 1]; ++i) {
|
||||
const int f = g->cell_faces[i];
|
||||
const int p = 1 - (g->face_cells[2*f + 0] == c);
|
||||
|
||||
f2hf_[2*f + p] = i;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
computeStaticGravity(const Grid* g ,
|
||||
const double* grav ,
|
||||
const double* trans) {
|
||||
const int d = g->dimensions;
|
||||
|
||||
for (int c = 0, i = 0; c < g->number_of_cells; ++c) {
|
||||
const double* cc = g->cell_centroids + (c * d);
|
||||
|
||||
for (; i < g->cell_facepos[c + 1]; ++i) {
|
||||
const int f = g->cell_faces[i];
|
||||
const double* fc = g->face_centroids + (f * d);
|
||||
|
||||
double dg = 0.0;
|
||||
for (int j = 0; j < d; ++j) {
|
||||
dg += grav[j] * (fc[j] - cc[j]);
|
||||
}
|
||||
|
||||
data_->dg(i) = trans[f] * dg;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
initResidual(const int c, double* F) const {
|
||||
(void) c; // Suppress 'unused' warning
|
||||
*F = 0.0;
|
||||
}
|
||||
|
||||
template <class ReservoirState,
|
||||
class Grid >
|
||||
void
|
||||
fluxConnection(const ReservoirState* state,
|
||||
const Grid* g ,
|
||||
const double dt ,
|
||||
const int c ,
|
||||
const int f ,
|
||||
double* J1 ,
|
||||
double* J2 ,
|
||||
double* F ) const {
|
||||
|
||||
const int *n = g->face_cells + (2 * f);
|
||||
int pix[2];
|
||||
double sgn, mt, f1, dflux, gflux;
|
||||
double m[2], dm[2], *J[2];
|
||||
|
||||
dflux = state->faceflux[f];
|
||||
gflux = gravityFlux(f);
|
||||
upwindMobility(dflux, gflux, n, pix, m, dm);
|
||||
|
||||
assert (! ((m[0] < 0) || (m[1] < 0)));
|
||||
|
||||
mt = m[0] + m[1];
|
||||
assert (mt > 0);
|
||||
|
||||
f1 = m[0] / mt;
|
||||
|
||||
sgn = 2.0*(c1 == c) - 1.0;
|
||||
dflux *= sgn;
|
||||
gflux *= sgn;
|
||||
|
||||
const double v1 = dflux + m[1]*gflux;
|
||||
|
||||
// Assemble residual contributions
|
||||
*F += dt * f1 * v1;
|
||||
|
||||
// Assemble Jacobian (J1 <-> c, J2 <-> other)
|
||||
if (n[0] == c) { J[0] = J1; J[1] = J2; }
|
||||
else { J[0] = J2; J[1] = J1; }
|
||||
|
||||
// dF/dm_1 \cdot dm_1/ds
|
||||
*J[ pix[0] ] += dt * (1 - f1) / mt * v1 * dm[0];
|
||||
|
||||
/* dF/dm_2 \cdot dm_2/ds */
|
||||
*J[ pix[1] ] -= dt * f1 / mt * v1 * dm[1];
|
||||
*J[ pix[1] ] += dt * f1 * gflux * dm[1];
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
accumulation(const Grid* g,
|
||||
const int c,
|
||||
double* J,
|
||||
double* F) const {
|
||||
(void) g;
|
||||
|
||||
const double pv = data_->porevol(c);
|
||||
|
||||
*J += pv;
|
||||
*F += pv * data_->ds(c);
|
||||
}
|
||||
|
||||
template <class Grid ,
|
||||
class SourceTerms>
|
||||
void
|
||||
sourceTerms(const Grid* g ,
|
||||
const SourceTerms* src,
|
||||
const int i ,
|
||||
const double dt ,
|
||||
double* J ,
|
||||
double* F ) const {
|
||||
|
||||
double dflux = -src->flux[i]; // ->flux[] is rate of *inflow*
|
||||
|
||||
if (dflux < 0) {
|
||||
// src -> cell, affects residual only.
|
||||
*F += dt * dflux * src->saturation[2*i + 0];
|
||||
} else {
|
||||
// cell -> src
|
||||
const int c = src->cell[i];
|
||||
const double* m = data_->mob (c);
|
||||
const double* dm = data_->dmob(c);
|
||||
|
||||
const double mt = m[0] + m[1];
|
||||
|
||||
assert (! ((m[0] < 0) || (m[1] < 0)));
|
||||
assert (mt > 0);
|
||||
|
||||
const double f = m[0] / mt;
|
||||
const double df = ((1 - f)*dm[0] - f*dm[1]) / mt;
|
||||
|
||||
*F += dt * dflux * f;
|
||||
*J += dt * dflux * df;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
void
|
||||
upwindMobility(const double dflux,
|
||||
const double gflux,
|
||||
const int* n ,
|
||||
int* pix ,
|
||||
double* m ,
|
||||
double* dm ) const {
|
||||
bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) ||
|
||||
( (! (dflux > 0)) && (! (gflux > 0)) );
|
||||
|
||||
if (equal_sign) {
|
||||
|
||||
if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
|
||||
else { pix[0] = 1; }
|
||||
|
||||
m[0] = data_->mob(n[ pix[0] ]) [ 0 ];
|
||||
|
||||
if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
|
||||
else { pix[1] = 1; }
|
||||
|
||||
m[1] = data_->mob(n[ pix[1] ]) [ 1 ];
|
||||
|
||||
} else {
|
||||
|
||||
if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
|
||||
else { pix[1] = 1; }
|
||||
|
||||
m[1] = data_->mob(n[ pix[1] ]) [ 1 ];
|
||||
|
||||
if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
|
||||
else { pix[0] = 1; }
|
||||
|
||||
m[0] = data_->mob(n[ pix[0] ]) [ 0 ];
|
||||
}
|
||||
|
||||
dm[0] = data_->dmob(n[ pix[0] ]) [ 0 ];
|
||||
dm[1] = data_->dmob(n[ pix[1] ]) [ 1 ];
|
||||
}
|
||||
|
||||
double
|
||||
gravityFlux(const int f) const {
|
||||
double gflux;
|
||||
|
||||
if (gravity_) {
|
||||
int i1 = f2hf_[2*f + 0];
|
||||
int i2 = f2hf_[2*f + 1];
|
||||
|
||||
assert ((i1 >= 0) && (i2 >= 0));
|
||||
|
||||
gflux = data_->dg(i1) - data_->dg(i2);
|
||||
gflux *= data_->drho();
|
||||
} else {
|
||||
gflux = 0.0;
|
||||
}
|
||||
|
||||
return gflux;
|
||||
}
|
||||
|
||||
bool gravity_;
|
||||
std::vector<int> f2hf_ ;
|
||||
const ModelParameterStorage* data_ ;
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
@@ -312,30 +102,252 @@ namespace Opm {
|
||||
class SinglePointUpwindTwoPhase : private TwophaseFluid {
|
||||
public:
|
||||
template <class Grid>
|
||||
SinglePointUpwindTwoPhase(const Grid* g ,
|
||||
SinglePointUpwindTwoPhase(const Grid& g ,
|
||||
const std::vector<double>& porevol ,
|
||||
const double* grav = 0,
|
||||
const double* trans = 0)
|
||||
: TwophaseFluid(),
|
||||
store_(g->number_of_cells,
|
||||
g->cell_facepos[ g->number_of_cells ]),
|
||||
model_(g, &store_, (grav == 0) || (trans == 0))
|
||||
: TwophaseFluid() ,
|
||||
gravity_ ((grav != 0) && (trans != 0)) ,
|
||||
f2hf_ (2 * g.number_of_faces, -1) ,
|
||||
store_ (g.number_of_cells,
|
||||
g.cell_facepos[ g.number_of_cells ]),
|
||||
{
|
||||
if ((grav != 0) && (trans != 0)) {
|
||||
store_.drho() = density(0) - density(1);
|
||||
if (gravity_) {
|
||||
store_.drho() = this->density(0) - this->density(1);
|
||||
|
||||
model_.computeStaticGravity(g, grav, trans);
|
||||
this->computeStaticGravity(g, grav, trans);
|
||||
}
|
||||
|
||||
std::copy(porevol.begin(), porevol.end(), model_.porevol());
|
||||
for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
|
||||
for (; i < g.cell_facepos[c + 1]; ++i) {
|
||||
const int f = g.cell_faces[i];
|
||||
const int p = 1 - (g.face_cells[2*f + 0] == c);
|
||||
|
||||
f2hf_[2*f + p] = i;
|
||||
}
|
||||
}
|
||||
|
||||
std::copy(porevol.begin(), porevol.end(), data_.porevol());
|
||||
}
|
||||
|
||||
// -----------------------------------------------------------------
|
||||
// System assembly innards
|
||||
// -----------------------------------------------------------------
|
||||
void
|
||||
initResidual(const int c, double* F) const {
|
||||
(void) c; // Suppress 'unused' warning
|
||||
*F = 0.0;
|
||||
}
|
||||
|
||||
template <class ReservoirState,
|
||||
class Grid >
|
||||
void
|
||||
fluxConnection(const ReservoirState& state,
|
||||
const Grid& g ,
|
||||
const double dt ,
|
||||
const int c ,
|
||||
const int f ,
|
||||
double* J1 ,
|
||||
double* J2 ,
|
||||
double* F ) const {
|
||||
|
||||
const int *n = g.face_cells + (2 * f);
|
||||
double dflux = state.faceflux[f];
|
||||
double gflux = gravityFlux(f);
|
||||
|
||||
int pix[2];
|
||||
double m[2], dm[2];
|
||||
upwindMobility(dflux, gflux, n, pix, m, dm);
|
||||
|
||||
assert (! ((m[0] < 0) || (m[1] < 0)));
|
||||
|
||||
double mt = m[0] + m[1];
|
||||
assert (mt > 0);
|
||||
|
||||
double sgn = 2.0*(c1 == c) - 1.0;
|
||||
dflux *= sgn;
|
||||
gflux *= sgn;
|
||||
|
||||
double f1 = m[0] / mt;
|
||||
const double v1 = dflux + m[1]*gflux;
|
||||
|
||||
// Assemble residual contributions
|
||||
*F += dt * f1 * v1;
|
||||
|
||||
// Assemble Jacobian (J1 <-> c, J2 <-> other)
|
||||
double *J[2];
|
||||
if (n[0] == c) { J[0] = J1; J[1] = J2; }
|
||||
else { J[0] = J2; J[1] = J1; }
|
||||
|
||||
// dF/dm_1 \cdot dm_1/ds
|
||||
*J[ pix[0] ] += dt * (1 - f1) / mt * v1 * dm[0];
|
||||
|
||||
/* dF/dm_2 \cdot dm_2/ds */
|
||||
*J[ pix[1] ] -= dt * f1 / mt * v1 * dm[1];
|
||||
*J[ pix[1] ] += dt * f1 * gflux * dm[1];
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
accumulation(const Grid& g,
|
||||
const int c,
|
||||
double* J,
|
||||
double* F) const {
|
||||
(void) g;
|
||||
|
||||
const double pv = data_.porevol(c);
|
||||
|
||||
*J += pv;
|
||||
*F += pv * data_.ds(c);
|
||||
}
|
||||
|
||||
template <class Grid ,
|
||||
class SourceTerms>
|
||||
void
|
||||
sourceTerms(const Grid& g ,
|
||||
const SourceTerms& src,
|
||||
const int i ,
|
||||
const double dt ,
|
||||
double* J ,
|
||||
double* F ) const {
|
||||
|
||||
double dflux = -src.flux[i]; // .flux[] is rate of *inflow*
|
||||
|
||||
if (dflux < 0) {
|
||||
// src -> cell, affects residual only.
|
||||
*F += dt * dflux * src.saturation[2*i + 0];
|
||||
} else {
|
||||
// cell -> src
|
||||
const int c = src.cell[i];
|
||||
const double* m = data_.mob (c);
|
||||
const double* dm = data_.dmob(c);
|
||||
|
||||
const double mt = m[0] + m[1];
|
||||
|
||||
assert (! ((m[0] < 0) || (m[1] < 0)));
|
||||
assert (mt > 0);
|
||||
|
||||
const double f = m[0] / mt;
|
||||
const double df = ((1 - f)*dm[0] - f*dm[1]) / mt;
|
||||
|
||||
*F += dt * dflux * f;
|
||||
*J += dt * dflux * df;
|
||||
}
|
||||
}
|
||||
|
||||
// -----------------------------------------------------------------
|
||||
// Newton control
|
||||
// -----------------------------------------------------------------
|
||||
template <class ReservoirState,
|
||||
class Grid ,
|
||||
class JacobianSystem>
|
||||
void
|
||||
initStep(const ReservoirState& state,
|
||||
const Grid& g ,
|
||||
JacobianSystem& sys ) {
|
||||
(void) state; (void) g; // Suppress 'unused'
|
||||
|
||||
sys.vector().solution().fill(0.0);
|
||||
}
|
||||
|
||||
template <class ReservoirState,
|
||||
class Grid ,
|
||||
class JacobianSystem>
|
||||
void
|
||||
initIteration(const ReservoirState& state,
|
||||
const Grid& g ,
|
||||
JacobianSystem& sys ) {
|
||||
}
|
||||
|
||||
typedef spu_2p::ModelParameters model_type;
|
||||
const model_type* model() const { return &model_; }
|
||||
|
||||
private:
|
||||
spu_2p::ModelParameterStorage store_;
|
||||
spu_2p::ModelParameters model_;
|
||||
void
|
||||
upwindMobility(const double dflux,
|
||||
const double gflux,
|
||||
const int* n ,
|
||||
int* pix ,
|
||||
double* m ,
|
||||
double* dm ) const {
|
||||
bool equal_sign = ( (! (dflux < 0)) && (! (gflux < 0)) ) ||
|
||||
( (! (dflux > 0)) && (! (gflux > 0)) );
|
||||
|
||||
if (equal_sign) {
|
||||
|
||||
if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
|
||||
else { pix[0] = 1; }
|
||||
|
||||
m[0] = data_.mob(n[ pix[0] ]) [ 0 ];
|
||||
|
||||
if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
|
||||
else { pix[1] = 1; }
|
||||
|
||||
m[1] = data_.mob(n[ pix[1] ]) [ 1 ];
|
||||
|
||||
} else {
|
||||
|
||||
if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
|
||||
else { pix[1] = 1; }
|
||||
|
||||
m[1] = data_.mob(n[ pix[1] ]) [ 1 ];
|
||||
|
||||
if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
|
||||
else { pix[0] = 1; }
|
||||
|
||||
m[0] = data_.mob(n[ pix[0] ]) [ 0 ];
|
||||
}
|
||||
|
||||
dm[0] = data_.dmob(n[ pix[0] ]) [ 0 ];
|
||||
dm[1] = data_.dmob(n[ pix[1] ]) [ 1 ];
|
||||
}
|
||||
|
||||
template <class Grid>
|
||||
void
|
||||
computeStaticGravity(const Grid& g ,
|
||||
const double* grav ,
|
||||
const double* trans) {
|
||||
const int d = g.dimensions;
|
||||
|
||||
for (int c = 0, i = 0; c < g.number_of_cells; ++c) {
|
||||
const double* cc = g.cell_centroids + (c * d);
|
||||
|
||||
for (; i < g.cell_facepos[c + 1]; ++i) {
|
||||
const int f = g.cell_faces[i];
|
||||
const double* fc = g.face_centroids + (f * d);
|
||||
|
||||
double dg = 0.0;
|
||||
for (int j = 0; j < d; ++j) {
|
||||
dg += grav[j] * (fc[j] - cc[j]);
|
||||
}
|
||||
|
||||
data_.dg(i) = trans[f] * dg;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
double
|
||||
gravityFlux(const int f) const {
|
||||
double gflux;
|
||||
|
||||
if (gravity_) {
|
||||
int i1 = f2hf_[2*f + 0];
|
||||
int i2 = f2hf_[2*f + 1];
|
||||
|
||||
assert ((i1 >= 0) && (i2 >= 0));
|
||||
|
||||
gflux = data_.dg(i1) - data_.dg(i2);
|
||||
gflux *= data_.drho();
|
||||
} else {
|
||||
gflux = 0.0;
|
||||
}
|
||||
|
||||
return gflux;
|
||||
}
|
||||
|
||||
TwophaseFluid fluid_ ;
|
||||
|
||||
bool gravity_;
|
||||
std::vector<int> f2hf_ ;
|
||||
spu_2p::ModelParameterStorage store_ ;
|
||||
};
|
||||
}
|
||||
#endif /* OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER */
|
||||
|
||||
Reference in New Issue
Block a user