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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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/28] 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 9237aa2443c27c7437d609e4f4a51031704d6886 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 25 Mar 2014 11:52:13 +0100 Subject: [PATCH 09/28] FullyImplicitBlackoilSolver: make constantState() and variableState() consistent With this, constantState() just calls variableState() and throws away the derivatives. This might be a bit slower than necessary, but it makes these two methods automatically consistent and constantState() is only called once per timestep anyway... thanks to @atgeirr for the suggestion! --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 95 +++----------------- 1 file changed, 12 insertions(+), 83 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index b0b53c44e..8b8cc7f8f 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -333,90 +333,19 @@ namespace { FullyImplicitBlackoilSolver::constantState(const BlackoilState& x, const WellStateFullyImplicitBlackoil& xw) { - const int nc = grid_.number_of_cells; - const int np = x.numPhases(); + auto state = variableState(x, xw); - // The block pattern assumes the following primary variables: - // pressure - // water saturation (if water present) - // gas saturation, Rv (vapor oil/gas ratio) or Rs (solution gas/oil ratio) depending on hydrocarbon state - // Gas only (undersaturated gas): Rv - // Gas and oil: Sg - // Oil only (undersaturated oil): Rs - // well rates per active phase and well - // well bottom-hole pressure - // Note that oil is assumed to always be present, but is never - // a primary variable. - assert(active_[ Oil ]); - std::vector bpat(np, nc); - bpat.push_back(xw.bhp().size() * np); - bpat.push_back(xw.bhp().size()); - - SolutionState state(np); - - // Pressure. - assert (not x.pressure().empty()); - const V p = Eigen::Map(& x.pressure()[0], nc, 1); - state.pressure = ADB::constant(p, bpat); - - // Saturation. - assert (not x.saturation().empty()); - const DataBlock s = Eigen::Map(& x.saturation()[0], nc, np); - const Opm::PhaseUsage pu = fluid_.phaseUsage(); - { - V so = V::Ones(nc, 1); - if (active_[ Water ]) { - const int pos = pu.phase_pos[ Water ]; - const V sw = s.col(pos); - so -= sw; - - state.saturation[pos] = ADB::constant(sw, bpat); - } - if (active_[ Gas ]) { - const int pos = pu.phase_pos[ Gas ]; - const V sg = s.col(pos); - so -= sg; - - state.saturation[pos] = ADB::constant(sg, bpat); - } - if (active_[ Oil ]) { - const int pos = pu.phase_pos[ Oil ]; - state.saturation[pos] = ADB::constant(so, bpat); - } - } - - // Solution Gas-oil ratio (rs). - if (active_[ Oil ] && active_[ Gas ]) { - const V rs = Eigen::Map(& x.gasoilratio()[0], x.gasoilratio().size()); - state.rs = ADB::constant(rs, bpat); - } else { - const V Rs = V::Zero(nc, 1); - state.rs = ADB::constant(Rs, bpat); - } - - // Vapor Oil-gas ratio (rv). - if (active_[ Oil ] && active_[ Gas ]) { - const V rv = Eigen::Map(& x.rv()[0], x.rv().size()); - state.rv = ADB::constant(rv, bpat); - } else { - const V rv = V::Zero(nc, 1); - state.rv = ADB::constant(rv, bpat); - } - - // Well rates. - assert (not xw.wellRates().empty()); - // Need to reshuffle well rates, from ordered by wells, then phase, - // to ordered by phase, then wells. - 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); - - // Well bottom-hole pressure. - assert (not xw.bhp().empty()); - const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); - state.bhp = ADB::constant(bhp, bpat); + // HACK: throw away the derivatives. this may not be the most + // performant way to do things, but it will make the state + // automatically consistent with variableState() (and doing + // things automatically is all the rage in this module ;) + state.pressure = ADB::constant(state.pressure.value()); + state.rs = ADB::constant(state.rs.value()); + state.rv = ADB::constant(state.rv.value()); + for (int phaseIdx= 0; phaseIdx < x.numPhases(); ++ phaseIdx) + state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value()); + state.qs = ADB::constant(state.qs.value()); + state.bhp = ADB::constant(state.bhp.value()); return state; } From 168b4379d5ab61e8b27b7f756ba4c49119be912b Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 20 Mar 2014 17:03:30 +0100 Subject: [PATCH 10/28] don't pass the eclipsewriter to the simulator anymore i.e., the simulator does not deal with any output operations anymore. This makes sense because report steps are handled by Opm::TimeMap and the constructor for the simulator already needs more arguments than appropriate even without this... --- examples/sim_fibo_ad.cpp | 10 ++------- .../SimulatorFullyImplicitBlackoil.cpp | 22 +++++-------------- .../SimulatorFullyImplicitBlackoil.hpp | 4 +--- 3 files changed, 8 insertions(+), 28 deletions(-) diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index 4f30a6ef0..196ab009d 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -242,16 +242,11 @@ try rock_comp->isActive() ? rock_comp.get() : 0, wells, linsolver, - grav, - outputWriter); + grav); SimulatorReport episodeReport = simulator.run(simtimer, state, well_state); outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState()); fullReport += episodeReport; - - if (output) { - episodeReport.reportParam(outStream); - } } std::cout << "\n\n================ End of simulation ===============\n\n"; @@ -318,8 +313,7 @@ try rock_comp->isActive() ? rock_comp.get() : 0, wells, linsolver, - grav, - outputWriter); + grav); outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState()); if (epoch == 0) { diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp index aae39d9ca..9416cb8ab 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.cpp @@ -38,7 +38,6 @@ #include #include #include -#include #include #include #include @@ -72,8 +71,7 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, LinearSolverInterface& linsolver, - const double* gravity, - EclipseWriter &writer); + const double* gravity); SimulatorReport run(SimulatorTimer& timer, BlackoilState& state, @@ -102,7 +100,6 @@ namespace Opm FullyImplicitBlackoilSolver solver_; // Misc. data std::vector allcells_; - EclipseWriter &eclipseWriter_; }; @@ -114,11 +111,10 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, LinearSolverInterface& linsolver, - const double* gravity, - EclipseWriter &eclipseWriter) + const double* gravity) { - pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity, eclipseWriter)); + pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, linsolver, gravity)); } @@ -262,8 +258,7 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, LinearSolverInterface& linsolver, - const double* gravity, - EclipseWriter &eclipseWriter) + const double* gravity) : grid_(grid), props_(props), rock_comp_props_(rock_comp_props), @@ -271,9 +266,7 @@ namespace Opm wells_(wells_manager.c_wells()), gravity_(gravity), geo_(grid_, props_, gravity_), - solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver), - eclipseWriter_(eclipseWriter) - + solver_(grid_, props_, geo_, rock_comp_props, *wells_manager.c_wells(), linsolver) /* param.getDefault("nl_pressure_residual_tolerance", 0.0), param.getDefault("nl_pressure_change_tolerance", 1.0), param.getDefault("nl_pressure_maxiter", 10), @@ -419,11 +412,6 @@ namespace Opm // advance to next timestep before reporting at this location ++timer; - - // write an output file for later inspection - if (output_) { - eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState()); - } } total_timer.stop(); diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp index 26b91b4e3..07e80fda2 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoil.hpp @@ -31,7 +31,6 @@ namespace Opm { namespace parameter { class ParameterGroup; } class BlackoilPropsAdInterface; - class EclipseWriter; class RockCompressibility; class WellsManager; class LinearSolverInterface; @@ -72,8 +71,7 @@ namespace Opm const RockCompressibility* rock_comp_props, WellsManager& wells_manager, LinearSolverInterface& linsolver, - const double* gravity, - EclipseWriter &writer); + const double* gravity); /// Run the simulation. /// This will run succesive timesteps until timer.done() is true. It will 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 11/28] 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 12/28] 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 13/28] 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 14/28] 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 2f23e72beafbc39b7de7c772784b29ac30233b71 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Tue, 25 Mar 2014 19:36:42 +0100 Subject: [PATCH 15/28] Removed calls to WellsManager() constructor which uses old EclipseGridParser. --- examples/sim_2p_comp_ad.cpp | 7 ++++++- examples/sim_2p_incomp_ad.cpp | 9 ++++++++- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/examples/sim_2p_comp_ad.cpp b/examples/sim_2p_comp_ad.cpp index be7c4c205..c7dfe8bfb 100644 --- a/examples/sim_2p_comp_ad.cpp +++ b/examples/sim_2p_comp_ad.cpp @@ -48,6 +48,8 @@ #include #include +#include + #include #include @@ -94,6 +96,7 @@ try boost::scoped_ptr grid; boost::scoped_ptr props; boost::scoped_ptr rock_comp; + EclipseStateConstPtr eclipseState; BlackoilState state; // bool check_well_controls = false; // int max_well_control_iterations = 0; @@ -103,6 +106,8 @@ try deck.reset(new EclipseGridParser(deck_filename)); Opm::ParserPtr newParser(new Opm::Parser() ); Opm::DeckConstPtr newParserDeck = newParser->parseFile( deck_filename ); + eclipseState.reset( new EclipseState(newParserDeck )); + // Grid init grid.reset(new GridManager(newParserDeck)); @@ -255,7 +260,7 @@ try << simtimer.numSteps() - step << ")\n\n" << std::flush; // Create new wells, well_state - WellsManager wells(*deck, *grid->c_grid(), props->permeability()); + WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability()); // @@@ HACK: we should really make a new well state and // properly transfer old well state to it every epoch, // since number of wells may change etc. diff --git a/examples/sim_2p_incomp_ad.cpp b/examples/sim_2p_incomp_ad.cpp index 7469bdca0..e0c417c14 100644 --- a/examples/sim_2p_incomp_ad.cpp +++ b/examples/sim_2p_incomp_ad.cpp @@ -46,6 +46,9 @@ #include #include +#include +#include + #include #include @@ -100,12 +103,16 @@ try boost::scoped_ptr grid; boost::scoped_ptr props; boost::scoped_ptr rock_comp; + EclipseStateConstPtr eclipseState; TwophaseState state; // bool check_well_controls = false; // int max_well_control_iterations = 0; double gravity[3] = { 0.0 }; if (use_deck) { std::string deck_filename = param.get("deck_filename"); + ParserPtr parser(new Opm::Parser()); + eclipseState.reset( new EclipseState(parser->parseFile(deck_filename))); + deck.reset(new EclipseGridParser(deck_filename)); // Grid init grid.reset(new GridManager(*deck)); @@ -262,7 +269,7 @@ try << simtimer.numSteps() - step << ")\n\n" << std::flush; // Create new wells, well_state - WellsManager wells(*deck, *grid->c_grid(), props->permeability()); + WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability()); // @@@ HACK: we should really make a new well state and // properly transfer old well state to it every epoch, // since number of wells may change etc. 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 16/28] 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 17/28] 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 18/28] 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 19/28] 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 2f84156b91e359cc174632dc3eaf5fbc66e87944 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 26 Mar 2014 11:53:48 +0100 Subject: [PATCH 20/28] FullyImplicitBlackoilSolver: Throw an exception uppon encountering a non-finite residual --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index b0b53c44e..fd51f3d36 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -31,6 +31,7 @@ #include #include #include +#include #include #include @@ -1272,7 +1273,12 @@ namespace { e = residual_.mass_balance.end(); b != e; ++b) { - r = std::max(r, (*b).value().matrix().norm()); + const double rowResid = (*b).value().matrix().norm(); + if (!std::isfinite(rowResid)) { + OPM_THROW(Opm::NumericalProblem, + "Encountered a non-finite residual"); + } + r = std::max(r, rowResid); } r = std::max(r, residual_.well_flux_eq.value().matrix().norm()); r = std::max(r, residual_.well_eq.value().matrix().norm()); From 8cd93bf70a41ef31d759afb12e458448bf2ff747 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 26 Mar 2014 13:28:49 +0100 Subject: [PATCH 21/28] FullyImplicitBlackoilSolver: cleanup the residualNorm() method --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 21 ++++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index fd51f3d36..9a0efed21 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -1267,23 +1267,22 @@ namespace { double FullyImplicitBlackoilSolver::residualNorm() const { - double r = 0; - for (std::vector::const_iterator - b = residual_.mass_balance.begin(), - e = residual_.mass_balance.end(); - b != e; ++b) + double globalNorm = 0; + std::vector::const_iterator quantityIt = residual_.mass_balance.begin(); + const std::vector::const_iterator endQuantityIt = residual_.mass_balance.end(); + for (; quantityIt != endQuantityIt; ++quantityIt) { - const double rowResid = (*b).value().matrix().norm(); - if (!std::isfinite(rowResid)) { + const double quantityResid = (*quantityIt).value().matrix().norm(); + if (!std::isfinite(quantityResid)) { OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual"); } - r = std::max(r, rowResid); + globalNorm = std::max(globalNorm, quantityResid); } - r = std::max(r, residual_.well_flux_eq.value().matrix().norm()); - r = std::max(r, residual_.well_eq.value().matrix().norm()); + globalNorm = std::max(globalNorm, residual_.well_flux_eq.value().matrix().norm()); + globalNorm = std::max(globalNorm, residual_.well_eq.value().matrix().norm()); - return r; + return globalNorm; } From 091484182350efbec4405df9096216960bef0a2f Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 27 Mar 2014 11:15:55 +0100 Subject: [PATCH 22/28] 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 23/28] 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 } From 5e7064b337386d38b2cd68581631889e47e931ee Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 31 Mar 2014 16:39:17 +0200 Subject: [PATCH 24/28] sim_fibo_ad: remove compatibility code for old parser this makes the simulator quite a bit more maintainable: setting USE_NEW_PARSER to 0 did not even compile after the the constructor for the wells manager which took the old deck was removed last week. Since, according to Atgeirr, SPE-1 is now producing exactly the same results as before, it also does no longer make too much sense to keep that code on life support... --- examples/sim_fibo_ad.cpp | 114 +-------------------------------------- 1 file changed, 2 insertions(+), 112 deletions(-) diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index a13d2d29f..dda83487b 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -99,51 +99,25 @@ try double gravity[3] = { 0.0 }; std::string deck_filename = param.get("deck_filename"); -#define USE_NEW_PARSER 1 -#if USE_NEW_PARSER Opm::ParserPtr newParser(new Opm::Parser() ); Opm::DeckConstPtr newParserDeck = newParser->parseFile( deck_filename ); -#else - std::shared_ptr deck; - deck.reset(new EclipseGridParser(deck_filename)); -#endif // Grid init -#if USE_NEW_PARSER grid.reset(new GridManager(newParserDeck)); -#else - grid.reset(new GridManager(*deck)); -#endif - -#if USE_NEW_PARSER Opm::EclipseWriter outputWriter(param, newParserDeck, share_obj(*grid->c_grid())); -#else - Opm::EclipseWriter outputWriter(param, deck, share_obj(*grid->c_grid())); -#endif // Rock and fluid init -#if USE_NEW_PARSER props.reset(new BlackoilPropertiesFromDeck(newParserDeck, *grid->c_grid(), param)); new_props.reset(new BlackoilPropsAdFromDeck(newParserDeck, *grid->c_grid())); -#else - props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid(), param)); - new_props.reset(new BlackoilPropsAdFromDeck(*deck, *grid->c_grid())); -#endif // check_well_controls = param.getDefault("check_well_controls", false); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Rock compressibility. -#if USE_NEW_PARSER rock_comp.reset(new RockCompressibility(newParserDeck)); -#else - rock_comp.reset(new RockCompressibility(*deck)); -#endif + // Gravity. -#if USE_NEW_PARSER gravity[2] = newParserDeck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; -#else - gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; -#endif + // Init state variables (saturation and pressure). if (param.has("init_saturation")) { initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); @@ -159,11 +133,7 @@ try } } } else { -#if USE_NEW_PARSER initBlackoilStateFromDeck(*grid->c_grid(), *props, newParserDeck, gravity[2], state); -#else - initBlackoilStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); -#endif } bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); @@ -195,7 +165,6 @@ try param.writeParam(output_dir + "/simulation.param"); } -#if USE_NEW_PARSER std::cout << "\n\n================ Starting main simulation loop ===============\n" << std::flush; @@ -260,85 +229,6 @@ try fullReport.reportParam(tot_os); warnIfUnusedParams(param); } -#else - std::cout << "\n\n================ Starting main simulation loop ===============\n" - << " (number of epochs: " - << (deck->numberOfEpochs()) << ")\n\n" << std::flush; - - SimulatorReport rep; - // With a deck, we may have more epochs etc. - WellStateFullyImplicitBlackoil well_state; - int step = 0; - SimulatorTimer simtimer; - // Use timer for last epoch to obtain total time. - deck->setCurrentEpoch(deck->numberOfEpochs() - 1); - simtimer.init(*deck); - - const double total_time = simtimer.totalTime(); - for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { - // Set epoch index. - deck->setCurrentEpoch(epoch); - - // Update the timer. - if (deck->hasField("TSTEP")) { - simtimer.init(*deck); - } else { - if (epoch != 0) { - OPM_THROW(std::runtime_error, "No TSTEP in deck for epoch " << epoch); - } - simtimer.init(param); - } - simtimer.setCurrentStepNum(step); - simtimer.setTotalTime(total_time); - - // Report on start of epoch. - std::cout << "\n\n-------------- Starting epoch " << epoch << " --------------" - << "\n (number of steps: " - << simtimer.numSteps() - step << ")\n\n" << std::flush; - - // Create new wells, well_state - WellsManager wells(*deck, *grid->c_grid(), props->permeability()); - // @@@ HACK: we should really make a new well state and - // properly transfer old well state to it every epoch, - // since number of wells may change etc. - if (epoch == 0) { - well_state.init(wells.c_wells(), state); - } - - if (epoch == 0) - outputWriter.writeInit(simtimer, state, well_state.basicWellState()); - - // Create and run simulator. - SimulatorFullyImplicitBlackoil simulator(param, - *grid->c_grid(), - *new_props, - rock_comp->isActive() ? rock_comp.get() : 0, - wells, - linsolver, - grav); - outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState()); - - if (epoch == 0) { - warnIfUnusedParams(param); - } - SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); - if (output) { - epoch_rep.reportParam(outStream); - } - // Update total timing report and remember step number. - rep += epoch_rep; - step = simtimer.currentStepNum(); - } - - std::cout << "\n\n================ End of simulation ===============\n\n"; - rep.report(std::cout); - - if (output) { - std::string filename = output_dir + "/walltime.param"; - std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); - rep.reportParam(tot_os); - } -#endif } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; From 6c4114ea691fca2a410ee4f29aed30d4f90c9a5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 1 Apr 2014 15:49:01 +0200 Subject: [PATCH 25/28] Add more options for Selector class. You can now choose which comparison to do for the indicator vector to find when to choose the left argument. Only >= before (still default). --- opm/autodiff/AutoDiffHelpers.hpp | 34 ++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/AutoDiffHelpers.hpp b/opm/autodiff/AutoDiffHelpers.hpp index 9d26a6d36..8d196b2cf 100644 --- a/opm/autodiff/AutoDiffHelpers.hpp +++ b/opm/autodiff/AutoDiffHelpers.hpp @@ -292,7 +292,10 @@ spdiag(const AutoDiffBlock::V& d) public: typedef AutoDiffBlock ADB; - Selector(const typename ADB::V& selection_basis) + enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero }; + + Selector(const typename ADB::V& selection_basis, + CriterionForLeftElement crit = GreaterEqualZero) { // Define selector structure. const int n = selection_basis.size(); @@ -300,10 +303,33 @@ spdiag(const AutoDiffBlock::V& d) left_elems_.reserve(n); right_elems_.reserve(n); for (int i = 0; i < n; ++i) { - if (selection_basis[i] < 0.0) { - right_elems_.push_back(i); - } else { + bool chooseleft = false; + switch (crit) { + case GreaterEqualZero: + chooseleft = selection_basis[i] >= 0.0; + break; + case GreaterZero: + chooseleft = selection_basis[i] > 0.0; + break; + case Zero: + chooseleft = selection_basis[i] == 0.0; + break; + case NotEqualZero: + chooseleft = selection_basis[i] != 0.0; + break; + case LessZero: + chooseleft = selection_basis[i] < 0.0; + break; + case LessEqualZero: + chooseleft = selection_basis[i] <= 0.0; + break; + default: + OPM_THROW(std::logic_error, "No such criterion: " << crit); + } + if (chooseleft) { left_elems_.push_back(i); + } else { + right_elems_.push_back(i); } } } From 8437051416bbab56879104c4aaaf5ba623dee883 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 2 Apr 2014 17:06:20 +0200 Subject: [PATCH 26/28] make everything compile with the newly refactored EclipseWriter the states needed to be removed from the call to writeInit()... --- examples/sim_fibo_ad.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/sim_fibo_ad.cpp b/examples/sim_fibo_ad.cpp index dda83487b..348f983cb 100644 --- a/examples/sim_fibo_ad.cpp +++ b/examples/sim_fibo_ad.cpp @@ -200,7 +200,7 @@ try simtimer.setCurrentStepNum(reportStepIdx); if (reportStepIdx == 0) { - outputWriter.writeInit(simtimer, state, well_state.basicWellState()); + outputWriter.writeInit(simtimer); outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState()); } From a905eade7a32b73ed816875a9d2019f78e79c8b5 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 4 Apr 2014 11:33:36 +0200 Subject: [PATCH 27/28] Use preferred phase to compute wellbore mix for dead wells The wellmore mix is set to preferred phase for dead wells, where the total volumetric rates are zero. This sets the phase of the flow out of perforations in dead wells and thus avoids zero volumerats for injecting preforations in dead wells. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index e2a1cab13..873148629 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -792,9 +792,11 @@ namespace { isNotDeadWells[c] = 0; } // compute wellbore mixture at std conds + Selector notDeadWells_selector(wbqt.value(),Selector::Zero); std::vector mix_s(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { - mix_s[phase] = isNotDeadWells * wbq[phase]/wbqt; + const int pos = pu.phase_pos[phase]; + mix_s[phase] = notDeadWells_selector.select(ADB::constant(compi.col(pos), state.bhp.blockPattern()),wbq[phase]/wbqt); } @@ -832,7 +834,7 @@ namespace { 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); + volRat += tmp / subset(rq_[phase].b,well_cells); } From e699f9e1357fb5234c2e2c06be3f4e93dc6586d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sat, 5 Apr 2014 00:12:03 +0200 Subject: [PATCH 28/28] Remove unnecessary argument, and minor formatting corrections. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 873148629..d36d47590 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -787,16 +787,17 @@ namespace { // check for dead wells V isNotDeadWells = V::Constant(nw,1.0); - for (int c = 0; c < nw; ++c){ - if (wbqt.value()[c] == 0) - isNotDeadWells[c] = 0; + for (int w = 0; w < nw; ++w) { + if (wbqt.value()[w] == 0) { + isNotDeadWells[w] = 0; + } } // compute wellbore mixture at std conds - Selector notDeadWells_selector(wbqt.value(),Selector::Zero); + Selector notDeadWells_selector(wbqt.value(), Selector::Zero); std::vector mix_s(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { const int pos = pu.phase_pos[phase]; - mix_s[phase] = notDeadWells_selector.select(ADB::constant(compi.col(pos), state.bhp.blockPattern()),wbq[phase]/wbqt); + mix_s[phase] = notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); } @@ -834,7 +835,7 @@ namespace { 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); + volRat += tmp / subset(rq_[phase].b,well_cells); }