From 57fe2cf237e9d2009d34fd1a070b2bf284d811f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 May 2012 14:48:09 +0200 Subject: [PATCH 1/4] Update comments in implementation file. --- opm/core/pressure/CompressibleTpfa.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 52acc35f..f2c910af 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -44,7 +44,10 @@ namespace Opm /// \param[in] grid A 2d or 3d grid. /// \param[in] props Rock and fluid properties. /// \param[in] linsolver Linear solver to use. - /// \param[in] gravity Gravity vector. If nonzero, the array should + /// \param[in] residual_tol Solution accepted if inf-norm of residual is smaller. + /// \param[in] change_tol Solution accepted if inf-norm of change is smaller. + /// \param[in] maxiter Maximum acceptable + /// \param[in] gravity Gravity vector. If non-null, the array should /// have D elements. /// \param[in] wells The wells argument. Will be used in solution, /// is ignored if NULL. From 2bce30d9f7fed5d097deb7d83a7460949da98a40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 May 2012 16:39:05 +0200 Subject: [PATCH 2/4] Initialize WellState::bhp() to cell pressures of first perforation. --- opm/core/WellState.hpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/opm/core/WellState.hpp b/opm/core/WellState.hpp index 442b7533..56ede664 100644 --- a/opm/core/WellState.hpp +++ b/opm/core/WellState.hpp @@ -29,11 +29,18 @@ namespace Opm class WellState { public: - void init(const Wells* wells) + template + void init(const Wells* wells, const State& state) { if (wells) { - bhp_.resize(wells->number_of_wells); - perfrates_.resize(wells->well_connpos[wells->number_of_wells]); + const int nw = wells->number_of_wells; + bhp_.resize(nw); + // Initialize bhp to be pressure in first perforation cell. + for (int w = 0; w < nw; ++w) { + const int cell = wells->well_cells[wells->well_connpos[w]]; + bhp_[w] = state.pressure()[cell]; + } + perfrates_.resize(wells->well_connpos[nw]); } } From 4245c975366dfbdf6d8bae0556af741bf17825b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 May 2012 16:39:39 +0200 Subject: [PATCH 3/4] Update to match WellState::init() interface change. --- examples/sim_wateroil.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 1e8df834..cd6b996b 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -363,7 +363,7 @@ main(int argc, char** argv) int num_wells = 0; if (wells->c_wells()) { num_wells = wells->c_wells()->number_of_wells; - well_state.init(wells->c_wells()); + well_state.init(wells->c_wells(), state); well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(num_wells), 0.0); wellreport.push(*props, *wells->c_wells(), state.pressure(), state.surfacevol(), state.saturation(), From b5d3615fa555914548b380080e63cf9b4b749494 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 May 2012 16:42:03 +0200 Subject: [PATCH 4/4] Some initial testing of CompressibleTpfa done. Fixed bug in pressure update from increments. Change injection wellperf_phasemob_ to be same as cell mobilities. Improve iteration reporting. --- opm/core/pressure/CompressibleTpfa.cpp | 36 +++++++++++++------------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index f2c910af..0ee8bf8b 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -114,6 +114,7 @@ namespace Opm WellState& well_state) { const int nc = grid_.number_of_cells; + const int nw = wells_->number_of_wells; // Set up dynamic data. computePerSolveDynamicData(dt, state, well_state); @@ -127,20 +128,29 @@ namespace Opm double res_norm = residualNorm(); std::cout << "\nIteration Residual Change in p\n" << std::setw(9) << iter - << std::setw(18) << res_norm << std::endl; - for (; (iter < maxiter_) && (res_norm > residual_tol_); ++iter) { + << std::setw(18) << res_norm + << std::setw(18) << '*' << std::endl; + for (; (iter < maxiter_) && (res_norm > residual_tol_); ) { // Solve for increment in Newton method: // incr = x_{n+1} - x_{n} = -J^{-1}F // (J is Jacobian matrix, F is residual) solveIncrement(); + ++iter; // Update pressure vars with increment. - std::copy(pressure_increment_.begin(), pressure_increment_.begin() + nc, state.pressure().begin()); - std::copy(pressure_increment_.begin() + nc, pressure_increment_.end(), well_state.bhp().begin()); + for (int c = 0; c < nc; ++c) { + state.pressure()[c] += pressure_increment_[c]; + } + for (int w = 0; w < nw; ++w) { + well_state.bhp()[w] += pressure_increment_[nc + w]; + } // Stop iterating if increment is small. inc_norm = incrementNorm(); if (inc_norm <= change_tol_) { + std::cout << std::setw(9) << iter + << std::setw(18) << '*' + << std::setw(18) << inc_norm << std::endl; break; } @@ -153,7 +163,7 @@ namespace Opm // Update residual norm. res_norm = residualNorm(); - std::cout << std::setw(9) << iter + 1 + std::cout << std::setw(9) << iter << std::setw(18) << res_norm << std::setw(18) << inc_norm << std::endl; } @@ -411,11 +421,9 @@ namespace Opm wellperf_phasemob_.resize(nperf*np); // The A matrix is set equal to the perforation grid cells' // matrix, for both injectors and producers. - // The mobilities are all set equal to the total mobility for the - // cell for injectors, and equal to individual phase mobilities - // for producers. + // The mobilities are set equal to the perforation grid cells' + // mobilities, for both injectors and producers. for (int w = 0; w < nw; ++w) { - bool is_injector = wells_->type[w] == INJECTOR; for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) { const int c = wells_->well_cells[j]; const double* cA = &cell_A_[np*np*c]; @@ -423,15 +431,7 @@ namespace Opm std::copy(cA, cA + np*np, wpA); const double* cM = &cell_phasemob_[np*c]; double* wpM = &wellperf_phasemob_[np*j]; - if (is_injector) { - double totmob = 0.0; - for (int phase = 0; phase < np; ++phase) { - totmob += cM[phase]; - } - std::fill(wpM, wpM + np, totmob); - } else { - std::copy(cM, cM + np, wpM); - } + std::copy(cM, cM + np, wpM); } } }