From 5e7064b337386d38b2cd68581631889e47e931ee Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 31 Mar 2014 16:39:17 +0200 Subject: [PATCH 1/6] 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 2/6] 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 3/6] 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 4/6] 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 5/6] 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); } From bccc7c1350810683122b12d791978d8ff42fa1e8 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 9 Apr 2014 12:16:33 +0200 Subject: [PATCH 6/6] Updates interface to fix issues not resolved in the latest merge. Namely the template parameter was missing for the added methods. In addition there were still various instances where the grid was used directly rather than via the functions in GridHelpers.hh. --- .../FullyImplicitBlackoilSolver_impl.hpp | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index 45ce2839a..59d801387 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -562,6 +562,7 @@ namespace { void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw) { + using namespace Opm::AutoDiffGrid; // 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. @@ -599,7 +600,7 @@ namespace { // 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 depth = cellCentroidsZ(grid_); const V pdepth = subset(depth, well_cells); std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); // Surface density. @@ -607,7 +608,7 @@ namespace { // Gravity double grav = 0.0; const double* g = geo_.gravity(); - const int dim = grid_.dimensions; + const int dim = dimensions(grid_); if (g) { // Guard against gravity in anything but last dimension. for (int dd = 0; dd < dim - 1; ++dd) { @@ -702,10 +703,12 @@ namespace { - void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state, + template + void FullyImplicitBlackoilSolver::addWellEq(const SolutionState& state, WellStateFullyImplicitBlackoil& xw) { - const int nc = grid_.number_of_cells; + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); const int np = wells_.number_of_phases; const int nw = wells_.number_of_wells; const int nperf = wells_.well_connpos[nw]; @@ -953,8 +956,8 @@ namespace { - - void FullyImplicitBlackoilSolver::updateWellControls(const ADB& bhp, + template + void FullyImplicitBlackoilSolver::updateWellControls(const ADB& bhp, const ADB& well_phase_flow_rate, WellStateFullyImplicitBlackoil& xw) const { @@ -1007,8 +1010,8 @@ namespace { - - void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state, + template + void FullyImplicitBlackoilSolver::addWellControlEq(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw) { // Handling BHP and SURFACE_RATE wells. @@ -1052,8 +1055,8 @@ namespace { - - void FullyImplicitBlackoilSolver::addOldWellEq(const SolutionState& state) + template + void FullyImplicitBlackoilSolver::addOldWellEq(const SolutionState& state) { // -------- Well equation, and well contributions to the mass balance equations --------