diff --git a/src/ImplicitAssembly.hpp b/src/ImplicitAssembly.hpp index 0d79c7d9..d3033976 100644 --- a/src/ImplicitAssembly.hpp +++ b/src/ImplicitAssembly.hpp @@ -36,24 +36,28 @@ #ifndef OPM_IMPLICITASSEMBLY_HPP_HEADER #define OPM_IMPLICITASSEMBLY_HPP_HEADER +#include #include #include namespace Opm { template - class ImplicitAssembly : private Model { + class ImplicitAssembly { enum { DofPerCell = Model::DofPerCell }; public: + ImplicitAssembly(Model& model) + : model_(model) + {} + template void createSystem(const Grid& g , JacobianSystem& sys) const { - typedef typename JacobianSystem::matrix_type::size_type sz_t; - sz_t m = g.number_of_cells; - sz_t nnz = g.number_of_cells + countConnections(g); + ::std::size_t m = g.number_of_cells; + ::std::size_t nnz = g.number_of_cells + countConnections(g); sys.setSize(DofPerCell, m, nnz); } @@ -65,12 +69,12 @@ namespace Opm { void assemble(const ReservoirState& state, const Grid& g , - const SourceTerms& src , + const SourceTerms* src , const double dt , JacobianSystem& sys ) { - for (int c = 0; c < g->number_of_cells; ++c) { - this->computeCellContrib(g, c, dt); + for (int c = 0; c < g.number_of_cells; ++c) { + this->computeCellContrib(state, g, c, dt); this->assembleCellContrib(g, c, sys); } @@ -122,17 +126,17 @@ namespace Opm { nconn_ = countConnections(g, c); connections_.resize (0); - connections_.reserve (nconn + 1); + connections_.reserve (nconn_ + 1); 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); - double* F = &asm_buffer_[(2*nconn + 1) * ndof2]; - double* J1 = &asm_buffer_[(0*nconn + 1) * ndof2]; - double* J2 = J1 + (1*nconn + 0) * ndof2 ; + double* F = &asm_buffer_[(2*nconn_ + 1) * ndof2]; + double* J1 = &asm_buffer_[(0*nconn_ + 1) * ndof2]; + double* J2 = J1 + (1*nconn_ + 0) * ndof2 ; - this->initResidual(c, F); + model_.initResidual(c, F); F += ndof; for (int i = g.cell_facepos[c + 0]; @@ -144,15 +148,18 @@ namespace Opm { if ((c1 >= 0) && (c2 >= 0)) { 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; } } - this->accumulation(g, c, &asm_buffer_[0], F); + model_.accumulation(g, c, &asm_buffer_[0], F); + + return 1; } template + void assembleCellContrib(const Grid& g , const int c , System& sys) const { @@ -164,7 +171,7 @@ namespace Opm { typedef std::vector::size_type sz_t; 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 sys.matasm().assembleBlock(ndof, c, c, J1); J1 += ndof2; @@ -173,8 +180,8 @@ namespace Opm { 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]; + int c1 = g.face_cells[2*f + 0]; + int c2 = g.face_cells[2*f + 1]; c2 = (c1 == c) ? c2 : c1; @@ -188,8 +195,8 @@ namespace Opm { } // Assemble residual - const double* F = &asm_buffer_[(2*nconn + 1) * ndof2]; - for (int conn = 0; conn < nconn + 2; ++conn, F += ndof) { + 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); } } @@ -197,27 +204,28 @@ namespace Opm { template void assembleSourceContrib(const Grid& g, - const SourceTerms& src, + 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) { + 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; - 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.vector().assembleBlock(ndof, c, F); } } + Model& model_ ; int nconn_ ; std::vector connections_; std::vector asm_buffer_ ; diff --git a/src/ImplicitTransport.hpp b/src/ImplicitTransport.hpp index 41c5a24e..e5ea065f 100644 --- a/src/ImplicitTransport.hpp +++ b/src/ImplicitTransport.hpp @@ -41,6 +41,13 @@ namespace Opm { namespace ImplicitTransportDetails { struct NRControl { + NRControl() + : max_it(1), + atol(1.0e-6), + rtol(5.0e-7), + dxtol(1.0e-8) + {} + int max_it; double atol ; double rtol ; @@ -55,21 +62,24 @@ namespace Opm { }; } - template VNorm , - template VNeg = VectorNegater, - template VZero = VectorZero> - class ImplicitTransport : private Model { + template class VNorm , + template class VNeg , + template class VZero > + class ImplicitTransport { public: - ImplicitTransport() {} + ImplicitTransport(Model& model) + : model_(model), + asm_ (model) + {} template void solve(const Grid& g , - const SourceTerms& src , + const SourceTerms* src , const double dt , const ImplicitTransportDetails::NRControl& ctrl , ReservoirState& state , @@ -80,10 +90,10 @@ namespace Opm { asm_.createSystem(g, sys_); - VZero::zero(sys_.writableResidual()); + VZero::zero(sys_.vector().writableResidual()); - this->initStep(state, g, sys_); - this->initIteration(state, g, sys_); + model_.initStep(state, g, sys_); + model_.initIteration(state, g, sys_); asm_.assemble(state, g, src, dt, sys_); @@ -103,12 +113,13 @@ namespace Opm { VNeg::negate(sys_.vector().writableIncrement()); - this->finishIteration(state, g, sys_.vector()); + model_.finishIteration(state, g, sys_.vector()); - nrm_dx = VNorm::norm(sys_.vector().increment()); + rpt.norm_dx = + VNorm::norm(sys_.vector().increment()); sys_.vector().addIncrement(); - this->initIteration(state, g, sys_); + model_.initIteration(state, g, sys_); asm_.assemble(state, g, src, dt, sys_); rpt.norm_res = @@ -118,13 +129,15 @@ namespace Opm { done = (rpt.norm_res < ctrl.atol) || (rpt.norm_res < ctrl.rtol * nrm_res0) || + (rpt.norm_dx < ctrl.dxtol) || (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; } 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; } } @@ -132,13 +145,16 @@ namespace Opm { ImplicitTransport (const ImplicitTransport&); ImplicitTransport& operator=(const ImplicitTransport&); +#if 0 using Model::initStep; using Model::initIteration; using Model::finishIteration; using Model::finishStep; +#endif + Model& model_; ImplicitAssembly asm_; - JacobianSystem& sys_; + JacobianSystem sys_; }; } #endif /* OPM_IMPLICITTRANSPORT_HPP_HEADER */ diff --git a/src/SinglePointUpwindTwoPhase.hpp b/src/SinglePointUpwindTwoPhase.hpp index b444a275..00453eb8 100644 --- a/src/SinglePointUpwindTwoPhase.hpp +++ b/src/SinglePointUpwindTwoPhase.hpp @@ -100,21 +100,22 @@ namespace Opm { template - class SinglePointUpwindTwoPhase : private TwophaseFluid { + class SinglePointUpwindTwoPhase { public: template - SinglePointUpwindTwoPhase(const Grid& g , + SinglePointUpwindTwoPhase(const TwophaseFluid& fluid , + const Grid& g , const std::vector& porevol , const double* grav = 0, const double* 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 ]), + : fluid_ (fluid) , + 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 (gravity_) { - store_.drho() = this->density(0) - this->density(1); + store_.drho() = fluid_.density(0) - fluid_.density(1); 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 // ----------------------------------------------------------------- + + enum { DofPerCell = 1 }; + void initResidual(const int c, double* F) const { (void) c; // Suppress 'unused' warning @@ -153,7 +157,7 @@ namespace Opm { double* F ) const { const int *n = g.face_cells + (2 * f); - double dflux = state.faceflux[f]; + double dflux = state.faceflux()[f]; double gflux = gravityFlux(f); int pix[2]; @@ -165,7 +169,7 @@ namespace Opm { double mt = m[0] + m[1]; assert (mt > 0); - double sgn = 2.0*(c1 == c) - 1.0; + double sgn = 2.0*(n[0] == c) - 1.0; dflux *= sgn; gflux *= sgn; @@ -196,32 +200,34 @@ namespace Opm { double* F) const { (void) g; - const double pv = data_.porevol(c); + const double pv = store_.porevol(c); *J += pv; - *F += pv * data_.ds(c); + *F += pv * store_.ds(c); } template void sourceTerms(const Grid& g , - const SourceTerms& src, + const SourceTerms* src, const int i , const double dt , double* J , 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) { // src -> cell, affects residual only. - *F += dt * dflux * src.saturation[2*i + 0]; + *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 int c = src->cell[i]; + const double* m = store_.mob (c); + const double* dm = store_.dmob(c); const double mt = m[0] + m[1]; @@ -248,7 +254,8 @@ namespace Opm { JacobianSystem& sys ) { (void) state; (void) g; // Suppress 'unused' - sys.vector().solution().fill(0.0); + // sys.vector().solution().fill(0.0); + (void) sys; } template s ; + std::array mob ; std::array dmob; + const ::std::vector& sat = state.saturation(); + 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]; - 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)[1] = -dmob[1*2 + 1]; } @@ -292,7 +304,7 @@ namespace Opm { const SolutionVector& x , 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) { s[0] += x[c] ; @@ -316,28 +328,28 @@ namespace Opm { if (! (dflux < 0) && ! (gflux < 0)) { pix[0] = 0; } 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; } else { pix[1] = 1; } - m[1] = data_.mob(n[ pix[1] ]) [ 1 ]; + m[1] = store_.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 ]; + m[1] = store_.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 ]; + m[0] = store_.mob(n[ pix[0] ]) [ 0 ]; } - dm[0] = data_.dmob(n[ pix[0] ]) [ 0 ]; - dm[1] = data_.dmob(n[ pix[1] ]) [ 1 ]; + dm[0] = store_.dmob(n[ pix[0] ]) [ 0 ]; + dm[1] = store_.dmob(n[ pix[1] ]) [ 1 ]; } template @@ -359,7 +371,7 @@ namespace Opm { 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)); - gflux = data_.dg(i1) - data_.dg(i2); - gflux *= data_.drho(); + gflux = store_.dg(i1) - store_.dg(i2); + gflux *= store_.drho(); } else { gflux = 0.0; } @@ -383,6 +395,7 @@ namespace Opm { return gflux; } + TwophaseFluid fluid_ ; bool gravity_; std::vector f2hf_ ; spu_2p::ModelParameterStorage store_ ;