OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Move to reference-storing rather than private inheritance model.

Browse Source 
Fix various build problems while here.

Finally, add default Newton-control parameters to the NRControl structure.
This commit is contained in:
Bård Skaflestad 2011-10-05 20:09:15 +02:00
parent 72953ac126
commit 5202631f34
3 changed files with 109 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
src/ImplicitAssembly.hpp
View File

@ -36,24 +36,28 @@
#ifndef OPM_IMPLICITASSEMBLY_HPP_HEADER #ifndef OPM_IMPLICITASSEMBLY_HPP_HEADER
#define OPM_IMPLICITASSEMBLY_HPP_HEADER #define OPM_IMPLICITASSEMBLY_HPP_HEADER
#include <cstddef>
#include <algorithm> #include <algorithm>
#include <vector> #include <vector>
namespace Opm { namespace Opm {
template <class Model> template <class Model>
class ImplicitAssembly : private Model { class ImplicitAssembly {
enum { DofPerCell = Model::DofPerCell }; enum { DofPerCell = Model::DofPerCell };
public: public:
ImplicitAssembly(Model& model)
: model_(model)
{}
template <class Grid , template <class Grid ,
class JacobianSystem> class JacobianSystem>
void void
createSystem(const Grid& g , createSystem(const Grid& g ,
JacobianSystem& sys) const { JacobianSystem& sys) const {
typedef typename JacobianSystem::matrix_type::size_type sz_t; ::std::size_t m = g.number_of_cells;
sz_t m = g.number_of_cells; ::std::size_t nnz = g.number_of_cells + countConnections(g);
sz_t nnz = g.number_of_cells + countConnections(g);
sys.setSize(DofPerCell, m, nnz); sys.setSize(DofPerCell, m, nnz);
} }
@ -65,12 +69,12 @@ namespace Opm {
void void
assemble(const ReservoirState& state, assemble(const ReservoirState& state,
const Grid& g , const Grid& g ,
const SourceTerms& src , const SourceTerms* src ,
const double dt , const double dt ,
JacobianSystem& sys ) { JacobianSystem& sys ) {
for (int c = 0; c < g->number_of_cells; ++c) { for (int c = 0; c < g.number_of_cells; ++c) {
this->computeCellContrib(g, c, dt); this->computeCellContrib(state, g, c, dt);
this->assembleCellContrib(g, c, sys); this->assembleCellContrib(g, c, sys);
} }
@ -122,17 +126,17 @@ namespace Opm {
nconn_ = countConnections(g, c); nconn_ = countConnections(g, c);
connections_.resize (0); connections_.resize (0);
connections_.reserve (nconn + 1); connections_.reserve (nconn_ + 1);
connections_.push_back(c); connections_.push_back(c);
asm_buffer_.resize((2*nconn + 1)*ndof2 + (nconn + 2)*ndof); asm_buffer_.resize((2*nconn_ + 1)*ndof2 + (nconn_ + 2)*ndof);
std::fill(asm_buffer_.begin(), asm_buffer_.end(), 0.0); std::fill(asm_buffer_.begin(), asm_buffer_.end(), 0.0);
double* F = &asm_buffer_[(2*nconn + 1) * ndof2]; double* F = &asm_buffer_[(2*nconn_ + 1) * ndof2];
double* J1 = &asm_buffer_[(0*nconn + 1) * ndof2]; double* J1 = &asm_buffer_[(0*nconn_ + 1) * ndof2];
double* J2 = J1 + (1*nconn + 0) * ndof2 ; double* J2 = J1 + (1*nconn_ + 0) * ndof2 ;
this->initResidual(c, F); model_.initResidual(c, F);
F += ndof; F += ndof;
for (int i = g.cell_facepos[c + 0]; for (int i = g.cell_facepos[c + 0];
@ -144,15 +148,18 @@ namespace Opm {
if ((c1 >= 0) && (c2 >= 0)) { if ((c1 >= 0) && (c2 >= 0)) {
connections_.push_back((c1 == c) ? c2 : c1); connections_.push_back((c1 == c) ? c2 : c1);
this->fluxConnection(state, g, dt, c, f, J1, J2, F); model_.fluxConnection(state, g, dt, c, f, J1, J2, F);
J1 += ndof2; J2 += ndof2; F += ndof; J1 += ndof2; J2 += ndof2; F += ndof;
} }
} }
this->accumulation(g, c, &asm_buffer_[0], F); model_.accumulation(g, c, &asm_buffer_[0], F);
return 1;
} }
template <class Grid, class System> template <class Grid, class System>
void
assembleCellContrib(const Grid& g , assembleCellContrib(const Grid& g ,
const int c , const int c ,
System& sys) const { System& sys) const {
@ -164,7 +171,7 @@ namespace Opm {
typedef std::vector<int>::size_type sz_t; typedef std::vector<int>::size_type sz_t;
const double* J1 = &asm_buffer_[0]; const double* J1 = &asm_buffer_[0];
const double* J2 = J1 + ((1*nconn + 1) * ndof2); const double* J2 = J1 + ((1*nconn_ + 1) * ndof2);
// Assemble contributions from accumulation term // Assemble contributions from accumulation term
sys.matasm().assembleBlock(ndof, c, c, J1); J1 += ndof2; sys.matasm().assembleBlock(ndof, c, c, J1); J1 += ndof2;
@ -173,8 +180,8 @@ namespace Opm {
for (int i = g.cell_facepos[c + 0]; for (int i = g.cell_facepos[c + 0];
i < g.cell_facepos[c + 1]; ++i) { i < g.cell_facepos[c + 1]; ++i) {
int f = g.cell_faces[i]; int f = g.cell_faces[i];
int c1 = g.face_cell[2*f + 0]; int c1 = g.face_cells[2*f + 0];
int c2 = g.face_cell[2*f + 1]; int c2 = g.face_cells[2*f + 1];
c2 = (c1 == c) ? c2 : c1; c2 = (c1 == c) ? c2 : c1;
@ -188,8 +195,8 @@ namespace Opm {
} }
// Assemble residual // Assemble residual
const double* F = &asm_buffer_[(2*nconn + 1) * ndof2]; const double* F = &asm_buffer_[(2*nconn_ + 1) * ndof2];
for (int conn = 0; conn < nconn + 2; ++conn, F += ndof) { for (int conn = 0; conn < nconn_ + 2; ++conn, F += ndof) {
sys.vector().assembleBlock(ndof, c, F); sys.vector().assembleBlock(ndof, c, F);
} }
} }
@ -197,27 +204,28 @@ namespace Opm {
template <class Grid, class SourceTerms, class System> template <class Grid, class SourceTerms, class System>
void void
assembleSourceContrib(const Grid& g, assembleSourceContrib(const Grid& g,
const SourceTerms& src, const SourceTerms* src,
const double dt, const double dt,
System& sys) { System& sys) {
const int ndof = DofPerCell; const int ndof = DofPerCell;
const int ndof2 = ndof * ndof; const int ndof2 = ndof * ndof;
for (int i = 0; i < src.nsrc; ++i) { for (int i = 0; i < src->nsrc; ++i) {
std::fill_n(asm_buffer_.begin(), ndof2 + ndof, 0.0); std::fill_n(asm_buffer_.begin(), ndof2 + ndof, 0.0);
double *J = &asm_buffer_[0]; double *J = &asm_buffer_[0];
double *F = J + ndof2; double *F = J + ndof2;
this->sourceTerms(g, src, i, dt, J, F); model_.sourceTerms(g, src, i, dt, J, F);
const int c = src.cell[i]; const int c = src->cell[i];
sys.matasm().assembleBlock(ndof, c, c, J); sys.matasm().assembleBlock(ndof, c, c, J);
sys.vector().assembleBlock(ndof, c, F); sys.vector().assembleBlock(ndof, c, F);
} }
} }
Model& model_ ;
int nconn_ ; int nconn_ ;
std::vector<int> connections_; std::vector<int> connections_;
std::vector<double> asm_buffer_ ; std::vector<double> asm_buffer_ ;

48
src/ImplicitTransport.hpp
View File

@ -41,6 +41,13 @@
namespace Opm { namespace Opm {
namespace ImplicitTransportDetails { namespace ImplicitTransportDetails {
struct NRControl { struct NRControl {
NRControl()
: max_it(1),
atol(1.0e-6),
rtol(5.0e-7),
dxtol(1.0e-8)
{}
int max_it; int max_it;
double atol ; double atol ;
double rtol ; double rtol ;
@ -55,21 +62,24 @@ namespace Opm {
}; };
} }
template <class Model , template <class Model ,
class JacobianSystem , class JacobianSystem ,
template <class> VNorm , template <class> class VNorm ,
template <class> VNeg = VectorNegater, template <class> class VNeg ,
template <class> VZero = VectorZero> template <class> class VZero >
class ImplicitTransport : private Model { class ImplicitTransport {
public: public:
ImplicitTransport() {} ImplicitTransport(Model& model)
: model_(model),
asm_ (model)
{}
template <class Grid , template <class Grid ,
class SourceTerms , class SourceTerms ,
class ReservoirState, class ReservoirState,
class LinearSolver > class LinearSolver >
void solve(const Grid& g , void solve(const Grid& g ,
const SourceTerms& src , const SourceTerms* src ,
const double dt , const double dt ,
const ImplicitTransportDetails::NRControl& ctrl , const ImplicitTransportDetails::NRControl& ctrl ,
ReservoirState& state , ReservoirState& state ,
@ -80,10 +90,10 @@ namespace Opm {
asm_.createSystem(g, sys_); asm_.createSystem(g, sys_);
VZero<vector_type>::zero(sys_.writableResidual()); VZero<vector_type>::zero(sys_.vector().writableResidual());
this->initStep(state, g, sys_); model_.initStep(state, g, sys_);
this->initIteration(state, g, sys_); model_.initIteration(state, g, sys_);
asm_.assemble(state, g, src, dt, sys_); asm_.assemble(state, g, src, dt, sys_);
@ -103,12 +113,13 @@ namespace Opm {
VNeg<vector_type>::negate(sys_.vector().writableIncrement()); VNeg<vector_type>::negate(sys_.vector().writableIncrement());
this->finishIteration(state, g, sys_.vector()); model_.finishIteration(state, g, sys_.vector());
nrm_dx = VNorm<vector_type>::norm(sys_.vector().increment()); rpt.norm_dx =
VNorm<vector_type>::norm(sys_.vector().increment());
sys_.vector().addIncrement(); sys_.vector().addIncrement();
this->initIteration(state, g, sys_); model_.initIteration(state, g, sys_);
asm_.assemble(state, g, src, dt, sys_); asm_.assemble(state, g, src, dt, sys_);
rpt.norm_res = rpt.norm_res =
@ -118,13 +129,15 @@ namespace Opm {
done = (rpt.norm_res < ctrl.atol) || done = (rpt.norm_res < ctrl.atol) ||
(rpt.norm_res < ctrl.rtol * nrm_res0) || (rpt.norm_res < ctrl.rtol * nrm_res0) ||
(rpt.norm_dx < ctrl.dxtol) ||
(rpt.nit == ctrl.max_it); (rpt.nit == ctrl.max_it);
} }
this->finisStep(g, sys_.vector().solution(), state); model_.finishStep(g, sys_.vector().solution(), state);
if (rpt.norm_res < ctrl.atol) { rpt.flag = 1; } if (rpt.norm_res < ctrl.atol) { rpt.flag = 1; }
else if (rpt.norm_res < ctrl.rtol * nrm_res0) { rpt.flag = 2; } else if (rpt.norm_res < ctrl.rtol * nrm_res0) { rpt.flag = 2; }
else if (rpt.norm_dx < ctrl.dxtol) { rpt.flag = 3; }
else { rpt.flag = -1; } else { rpt.flag = -1; }
} }
@ -132,13 +145,16 @@ namespace Opm {
ImplicitTransport (const ImplicitTransport&); ImplicitTransport (const ImplicitTransport&);
ImplicitTransport& operator=(const ImplicitTransport&); ImplicitTransport& operator=(const ImplicitTransport&);
#if 0
using Model::initStep; using Model::initStep;
using Model::initIteration; using Model::initIteration;
using Model::finishIteration; using Model::finishIteration;
using Model::finishStep; using Model::finishStep;
#endif
Model& model_;
ImplicitAssembly<Model> asm_; ImplicitAssembly<Model> asm_;
JacobianSystem& sys_; JacobianSystem sys_;
}; };
} }
#endif /* OPM_IMPLICITTRANSPORT_HPP_HEADER */ #endif /* OPM_IMPLICITTRANSPORT_HPP_HEADER */

77
src/SinglePointUpwindTwoPhase.hpp
View File

@ -100,21 +100,22 @@ namespace Opm {
template <class TwophaseFluid> template <class TwophaseFluid>
class SinglePointUpwindTwoPhase : private TwophaseFluid { class SinglePointUpwindTwoPhase {
public: public:
template <class Grid> template <class Grid>
SinglePointUpwindTwoPhase(const Grid& g , SinglePointUpwindTwoPhase(const TwophaseFluid& fluid ,
const Grid& g ,
const std::vector<double>& porevol , const std::vector<double>& porevol ,
const double* grav = 0, const double* grav = 0,
const double* trans = 0) const double* trans = 0)
: TwophaseFluid() , : fluid_ (fluid) ,
gravity_ ((grav != 0) && (trans != 0)) , gravity_((grav != 0) && (trans != 0)) ,
f2hf_ (2 * g.number_of_faces, -1) , f2hf_ (2 * g.number_of_faces, -1) ,
store_ (g.number_of_cells, store_ (g.number_of_cells,
g.cell_facepos[ g.number_of_cells ]), g.cell_facepos[ g.number_of_cells ])
{ {
if (gravity_) { if (gravity_) {
store_.drho() = this->density(0) - this->density(1); store_.drho() = fluid_.density(0) - fluid_.density(1);
this->computeStaticGravity(g, grav, trans); this->computeStaticGravity(g, grav, trans);
} }
@ -128,12 +129,15 @@ namespace Opm {
} }
} }
std::copy(porevol.begin(), porevol.end(), data_.porevol()); std::copy(porevol.begin(), porevol.end(), store_.porevol());
} }
// ----------------------------------------------------------------- // -----------------------------------------------------------------
// System assembly innards // System assembly innards
// ----------------------------------------------------------------- // -----------------------------------------------------------------
enum { DofPerCell = 1 };
void void
initResidual(const int c, double* F) const { initResidual(const int c, double* F) const {
(void) c; // Suppress 'unused' warning (void) c; // Suppress 'unused' warning
@ -153,7 +157,7 @@ namespace Opm {
double* F ) const { double* F ) const {
const int *n = g.face_cells + (2 * f); const int *n = g.face_cells + (2 * f);
double dflux = state.faceflux[f]; double dflux = state.faceflux()[f];
double gflux = gravityFlux(f); double gflux = gravityFlux(f);
int pix[2]; int pix[2];
@ -165,7 +169,7 @@ namespace Opm {
double mt = m[0] + m[1]; double mt = m[0] + m[1];
assert (mt > 0); assert (mt > 0);
double sgn = 2.0*(c1 == c) - 1.0; double sgn = 2.0*(n[0] == c) - 1.0;
dflux *= sgn; dflux *= sgn;
gflux *= sgn; gflux *= sgn;
@ -196,32 +200,34 @@ namespace Opm {
double* F) const { double* F) const {
(void) g; (void) g;
const double pv = data_.porevol(c); const double pv = store_.porevol(c);
*J += pv; *J += pv;
*F += pv * data_.ds(c); *F += pv * store_.ds(c);
} }
template <class Grid , template <class Grid ,
class SourceTerms> class SourceTerms>
void void
sourceTerms(const Grid& g , sourceTerms(const Grid& g ,
const SourceTerms& src, const SourceTerms* src,
const int i , const int i ,
const double dt , const double dt ,
double* J , double* J ,
double* F ) const { double* F ) const {
double dflux = -src.flux[i]; // .flux[] is rate of *inflow* (void) g;
double dflux = -src->flux[i]; // ->flux[] is rate of *inflow*
if (dflux < 0) { if (dflux < 0) {
// src -> cell, affects residual only. // src -> cell, affects residual only.
*F += dt * dflux * src.saturation[2*i + 0]; *F += dt * dflux * src->saturation[2*i + 0];
} else { } else {
// cell -> src // cell -> src
const int c = src.cell[i]; const int c = src->cell[i];
const double* m = data_.mob (c); const double* m = store_.mob (c);
const double* dm = data_.dmob(c); const double* dm = store_.dmob(c);
const double mt = m[0] + m[1]; const double mt = m[0] + m[1];
@ -248,7 +254,8 @@ namespace Opm {
JacobianSystem& sys ) { JacobianSystem& sys ) {
(void) state; (void) g; // Suppress 'unused' (void) state; (void) g; // Suppress 'unused'
sys.vector().solution().fill(0.0); // sys.vector().solution().fill(0.0);
(void) sys;
} }
template <class ReservoirState, template <class ReservoirState,
@ -260,14 +267,19 @@ namespace Opm {
JacobianSystem& sys ) { JacobianSystem& sys ) {
std::array<double, 2 > s ; std::array<double, 2 > s ;
std::array<double, 2 > mob ;
std::array<double, 2*2> dmob; std::array<double, 2*2> dmob;
const ::std::vector<double>& sat = state.saturation();
for (int c = 0; c < g.number_of_cells; ++c) { for (int c = 0; c < g.number_of_cells; ++c) {
s[0] = state.saturation[c*2 + 0] + sys.vector().solution()[c]; s[0] = sat[c*2 + 0] + sys.vector().solution()[c];
s[1] = 1 - s[0]; s[1] = 1 - s[0];
this->mobility(c, s, store_.mob(c), dmob); fluid_.mobility(c, s, mob, dmob);
store_.mob (c)[0] = mob [0];
store_.mob (c)[1] = mob [1];
store_.dmob(c)[0] = dmob[0*2 + 0]; store_.dmob(c)[0] = dmob[0*2 + 0];
store_.dmob(c)[1] = -dmob[1*2 + 1]; store_.dmob(c)[1] = -dmob[1*2 + 1];
} }
@ -292,7 +304,7 @@ namespace Opm {
const SolutionVector& x , const SolutionVector& x ,
ReservoirState& state) { ReservoirState& state) {
double *s = &state.saturation[0*2 + 0]; double *s = &state.saturation()[0*2 + 0];
for (int c = 0; c < g.number_of_cells; ++c, s += 2) { for (int c = 0; c < g.number_of_cells; ++c, s += 2) {
s[0] += x[c] ; s[0] += x[c] ;
@ -316,28 +328,28 @@ namespace Opm {
if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; } if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; }
else { pix[0] = 1; } else { pix[0] = 1; }
m[0] = data_.mob(n[ pix[0] ]) [ 0 ]; m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; } if (! (dflux - m[0]*gflux < 0)) { pix[1] = 0; }
else { pix[1] = 1; } else { pix[1] = 1; }
m[1] = data_.mob(n[ pix[1] ]) [ 1 ]; m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
} else { } else {
if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; } if (! (dflux < 0) && ! (gflux > 0)) { pix[1] = 0; }
else { pix[1] = 1; } else { pix[1] = 1; }
m[1] = data_.mob(n[ pix[1] ]) [ 1 ]; m[1] = store_.mob(n[ pix[1] ]) [ 1 ];
if (dflux + m[1]*gflux > 0) { pix[0] = 0; } if (dflux + m[1]*gflux > 0) { pix[0] = 0; }
else { pix[0] = 1; } else { pix[0] = 1; }
m[0] = data_.mob(n[ pix[0] ]) [ 0 ]; m[0] = store_.mob(n[ pix[0] ]) [ 0 ];
} }
dm[0] = data_.dmob(n[ pix[0] ]) [ 0 ]; dm[0] = store_.dmob(n[ pix[0] ]) [ 0 ];
dm[1] = data_.dmob(n[ pix[1] ]) [ 1 ]; dm[1] = store_.dmob(n[ pix[1] ]) [ 1 ];
} }
template <class Grid> template <class Grid>
@ -359,7 +371,7 @@ namespace Opm {
dg += grav[j] * (fc[j] - cc[j]); dg += grav[j] * (fc[j] - cc[j]);
} }
data_.dg(i) = trans[f] * dg; store_.dg(i) = trans[f] * dg;
} }
} }
} }
@ -374,8 +386,8 @@ namespace Opm {
assert ((i1 >= 0) && (i2 >= 0)); assert ((i1 >= 0) && (i2 >= 0));
gflux = data_.dg(i1) - data_.dg(i2); gflux = store_.dg(i1) - store_.dg(i2);
gflux *= data_.drho(); gflux *= store_.drho();
} else { } else {
gflux = 0.0; gflux = 0.0;
} }
@ -383,6 +395,7 @@ namespace Opm {
return gflux; return gflux;
} }
TwophaseFluid fluid_ ;
bool gravity_; bool gravity_;
std::vector<int> f2hf_ ; std::vector<int> f2hf_ ;
spu_2p::ModelParameterStorage store_ ; spu_2p::ModelParameterStorage store_ ;