From 9cb7e8635e7f275a5fc95f482352f0ff39ee68a2 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 11 Mar 2014 13:21:51 +0100 Subject: [PATCH 01/18] Adds -= operator An elementwise -= operator is added to the autodiff class. --- opm/autodiff/AutoDiffBlock.hpp | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 927927826..a98c014be 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -194,6 +194,28 @@ namespace Opm return *this; } + /// Elementwise operator -= + AutoDiffBlock& operator-=(const AutoDiffBlock& rhs) + { + if (jac_.empty()) { + jac_ = rhs.jac_; + } else if (!rhs.jac_.empty()) { + assert (numBlocks() == rhs.numBlocks()); + assert (value().size() == rhs.value().size()); + + const int num_blocks = numBlocks(); + for (int block = 0; block < num_blocks; ++block) { + assert(jac_[block].rows() == rhs.jac_[block].rows()); + assert(jac_[block].cols() == rhs.jac_[block].cols()); + jac_[block] -= rhs.jac_[block]; + } + } + + val_ -= rhs.val_; + + return *this; + } + /// Elementwise operator + AutoDiffBlock operator+(const AutoDiffBlock& rhs) const { From a19aff63e7a9cc4947f4666621fb8cece96e406c Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 11 Mar 2014 14:28:16 +0100 Subject: [PATCH 02/18] Adds new well formulation Todo: incorporate WellDensitySegment. Currently values of the pressure drop is hardcoded to make the rest of the code work Todo: make it possible to shut perforation with crossflow. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 254 ++++++++++++++++--- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 9 + 2 files changed, 233 insertions(+), 30 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index f3f686712..c055377ad 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -410,6 +410,7 @@ namespace { const int nw = wells_.number_of_wells; // The transpose() below switches the ordering. const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); + const V qs = Eigen::Map(wrates.data(), nw*np); state.qs = ADB::constant(qs, bpat); @@ -676,6 +677,228 @@ namespace { } + addWellEq(state); + addWellControlEq(state); + + + } + + + void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state){ + + const int nc = grid_.number_of_cells; + const int np = wells_.number_of_phases; + const int nw = wells_.number_of_wells; + const int nperf = wells_.well_connpos[nw]; + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + V Tw = Eigen::Map(wells_.WI, nperf); + const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + + // pressure drop from solution + double tmp[2] = {-1.1319e-09, -2.0926e-09 }; + V cdp = Eigen::Map(tmp, 2); + + // Extract variables for perforation cell pressures + // and corresponding perforation well pressures. + const ADB p_perfcell = subset(state.pressure, well_cells); + + // Pressure drawdown (also used to determine direction of flow) + const ADB drawdown = p_perfcell - (wops_.w2p * state.bhp + cdp); + + // current injecting connections + auto connInjInx = drawdown.value() < 0; + + // injector == 1, producer == 0 + V isInj = V::Zero(nw); + for (int w = 0; w < nw; ++w) { + if (wells_.type[w] == INJECTOR) { + isInj[w] = 1; + } + } + +// // A cross-flow connection is defined as a connection which has opposite +// // flow-direction to the well total flow +// V isInjPerf = (wops_.w2p * isInj); +// auto crossFlowConns = (connInjInx != isInjPerf); + +// bool allowCrossFlow = true; + +// if (not allowCrossFlow) { +// auto closedConns = crossFlowConns; +// for (int c = 0; c < nperf; ++c) { +// if (closedConns[c]) { +// Tw[c] = 0; +// } +// } +// connInjInx = !closedConns; +// } +// TODO: not allow for crossflow + + + V isInjInx = V::Zero(nperf); + V isNotInjInx = V::Zero(nperf); + for (int c = 0; c < nperf; ++c){ + if (connInjInx[c]) + isInjInx[c] = 1; + else + isNotInjInx[c] = 1; + } + + + // HANDLE FLOW INTO WELLBORE + + // compute phase volumerates standard conditions + std::vector cq_ps(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const ADB& wellcell_mob = subset ( rq_[phase].mob, well_cells); + const ADB cq_p = -(isNotInjInx * Tw) * (wellcell_mob * drawdown); + cq_ps[phase] = subset(rq_[phase].b,well_cells) * cq_p; + } + if (active_[Oil] && active_[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + ADB cq_psOil = cq_ps[oilpos]; + ADB cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += subset(state.rs,well_cells) * cq_psOil; + cq_ps[oilpos] += subset(state.rv,well_cells) * cq_psGas; + } + + // phase rates at std. condtions + std::vector q_ps(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + q_ps[phase] = wops_.w2p * cq_ps[phase]; + } + + // total rates at std + ADB qt_s = ADB::constant(V::Zero(nw), state.bhp.blockPattern()); + for (int phase = 0; phase < np; ++phase) { + qt_s += subset(state.qs, Span(nw, 1, phase*nw)); + } + + // compute avg. and total wellbore phase volumetric rates at std. conds + const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); + std::vector wbq(np, ADB::null()); + ADB wbqt = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); + for (int phase = 0; phase < np; ++phase) { + const int pos = pu.phase_pos[phase]; + wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase]; + wbqt += wbq[phase]; + } + + // check for dead wells + V isNotDeadWells = V::Constant(nperf,1.0); + for (int c = 0; c < nperf; ++c){ + if (wbqt.value()[c] == 0) + isNotDeadWells[c] = 0; + } + // compute wellbore mixture at std conds + std::vector mix_s(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + mix_s[phase] = isNotDeadWells * wbq[phase]/wbqt; + } + + + // HANDLE FLOW OUT FROM WELLBORE + + // Total mobilities + ADB mt = subset(rq_[0].mob,well_cells); + for (int phase = 1; phase < np; ++phase) { + mt += subset(rq_[phase].mob,well_cells); + + } + + // injection connections total volumerates + ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown); + + + // compute volume ratio between connection at standard conditions + ADB volRat = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); + std::vector cmix_s(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + cmix_s[phase] = wops_.w2p * mix_s[phase]; + } + + ADB well_rv = subset(state.rv,well_cells); + ADB well_rs = subset(state.rs,well_cells); + ADB d = V::Constant(nperf,1.0) - well_rv * well_rs; + + for (int phase = 0; phase < np; ++phase) { + ADB tmp = cmix_s[phase]; + + if (phase == Oil && active_[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + tmp = tmp - subset(state.rv,well_cells) * cmix_s[gaspos] / d; + } + if (phase == Gas && active_[Oil]) { + const int oilpos = pu.phase_pos[Oil]; + tmp = tmp - subset(state.rs,well_cells) * cmix_s[oilpos] / d; + } + volRat += tmp / subset(rq_[phase].b,well_cells); + + } + + // injecting connections total volumerates at std cond + ADB cqt_is = cqt_i/volRat; + + + // connection phase volumerates at std cond + std::vector cq_s(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + cq_s[phase] = cq_ps[phase] + (wops_.w2p * mix_s[phase])*cqt_is; + } + + // Add well contributions to mass balance equations + for (int phase = 0; phase < np; ++phase) { + residual_.mass_balance[phase] += superset(cq_s[phase],well_cells,nc); + } + + // Add WELL EQUATIONS + ADB qs = state.qs; + for (int phase = 0; phase < np; ++phase) { + qs -= superset(wops_.w2p * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + + } + residual_.well_flux_eq = qs; + + } + + void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state){ + // Handling BHP and SURFACE_RATE wells. + + const int np = wells_.number_of_phases; + const int nw = wells_.number_of_wells; + + V bhp_targets(nw); + V rate_targets(nw); + M rate_distr(nw, np*nw); + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells_.ctrls[w]; + if (well_controls_get_current_type(wc) == BHP) { + bhp_targets[w] = well_controls_get_current_target(wc); + rate_targets[w] = -1e100; + } else if (well_controls_get_current_type( wc ) == SURFACE_RATE) { + bhp_targets[w] = -1e100; + rate_targets[w] = well_controls_get_current_target(wc); + { + const double * distr = well_controls_get_current_distr( wc ); + for (int phase = 0; phase < np; ++phase) { + rate_distr.insert(w, phase*nw + w) = distr[phase]; + } + } + } else { + OPM_THROW(std::runtime_error, "Can only handle BHP and SURFACE_RATE type controls."); + } + } + const ADB bhp_residual = state.bhp - bhp_targets; + const ADB rate_residual = rate_distr * state.qs - rate_targets; + // Choose bhp residual for positive bhp targets. + Selector bhp_selector(bhp_targets); + residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual); + // DUMP(residual_.well_eq); + } + + void FullyImplicitBlackoilSolver::addOldWellEq(const SolutionState& state) + { // -------- Well equation, and well contributions to the mass balance equations -------- // Contribution to mass balance will have to wait. @@ -790,40 +1013,11 @@ namespace { residual_.well_flux_eq = state.qs + well_rates_all; // DUMP(residual_.well_flux_eq); - // Handling BHP and SURFACE_RATE wells. - V bhp_targets(nw); - V rate_targets(nw); - M rate_distr(nw, np*nw); - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells_.ctrls[w]; - if (well_controls_get_current_type(wc) == BHP) { - bhp_targets[w] = well_controls_get_current_target(wc); - rate_targets[w] = -1e100; - } else if (well_controls_get_current_type( wc ) == SURFACE_RATE) { - bhp_targets[w] = -1e100; - rate_targets[w] = well_controls_get_current_target(wc); - { - const double * distr = well_controls_get_current_distr( wc ); - for (int phase = 0; phase < np; ++phase) { - rate_distr.insert(w, phase*nw + w) = distr[phase]; - } - } - } else { - OPM_THROW(std::runtime_error, "Can only handle BHP and SURFACE_RATE type controls."); - } - } - const ADB bhp_residual = bhp - bhp_targets; - const ADB rate_residual = rate_distr * state.qs - rate_targets; - // Choose bhp residual for positive bhp targets. - Selector bhp_selector(bhp_targets); - residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual); - // DUMP(residual_.well_eq); + } - - V FullyImplicitBlackoilSolver::solveJacobianSystem() const { const int np = fluid_.numPhases(); diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 65d888026..0243c3274 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -158,6 +158,15 @@ namespace Opm { computeAccum(const SolutionState& state, const int aix ); + void + addOldWellEq(const SolutionState& state); + + void + addWellControlEq(const SolutionState& state); + + void + addWellEq(const SolutionState& state); + void assemble(const V& dtpv, const BlackoilState& x , From d821afe11f96ef09b962f51296531aea1573b172 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 17 Mar 2014 10:14:45 +0100 Subject: [PATCH 03/18] Whitespace fix. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index c055377ad..1cc8fcd5d 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -684,6 +684,9 @@ namespace { } + + + void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state){ const int nc = grid_.number_of_cells; @@ -862,6 +865,10 @@ namespace { } + + + + void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state){ // Handling BHP and SURFACE_RATE wells. @@ -897,6 +904,10 @@ namespace { // DUMP(residual_.well_eq); } + + + + void FullyImplicitBlackoilSolver::addOldWellEq(const SolutionState& state) { // -------- Well equation, and well contributions to the mass balance equations -------- @@ -1018,6 +1029,8 @@ namespace { + + V FullyImplicitBlackoilSolver::solveJacobianSystem() const { const int np = fluid_.numPhases(); @@ -1088,7 +1101,7 @@ namespace { V isSg = V::Zero(nc,1); bool disgas = false; - bool vapoil = false; + bool vapoil = false; // this is a temporary hack to find if vapoil or disgas // is a active component. Should be given directly from From 737affb077ae45cc78645a38873407a21540d8a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 17 Mar 2014 10:23:40 +0100 Subject: [PATCH 04/18] Minor whitespace adjustments. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 1cc8fcd5d..899358ec1 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -679,16 +679,14 @@ namespace { addWellEq(state); addWellControlEq(state); - - } - void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state){ - + void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state) + { const int nc = grid_.number_of_cells; const int np = wells_.number_of_phases; const int nw = wells_.number_of_wells; @@ -813,7 +811,6 @@ namespace { // injection connections total volumerates ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown); - // compute volume ratio between connection at standard conditions ADB volRat = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); std::vector cmix_s(np, ADB::null()); @@ -843,7 +840,6 @@ namespace { // injecting connections total volumerates at std cond ADB cqt_is = cqt_i/volRat; - // connection phase volumerates at std cond std::vector cq_s(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { @@ -862,14 +858,14 @@ namespace { } residual_.well_flux_eq = qs; - } - void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state){ + void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state) + { // Handling BHP and SURFACE_RATE wells. const int np = wells_.number_of_phases; @@ -1023,8 +1019,6 @@ namespace { // Set the well flux equation residual_.well_flux_eq = state.qs + well_rates_all; // DUMP(residual_.well_flux_eq); - - } From e8ee80571727d32e8bf14fc32dfd395dc12be270 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 18 Mar 2014 08:48:34 +0100 Subject: [PATCH 05/18] Work in progress: use WellDensitySegmented class. This is work in progress since it is non-working: phaseRates() is not used correctly, and changes must be made to WellState or equivalent. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 77 +++++++++++++++++++- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 4 + 2 files changed, 77 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 899358ec1..bd2982238 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -19,11 +19,11 @@ #include - #include #include #include #include +#include #include #include @@ -231,6 +231,7 @@ namespace { { const SolutionState state = constantState(x, xw); computeAccum(state, 0); + computeWellConnectionPressures(state, xw); } const double atol = 1.0e-12; @@ -613,6 +614,75 @@ namespace { + void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state, + const WellState& xw) + { + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + const int nperf = wells_.well_connpos[wells_.number_of_wells]; + const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + // Compute b, rsmax, rvmax values for perforations. + const ADB perf_press = subset(state.pressure, well_cells); + std::vector perf_cond(nperf); + const std::vector& pc = phaseCondition(); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = pc[well_cells[perf]]; + } + const PhaseUsage& pu = fluid_.phaseUsage(); + DataBlock b(nperf, pu.num_phases); + std::vector rssat_perf(nperf, 0.0); + std::vector rvsat_perf(nperf, 0.0); + if (pu.phase_used[BlackoilPhases::Aqua]) { + const ADB bw = fluid_.bWat(perf_press, well_cells); + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value(); + } + if (pu.phase_used[BlackoilPhases::Liquid]) { + const ADB perf_rs = subset(state.rs, well_cells); + const ADB bo = fluid_.bOil(perf_press, perf_rs, perf_cond, well_cells); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value(); + const V rssat = fluidRsSat(perf_press.value(), well_cells); + rssat_perf.assign(rssat.data(), rssat.data() + nperf); + } + if (pu.phase_used[BlackoilPhases::Vapour]) { + const ADB perf_rv = subset(state.rv, well_cells); + const ADB bg = fluid_.bGas(perf_press, perf_rv, perf_cond, well_cells); + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value(); + const V rvsat = fluidRvSat(perf_press.value(), well_cells); + rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf); + } + // b is row major, so can just copy data. + std::vector b_perf(b.data(), b.data() + nperf * pu.num_phases); + // Extract well connection depths. + const V depth = Eigen::Map(grid_.cell_centroids, grid_.number_of_cells, grid_.dimensions).rightCols<1>(); + const V pdepth = subset(depth, well_cells); + std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); + // Surface density. + std::vector surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases); + // Gravity + double grav = 0.0; + const double* g = geo_.gravity(); + const int dim = grid_.dimensions; + if (g) { + // Guard against gravity in anything but last dimension. + for (int dd = 0; dd < dim - 1; ++dd) { + assert(g[dd] == 0.0); + } + grav = g[dim - 1]; + } + + // 2. Compute pressure deltas, and store the results. + std::vector cdp = WellDensitySegmented + ::computeConnectionPressureDelta(wells_, xw, fluid_.phaseUsage(), + b_perf, rssat_perf, rvsat_perf, perf_depth, + surf_dens, grav); + well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + } + + + + + void FullyImplicitBlackoilSolver:: assemble(const V& pvdt, @@ -695,9 +765,8 @@ namespace { V Tw = Eigen::Map(wells_.WI, nperf); const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); - // pressure drop from solution - double tmp[2] = {-1.1319e-09, -2.0926e-09 }; - V cdp = Eigen::Map(tmp, 2); + // pressure diffs computed already (once per step, not changing per iteration) + const V& cdp = well_perforation_pressure_diffs_; // Extract variables for perforation cell pressures // and corresponding perforation well pressures. diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 0243c3274..671169bb3 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -135,6 +135,7 @@ namespace Opm { std::vector rq_; std::vector phaseCondition_; + V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. // The mass_balance vector has one element for each active phase, // each of which has size equal to the number of cells. @@ -158,6 +159,9 @@ namespace Opm { computeAccum(const SolutionState& state, const int aix ); + void computeWellConnectionPressures(const SolutionState& state, + const WellState& xw); + void addOldWellEq(const SolutionState& state); From 33d58792dece9ca85b7ef7b5f926c8cbcda6ccca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 24 Mar 2014 15:28:46 +0100 Subject: [PATCH 06/18] Fix sign bug in well contribution to mass balance. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 4052b2417..1f52070c4 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -917,7 +917,7 @@ namespace { // Add well contributions to mass balance equations for (int phase = 0; phase < np; ++phase) { - residual_.mass_balance[phase] += superset(cq_s[phase],well_cells,nc); + residual_.mass_balance[phase] -= superset(cq_s[phase],well_cells,nc); } // Add WELL EQUATIONS From bad64de4f2656b7f974e3085a36d65c042a7f5be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Mar 2014 11:11:19 +0100 Subject: [PATCH 07/18] Add currentControls() field to well state class. This is done since the solver will need to be able to switch well controls during Newton iterations. The current control specified in the Wells struct will be used as default and initial value for currentControls(). --- opm/autodiff/WellStateFullyImplicitBlackoil.hpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index f2430b1b8..ea74f4e8a 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -70,6 +70,14 @@ namespace Opm } } } + + // Initialize current_controls_. + // The controls set in the Wells object are treated as defaults, + // and also used for initial values. + current_controls_.resize(nw); + for (int w = 0; w < nw; ++w) { + current_controls_[w] = well_controls_get_current(wells->ctrls[w]); + } } /// One bhp pressure per well. @@ -84,6 +92,10 @@ namespace Opm std::vector& perfPhaseRates() { return perfphaserates_; } const std::vector& perfPhaseRates() const { return perfphaserates_; } + /// One current control per well. + std::vector& currentControls() { return current_controls_; } + const std::vector& currentControls() const { return current_controls_; } + /// For interfacing with functions that take a WellState. const WellState& basicWellState() const { @@ -93,6 +105,7 @@ namespace Opm private: WellState basic_well_state_; std::vector perfphaserates_; + std::vector current_controls_; }; } // namespace Opm From 20ecd37709bd5cc11817b71a57f7bd517211d8be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Mar 2014 11:14:11 +0100 Subject: [PATCH 08/18] Use currentControls() field of well state class. Now the well state is consulted for the controls to use when assembling the well control equations, instead of the (immutable) wells_ struct. Also, added dummy implementation of updateWellControls(). --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 31 +++++++++++++++----- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 7 +++-- 2 files changed, 28 insertions(+), 10 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 1f52070c4..55d460733 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -687,7 +687,7 @@ namespace { FullyImplicitBlackoilSolver:: assemble(const V& pvdt, const BlackoilState& x , - const WellStateFullyImplicitBlackoil& xw ) + WellStateFullyImplicitBlackoil& xw ) { // Create the primary variables. const SolutionState state = variableState(x, xw); @@ -748,7 +748,8 @@ namespace { } addWellEq(state); - addWellControlEq(state); + updateWellControls(xw); + addWellControlEq(state, xw); } @@ -933,7 +934,17 @@ namespace { - void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state) + void FullyImplicitBlackoilSolver::updateWellControls(WellStateFullyImplicitBlackoil& /*xw*/) const + { + // Do stuff. + } + + + + + + void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state, + const WellStateFullyImplicitBlackoil& xw) { // Handling BHP and SURFACE_RATE wells. @@ -945,14 +956,18 @@ namespace { M rate_distr(nw, np*nw); for (int w = 0; w < nw; ++w) { const WellControls* wc = wells_.ctrls[w]; - if (well_controls_get_current_type(wc) == BHP) { - bhp_targets[w] = well_controls_get_current_target(wc); + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = xw.currentControls()[w]; + if (well_controls_iget_type(wc, current) == BHP) { + bhp_targets[w] = well_controls_iget_target(wc, current); rate_targets[w] = -1e100; - } else if (well_controls_get_current_type( wc ) == SURFACE_RATE) { + } else if (well_controls_iget_type(wc, current) == SURFACE_RATE) { bhp_targets[w] = -1e100; - rate_targets[w] = well_controls_get_current_target(wc); + rate_targets[w] = well_controls_iget_target(wc, current); { - const double * distr = well_controls_get_current_distr( wc ); + const double * distr = well_controls_iget_distr(wc, current); for (int phase = 0; phase < np; ++phase) { rate_distr.insert(w, phase*nw + w) = distr[phase]; } diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index c9c8ac487..4cd7f1528 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -166,15 +166,18 @@ namespace Opm { addOldWellEq(const SolutionState& state); void - addWellControlEq(const SolutionState& state); + addWellControlEq(const SolutionState& state, + const WellStateFullyImplicitBlackoil& xw); void addWellEq(const SolutionState& state); + void updateWellControls(WellStateFullyImplicitBlackoil& xw) const; + void assemble(const V& dtpv, const BlackoilState& x, - const WellStateFullyImplicitBlackoil& xw); + WellStateFullyImplicitBlackoil& xw); V solveJacobianSystem() const; From 4a22c560c7b987e94ecfa2da0e06ce8fe399a144 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Mar 2014 14:31:06 +0100 Subject: [PATCH 09/18] Implemented updateWellControls(). --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 122 ++++++++++++++++++- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 7 +- 2 files changed, 122 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 55d460733..f1931765c 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -747,8 +747,7 @@ namespace { } - addWellEq(state); - updateWellControls(xw); + addWellEq(state, xw); // Note: this can change xw (switching well controls). addWellControlEq(state, xw); } @@ -756,7 +755,8 @@ namespace { - void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state) + void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state, + WellStateFullyImplicitBlackoil& xw) { const int nc = grid_.number_of_cells; const int np = wells_.number_of_phases; @@ -928,15 +928,127 @@ namespace { } residual_.well_flux_eq = qs; + + // Updating the well controls is done from here, since we have + // access to the necessary bhp and rate data in this method. + updateWellControls(state.bhp, state.qs, xw); } - void FullyImplicitBlackoilSolver::updateWellControls(WellStateFullyImplicitBlackoil& /*xw*/) const + namespace { - // Do stuff. + double rateToCompare(const ADB& well_phase_flow_rate, + const int well, + const int num_phases, + const double* distr) + { + double rate = 0.0; + for (int phase = 0; phase < num_phases; ++phase) { + rate += well_phase_flow_rate.value()[well*num_phases + phase] * distr[phase]; + } + return rate; + } + + bool constraintBroken(const ADB& bhp, + const ADB& well_phase_flow_rate, + const int well, + const int num_phases, + const WellType& well_type, + const WellControls* wc, + const int ctrl_index) + { + const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index); + const double target = well_controls_iget_target(wc, ctrl_index); + const double* distr = well_controls_iget_distr(wc, ctrl_index); + switch (well_type) { + case INJECTOR: + switch (ctrl_type) { + case BHP: + return bhp.value()[well] > target; + case SURFACE_RATE: + return rateToCompare(well_phase_flow_rate, well, num_phases, distr) > target; + case RESERVOIR_RATE: + // Intentional fall-through. + default: + OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls."); + } + break; + case PRODUCER: + switch (ctrl_type) { + case BHP: + return bhp.value()[well] < target; + case SURFACE_RATE: + // Note that the rates compared below are negative, + // so breaking the constraints means: too high flow rate + // (as for injection). + return rateToCompare(well_phase_flow_rate, well, num_phases, distr) < target; + case RESERVOIR_RATE: + // Intentional fall-through. + default: + OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls."); + } + break; + default: + OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells."); + } + } + } // anonymous namespace + + + + + + void FullyImplicitBlackoilSolver::updateWellControls(const ADB& bhp, + const ADB& well_phase_flow_rate, + WellStateFullyImplicitBlackoil& xw) const + { + std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" }; + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. + const int np = wells_.number_of_phases; + const int nw = wells_.number_of_wells; + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells_.ctrls[w]; + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = xw.currentControls()[w]; + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (well_controls_iget_type(wc, ctrl_index) == RESERVOIR_RATE) { + // We cannot handle this yet. +#ifdef OPM_VERBOSE + std::cout << "Warning: a RESERVOIR_RATE well control exists for well " + << wells_.name[w] << ", but will never be checked." << std::endl; +#endif + continue; + } + if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells_.type[w], wc, ctrl_index)) { + // ctrl_index will be the index of the broken constraint after the loop. + break; + } + } + if (ctrl_index != nwc) { + // Constraint number ctrl_index was broken, switch to it. + std::cout << "Switching control mode for well " << wells_.name[w] + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)]; + xw.currentControls()[w] = ctrl_index; + } + } } diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 4cd7f1528..15c27e35a 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -170,9 +170,12 @@ namespace Opm { const WellStateFullyImplicitBlackoil& xw); void - addWellEq(const SolutionState& state); + addWellEq(const SolutionState& state, + WellStateFullyImplicitBlackoil& xw); - void updateWellControls(WellStateFullyImplicitBlackoil& xw) const; + void updateWellControls(const ADB& bhp, + const ADB& well_phase_flow_rate, + WellStateFullyImplicitBlackoil& xw) const; void assemble(const V& dtpv, From cb50728ac27ad884a1f5d2b0cd3ac078604ebbee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Mar 2014 17:41:41 +0100 Subject: [PATCH 10/18] Add end-of-line to switching message. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index f1931765c..76951061c 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -1045,7 +1045,7 @@ namespace { // Constraint number ctrl_index was broken, switch to it. std::cout << "Switching control mode for well " << wells_.name[w] << " from " << modestring[well_controls_iget_type(wc, current)] - << " to " << modestring[well_controls_iget_type(wc, ctrl_index)]; + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; xw.currentControls()[w] = ctrl_index; } } From dcb59a2fab361f027d34499a033b8e3d95141623 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Mar 2014 17:42:14 +0100 Subject: [PATCH 11/18] Fix bug in access order for state.qs. State.qs is called well_phase_flow_rate in this function. It is stored with wells running fastest, not phases. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 76951061c..439c5b05b 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -945,9 +945,12 @@ namespace { const int num_phases, const double* distr) { + const int num_wells = well_phase_flow_rate.size() / num_phases; double rate = 0.0; for (int phase = 0; phase < num_phases; ++phase) { - rate += well_phase_flow_rate.value()[well*num_phases + phase] * distr[phase]; + // Important: well_phase_flow_rate is ordered with all rates for first + // phase coming first, then all for second phase etc. + rate += well_phase_flow_rate.value()[well + phase*num_wells] * distr[phase]; } return rate; } From 41fd953ba20f9049296e0f4bc2c409aeb280ae28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Mar 2014 18:44:37 +0100 Subject: [PATCH 12/18] Fix bugs relating to wells versus perforations. Two kinds of bugs: correct usage of wops_.w2p vs. wops_.p2w, and correct sizes for some variables (nw vs. nperf). --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 439c5b05b..6c36c98a5 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -837,7 +837,7 @@ namespace { // phase rates at std. condtions std::vector q_ps(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { - q_ps[phase] = wops_.w2p * cq_ps[phase]; + q_ps[phase] = wops_.p2w * cq_ps[phase]; } // total rates at std @@ -849,7 +849,7 @@ namespace { // compute avg. and total wellbore phase volumetric rates at std. conds const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); std::vector wbq(np, ADB::null()); - ADB wbqt = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); + ADB wbqt = ADB::constant(V::Zero(nw), state.pressure.blockPattern()); for (int phase = 0; phase < np; ++phase) { const int pos = pu.phase_pos[phase]; wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase]; @@ -857,8 +857,8 @@ namespace { } // check for dead wells - V isNotDeadWells = V::Constant(nperf,1.0); - for (int c = 0; c < nperf; ++c){ + V isNotDeadWells = V::Constant(nw,1.0); + for (int c = 0; c < nw; ++c){ if (wbqt.value()[c] == 0) isNotDeadWells[c] = 0; } @@ -924,7 +924,7 @@ namespace { // Add WELL EQUATIONS ADB qs = state.qs; for (int phase = 0; phase < np; ++phase) { - qs -= superset(wops_.w2p * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); } residual_.well_flux_eq = qs; From 3bc57780aaa53889efdbc9751c6ada3ece515160 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Mar 2014 11:44:49 +0100 Subject: [PATCH 13/18] Fix bug in -= operator. --- opm/autodiff/AutoDiffBlock.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index a98c014be..7255b8dd1 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -198,7 +198,11 @@ namespace Opm AutoDiffBlock& operator-=(const AutoDiffBlock& rhs) { if (jac_.empty()) { - jac_ = rhs.jac_; + const int num_blocks = rhs.numBlocks(); + jac_.resize(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jac_[block] = -rhs.jac_[block]; + } } else if (!rhs.jac_.empty()) { assert (numBlocks() == rhs.numBlocks()); assert (value().size() == rhs.value().size()); From 37a23825c4cc90723afb59380b3d9b6cc3bc4a0b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Mar 2014 11:48:13 +0100 Subject: [PATCH 14/18] Removed dead code in test_block.cpp. --- tests/test_block.cpp | 51 -------------------------------------------- 1 file changed, 51 deletions(-) diff --git a/tests/test_block.cpp b/tests/test_block.cpp index 6a8936e3c..c0dceefe0 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -190,54 +190,3 @@ BOOST_AUTO_TEST_CASE(Addition) BOOST_CHECK(*j3b == ADB::M((*j1b) * 2)); } } - -#if 0 -#include - -int main() -try -{ - typedef AutoDiffBlock ADB; - std::vector blocksizes = { 3, 1, 2 }; - int num_blocks = blocksizes.size(); - ADB::V v1(3); - v1 << 0.2, 1.2, 13.4; - ADB::V v2(3); - v2 << 1.0, 2.2, 3.4; - enum { FirstVar = 0, SecondVar = 1, ThirdVar = 2 }; - ADB a = ADB::constant(v1, blocksizes); - ADB x = ADB::variable(FirstVar, v2, blocksizes); - std::vector jacs(num_blocks); - for (int i = 0; i < num_blocks; ++i) { - jacs[i] = ADB::M(blocksizes[FirstVar], blocksizes[i]); - jacs[i].insert(0,0) = -1.0; - } - ADB f = ADB::function(v2, jacs); - - ADB xpx = x + x; - std::cout << xpx; - ADB xpxpa = x + x + a; - std::cout << xpxpa; - - - std::cout << xpxpa - xpx; - - ADB sqx = x * x; - - std::cout << sqx; - - ADB sqxdx = sqx / x; - - std::cout << sqxdx; - - ADB::M m(2,3); - m.insert(0,0) = 4; - m.insert(0,1) = 3; - m.insert(1,1) = 1; - std::cout << m*sqx; -} -catch (const std::exception &e) { - std::cerr << "Program threw an exception: " << e.what() << "\n"; - throw; -} -#endif From 60ff38adc6ea2885dc4bb1844e9b717e2c4f229e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Mar 2014 12:03:44 +0100 Subject: [PATCH 15/18] Added test for += and -= operators. --- tests/test_block.cpp | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/tests/test_block.cpp b/tests/test_block.cpp index c0dceefe0..f92a9c85f 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -190,3 +190,32 @@ BOOST_AUTO_TEST_CASE(Addition) BOOST_CHECK(*j3b == ADB::M((*j1b) * 2)); } } + +BOOST_AUTO_TEST_CASE(AssignAddSubtractOperators) +{ + typedef AutoDiffBlock ADB; + + ADB::V vx(3); + vx << 0.2, 1.2, 13.4; + + ADB::V vy(3); + vy << 1.0, 2.2, 3.4; + + std::vector vals{ vx, vy }; + std::vector vars = ADB::variables(vals); + + const ADB x = vars[0]; + const ADB y = vars[1]; + + ADB z = x; + z += y; + ADB sum = x + y; + const double tolerance = 1e-14; + BOOST_CHECK(z.value().isApprox(sum.value(), tolerance)); + BOOST_CHECK(z.derivative()[0].isApprox(sum.derivative()[0], tolerance)); + BOOST_CHECK(z.derivative()[1].isApprox(sum.derivative()[1], tolerance)); + z -= y; + BOOST_CHECK(z.value().isApprox(x.value(), tolerance)); + BOOST_CHECK(z.derivative()[0].isApprox(x.derivative()[0], tolerance)); + BOOST_CHECK(z.derivative()[1].isApprox(x.derivative()[1], tolerance)); +} From 347016a0667d01f420f5f34d6034067ce3e318fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Mar 2014 12:25:11 +0100 Subject: [PATCH 16/18] Add more checks to cover bug that was in operator -=. --- tests/test_block.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/test_block.cpp b/tests/test_block.cpp index f92a9c85f..102237924 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -195,6 +195,7 @@ BOOST_AUTO_TEST_CASE(AssignAddSubtractOperators) { typedef AutoDiffBlock ADB; + // Basic testing of += and -=. ADB::V vx(3); vx << 0.2, 1.2, 13.4; @@ -218,4 +219,17 @@ BOOST_AUTO_TEST_CASE(AssignAddSubtractOperators) BOOST_CHECK(z.value().isApprox(x.value(), tolerance)); BOOST_CHECK(z.derivative()[0].isApprox(x.derivative()[0], tolerance)); BOOST_CHECK(z.derivative()[1].isApprox(x.derivative()[1], tolerance)); + + // Testing the case when the left hand side has empty() jacobian. + ADB yconst = ADB::constant(vy); + z = yconst; + z -= x; + ADB diff = yconst - x; + BOOST_CHECK(z.value().isApprox(diff.value(), tolerance)); + BOOST_CHECK(z.derivative()[0].isApprox(diff.derivative()[0], tolerance)); + BOOST_CHECK(z.derivative()[1].isApprox(diff.derivative()[1], tolerance)); + z += x; + BOOST_CHECK(z.value().isApprox(yconst.value(), tolerance)); + BOOST_CHECK(z.derivative()[0].isApprox(Eigen::Matrix::Zero())); + BOOST_CHECK(z.derivative()[1].isApprox(Eigen::Matrix::Zero())); } From 091484182350efbec4405df9096216960bef0a2f Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 27 Mar 2014 11:15:55 +0100 Subject: [PATCH 17/18] A small fixing to make the sim_fibo_ad using the new timer. --- examples/sim_fibo_ad.cpp | 6 ++---- opm/autodiff/SimulatorFullyImplicitBlackoil.cpp | 1 + 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index 196ab009d..76777d784 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -205,7 +205,7 @@ try std::shared_ptr eclipseState(new EclipseState(newParserDeck)); // initialize variables - simtimer.init(timeMap, /*beginReportStepIdx=*/0, /*endReportStepIdx=*/0); + simtimer.init(timeMap); SimulatorReport fullReport; for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) { @@ -228,9 +228,7 @@ try well_state.init(wells.c_wells(), state); } - simtimer.init(timeMap, - /*beginReportStepIdx=*/reportStepIdx, - /*endReportStepIdx=*/reportStepIdx + 1); + simtimer.setCurrentStepNum(reportStepIdx); if (reportStepIdx == 0) outputWriter.writeInit(simtimer, state, well_state.basicWellState()); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp index 9416cb8ab..9d9b022e5 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -412,6 +412,7 @@ namespace Opm // advance to next timestep before reporting at this location ++timer; + break; // this is a temporary measure } total_timer.stop(); From 4b53f8868509761a27f8ad56810ff5d564fc13c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 28 Mar 2014 14:31:30 +0100 Subject: [PATCH 18/18] Fix timer usage. This makes the simulator produce proper summary output again. --- examples/sim_fibo_ad.cpp | 6 +++++- opm/autodiff/SimulatorFullyImplicitBlackoil.cpp | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index 76777d784..a13d2d29f 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -230,8 +230,10 @@ try simtimer.setCurrentStepNum(reportStepIdx); - if (reportStepIdx == 0) + if (reportStepIdx == 0) { outputWriter.writeInit(simtimer, state, well_state.basicWellState()); + outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState()); + } // Create and run simulator. SimulatorFullyImplicitBlackoil simulator(param, @@ -243,6 +245,8 @@ try grav); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); + ++simtimer; + outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState()); fullReport += episodeReport; } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp index 9d9b022e5..a36ae887e 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -411,7 +411,7 @@ namespace Opm } // advance to next timestep before reporting at this location - ++timer; + // ++timer; // Commented out since this has temporarily moved to the main() function. break; // this is a temporary measure }