Add trivial, partial implementation of an implicit solver.
This commit is contained in:
Bård Skaflestad
210
src/ImplicitAssembly.hpp
Normal file
210
src/ImplicitAssembly.hpp
Normal file
@@ -0,0 +1,210 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: ImplicitAssembly.hpp
|
||||
//
|
||||
// Created: 2011-09-28 10:00:46+0200
|
||||
//
|
||||
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
|
||||
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
|
||||
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
|
||||
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
|
||||
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
|
||||
//
|
||||
//===========================================================================
|
||||
|
||||
|
||||
/*
|
||||
Copyright 2011 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2011 Statoil ASA.
|
||||
|
||||
This file is part of the Open Porous Media (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_IMPLICITASSEMBLY_HPP_HEADER
|
||||
#define OPM_IMPLICITASSEMBLY_HPP_HEADER
|
||||
|
||||
#include <algorithm>
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
template <int DofPerCell = 1>
|
||||
class ImplicitAssembly {
|
||||
public:
|
||||
template <class ReservoirState,
|
||||
class Grid ,
|
||||
class Model ,
|
||||
class SourceTerms ,
|
||||
class JacobianSystem>
|
||||
void
|
||||
assemble(const ReservoirState* state,
|
||||
const Grid* g ,
|
||||
const Model* model,
|
||||
const SourceTerms* src ,
|
||||
const double dt ,
|
||||
JacobianSystem* sys ) {
|
||||
|
||||
size_t m = g->number_of_cells * DofPerCell;
|
||||
|
||||
sys->matrix()->setSize(m, m);
|
||||
sys->vector()->setSize(m);
|
||||
|
||||
for (int c = 0; c < g->number_of_cells; ++c) {
|
||||
const int nconn = computeCellContrib(g, model, c, dt);
|
||||
assembleCellContrib(g, c, nconn, sys);
|
||||
}
|
||||
|
||||
sys->matrix()->finalizeStructure();
|
||||
|
||||
if (src != 0) {
|
||||
assembleSourceContrib(g, model, src, dt, sys);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
template <class Grid>
|
||||
int
|
||||
countConnections(const Grid* g, int c) {
|
||||
int i, f, c1, c2, n;
|
||||
|
||||
for (i = g->cell_facepos[c + 0], n = 0;
|
||||
i < g->cell_facepos[c + 1]; ++i) {
|
||||
f = g->cell_faces[i];
|
||||
c1 = g->face_cells[2*f + 0];
|
||||
c2 = g->face_cells[2*f + 1];
|
||||
|
||||
n += (c1 >= 0) && (c2 >= 0);
|
||||
}
|
||||
|
||||
return n;
|
||||
}
|
||||
|
||||
template <class ReservoirState, class Grid, class Model>
|
||||
int
|
||||
computeCellContrib(const ReservoirState* state,
|
||||
const Grid* g ,
|
||||
const Model* model,
|
||||
const double dt ,
|
||||
const int c ) {
|
||||
const int ndof = DofPerCell;
|
||||
const int ndof2 = ndof * ndof;
|
||||
const int nconn = countConnections(g, c);
|
||||
|
||||
row_structure_.resize (0);
|
||||
row_structure_.reserve (nconn + 1);
|
||||
row_structure_.push_back(c);
|
||||
|
||||
asm_buffer_.resize((2*nconn + 1)*ndof2 + (nconn + 2)*ndof);
|
||||
std::fill(asm_buffer_.begin(), asm_buffer_.end(), 0.0);
|
||||
|
||||
double* F = &asm_buffer_[(2*nconn + 1) * ndof2];
|
||||
double* J1 = &asm_buffer_[(0*nconn + 1) * ndof2];
|
||||
double* J2 = J1 + (1*nconn + 0) * ndof2 ;
|
||||
|
||||
model->initResidual(c, F);
|
||||
F += ndof;
|
||||
|
||||
for (int i = g->cell_facepos[c + 0];
|
||||
i < g->cell_facepos[c + 1]; ++i) {
|
||||
int f = g->cell_faces[i];
|
||||
int c1 = g->face_cells[2*f + 0];
|
||||
int c2 = g->face_cells[2*f + 1];
|
||||
|
||||
if ((c1 >= 0) && (c2 >= 0)) {
|
||||
row_structure_.push_back((c1 == c) ? c2 : c1);
|
||||
|
||||
model->fluxConnection(state, g, dt, c, f, J1, J2, F);
|
||||
J1 += ndof2; J2 += ndof2; F += ndof;
|
||||
}
|
||||
}
|
||||
|
||||
model->accumulation(g, c, &asm_buffer_[0], F);
|
||||
|
||||
return nconn;
|
||||
}
|
||||
|
||||
template <class Grid, class System>
|
||||
assembleCellContrib(const Grid* g,
|
||||
const int c,
|
||||
const int nconn,
|
||||
System* sys) const {
|
||||
const int ndof = DofPerCell;
|
||||
const int ndof2 = ndof * ndof;
|
||||
|
||||
sys->matrix()->createBlockRow(c, connections_, ndof);
|
||||
|
||||
typedef std::vector<int>::size_type sz_t;
|
||||
|
||||
const double* J1 = &asm_buffer_[0];
|
||||
const double* J2 = J1 + ((1*nconn + 1) * ndof2);
|
||||
|
||||
// Assemble contributions from accumulation term
|
||||
sys->matrix()->assembleBlock(ndof, c, c, J1); J1 += ndof2;
|
||||
|
||||
// Assemble connection contributions.
|
||||
for (int i = g->cell_facepos[c + 0];
|
||||
i < g->cell_facepos[c + 1]; ++i) {
|
||||
int f = g->cell_faces[i];
|
||||
int c1 = g->face_cell[2*f + 0];
|
||||
int c2 = g->face_cell[2*f + 1];
|
||||
|
||||
c2 = (c1 == c) ? c2 : c1;
|
||||
|
||||
if (c2 >= 0) {
|
||||
sys->matrix()->assembleBlock(ndof, c, c , J1);
|
||||
sys->matrix()->assembleBlock(ndof, c, c2, J2);
|
||||
|
||||
J1 += ndof2;
|
||||
J2 += ndof2;
|
||||
}
|
||||
}
|
||||
|
||||
// Assemble residual
|
||||
const double* F = &asm_buffer_[(2*nconn + 1) * ndof2];
|
||||
for (int conn = 0; conn < nconn + 2; ++conn, F += ndof) {
|
||||
sys->vector.assembleBlock(ndof, c, F);
|
||||
}
|
||||
}
|
||||
|
||||
template <class Grid, class SourceTerms, class System>
|
||||
void
|
||||
assembleSourceContrib(const Grid* g,
|
||||
const Model* model,
|
||||
const SourceTerms* src,
|
||||
const double dt,
|
||||
System* sys) {
|
||||
const int ndof = DofPerCell;
|
||||
const int ndof2 = ndof * ndof;
|
||||
|
||||
for (int i = 0; i < src->nsrc; ++i) {
|
||||
std::fill_n(asm_buffer_.begin(), ndof2 + ndof, 0.0);
|
||||
|
||||
double *J = &asm_buffer_[0];
|
||||
double *F = J + ndof2;
|
||||
|
||||
model->sourceTerms(g, src, i, dt, J, F);
|
||||
|
||||
const int c = src->cell[i];
|
||||
|
||||
sys->matrix()->assembleBlock(ndof, c, c, J);
|
||||
sys->vector()->assembleBlock(ndof, c, F);
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<int> connections_;
|
||||
std::vector<double> asm_buffer_;
|
||||
};
|
||||
}
|
||||
#endif // OPM_IMPLICITASSEMBLY_HPP_HEADER
|
||||
341
src/SinglePointUpwindTwoPhase.hpp
Normal file
341
src/SinglePointUpwindTwoPhase.hpp
Normal file
@@ -0,0 +1,341 @@
|
||||
/*===========================================================================
|
||||
//
|
||||
// File: SinglePointUpwindTwoPhase.hpp
|
||||
//
|
||||
// Created: 2011-09-28 14:21:34+0200
|
||||
//
|
||||
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
|
||||
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
|
||||
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
|
||||
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
|
||||
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
|
||||
//
|
||||
//==========================================================================*/
|
||||
|
||||
|
||||
/*
|
||||
Copyright 2011 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2011 Statoil ASA.
|
||||
|
||||
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_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER
|
||||
#define OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER
|
||||
|
||||
#include <cassert>
|
||||
|
||||
#include <algorithm>
|
||||
#include <vector>
|
||||
|
||||
namespace Opm {
|
||||
namespace spu_2p {
|
||||
class ModelParameterStorage {
|
||||
public:
|
||||
ModelParameterStorage(int nc, int totconn, double drho = 0.0)
|
||||
: drho_(dhro), data_()
|
||||
{
|
||||
size_t alloc_sz;
|
||||
|
||||
alloc_sz = 2 * nc; // mob_
|
||||
alloc_sz += 2 * nc; // dmob_
|
||||
alloc_sz += 1 * nc; // porevol_
|
||||
alloc_sz += 1 * totconn; // dg_
|
||||
alloc_sz += 1 * nc; // ds_
|
||||
|
||||
data_.resize(alloc_sz);
|
||||
|
||||
mob_ = &data_[0];
|
||||
dmob_ = mob_ + (2 * nc );
|
||||
porevol_ = dmob_ + (2 * nc );
|
||||
dg_ = porevol_ + (1 * nc );
|
||||
ds_ = dg_ + (1 * totconn);
|
||||
}
|
||||
|
||||
double& drho () { return drho_ ; }
|
||||
double drho () const { return drho_ ; }
|
||||
|
||||
double* mob (int c) { return mob_ + (2*c + 0); }
|
||||
const double* mob (int c) const { return mob_ + (2*c + 0); }
|
||||
|
||||
double* dmob (int c) { return dmob_ + (2*c + 0); }
|
||||
const double* dmob (int c) const { return dmob_ + (2*c + 0); }
|
||||
|
||||
double* porevol() { return porevol_ ; }
|
||||
double porevol(int c) const { return porevol_[c] ; }
|
||||
|
||||
double& dg(int i) { return dg_[i] ; }
|
||||
double dg(int i) const { return dg_[i] ; }
|
||||
|
||||
double& ds(int c) { return ds_[c] ; }
|
||||
double ds(int c) const { return ds_[c] ; }
|
||||
|
||||
private:
|
||||
double drho_ ;
|
||||
double *mob_ ;
|
||||
double *dmob_ ;
|
||||
double *porevol_;
|
||||
double *dg_ ;
|
||||
double *ds_ ;
|
||||
|
||||
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_ ;
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
template <class TwophaseFluid>
|
||||
class SinglePointUpwindTwoPhase : private TwophaseFluid {
|
||||
public:
|
||||
template <class Grid>
|
||||
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))
|
||||
{
|
||||
if ((grav != 0) && (trans != 0)) {
|
||||
store_.drho() = density(0) - density(1);
|
||||
|
||||
model_.computeStaticGravity(g, grav, trans);
|
||||
}
|
||||
|
||||
std::copy(porevol.begin(), porevol.end(), model_.porevol());
|
||||
}
|
||||
|
||||
typedef spu_2p::ModelParameters model_type;
|
||||
const model_type* model() const { return &model_; }
|
||||
|
||||
private:
|
||||
spu_2p::ModelParameterStorage store_;
|
||||
spu_2p::ModelParameters model_;
|
||||
};
|
||||
}
|
||||
#endif /* OPM_SINGLEPOINTUPWINDTWOPHASE_HPP_HEADER */
|
||||
Reference in New Issue
Block a user