diff --git a/examples/sim_2p_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp index 822b8e18..35c49146 100644 --- a/examples/sim_2p_incomp_reorder.cpp +++ b/examples/sim_2p_incomp_reorder.cpp @@ -22,7 +22,6 @@ #include "config.h" #endif // HAVE_CONFIG_H -#include #include #include @@ -32,9 +31,6 @@ #include #include #include -#include -#include -#include #include #include @@ -49,18 +45,9 @@ #include #include -#include - -#include -#include #include -#include -#include #include -#include -#include -#include #include #include @@ -74,34 +61,34 @@ main(int argc, char** argv) using namespace Opm; std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; - Opm::parameter::ParameterGroup param(argc, argv, false); + parameter::ParameterGroup param(argc, argv, false); std::cout << "--------------- Reading parameters ---------------" << std::endl; // If we have a "deck_filename", grid and props will be read from that. bool use_deck = param.has("deck_filename"); - boost::scoped_ptr deck; - boost::scoped_ptr grid; - boost::scoped_ptr props; - boost::scoped_ptr rock_comp; - Opm::TwophaseState state; + boost::scoped_ptr deck; + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr rock_comp; + 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"); - deck.reset(new Opm::EclipseGridParser(deck_filename)); + deck.reset(new EclipseGridParser(deck_filename)); // Grid init - grid.reset(new Opm::GridManager(*deck)); + grid.reset(new GridManager(*deck)); // Rock and fluid init const int* gc = grid->c_grid()->global_cell; std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); - props.reset(new Opm::IncompPropertiesFromDeck(*deck, global_cell)); + props.reset(new IncompPropertiesFromDeck(*deck, global_cell)); // check_well_controls = param.getDefault("check_well_controls", false); // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Rock compressibility. - rock_comp.reset(new Opm::RockCompressibility(*deck)); + rock_comp.reset(new RockCompressibility(*deck)); // Gravity. - gravity[2] = deck->hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; + gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; // Init state variables (saturation and pressure). if (param.has("init_saturation")) { initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); @@ -116,11 +103,11 @@ main(int argc, char** argv) const double dx = param.getDefault("dx", 1.0); const double dy = param.getDefault("dy", 1.0); const double dz = param.getDefault("dz", 1.0); - grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); + grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); // Rock and fluid init. - props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); // Rock compressibility. - rock_comp.reset(new Opm::RockCompressibility(param)); + rock_comp.reset(new RockCompressibility(param)); // Gravity. gravity[2] = param.getDefault("gravity", 0.0); // Init state variables (saturation and pressure). @@ -134,47 +121,40 @@ main(int argc, char** argv) std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; } } - - // Source-related variables init. - int num_cells = grid->c_grid()->number_of_cells; - std::vector totmob; - std::vector omega; // Will remain empty if no gravity. - std::vector rc; // Will remain empty if no rock compressibility. - - // Extra rock init. - std::vector porevol; - if (rock_comp->isActive()) { - computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); - } else { - computePorevolume(*grid->c_grid(), props->porosity(), porevol); - } - double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + const double *grav = use_gravity ? &gravity[0] : 0; // Initialising src + int num_cells = grid->c_grid()->number_of_cells; std::vector src(num_cells, 0.0); if (use_deck) { // Do nothing, wells will be the driving force, not source terms. } else { + // Compute pore volumes, in order to enable specifying injection rate + // terms of total pore volume. + std::vector porevol; + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } else { + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); const double default_injection = use_gravity ? 0.0 : 0.1; const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) - *tot_porevol_init/Opm::unit::day; + *tot_porevol_init/unit::day; src[0] = flow_per_sec; src[num_cells - 1] = -flow_per_sec; } // Boundary conditions. - Opm::FlowBCManager bcs; + FlowBCManager bcs; if (param.getDefault("use_pside", false)) { int pside = param.get("pside"); double pside_pressure = param.get("pside_pressure"); - bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure); + bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure); } // Linear solver. - Opm::LinearSolverFactory linsolver(param); - - const double *grav = use_gravity ? &gravity[0] : 0; - + LinearSolverFactory linsolver(param); // Warn if any parameters are unused. // if (param.anyUnused()) { @@ -189,45 +169,41 @@ main(int argc, char** argv) // } + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << " (number of epochs: " + << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + + SimulatorTwophase::SimulatorReport rep; if (!use_deck) { // Simple simulation without a deck. - Opm::SimulatorTwophase simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - 0, // wells - src, - bcs.c_bcs(), - linsolver, - grav); - Opm::SimulatorTimer simtimer; + SimulatorTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + 0, // wells + src, + bcs.c_bcs(), + linsolver, + grav); + SimulatorTimer simtimer; simtimer.init(param); WellState well_state; well_state.init(0, state); - simulator.run(simtimer, state, well_state); + rep = simulator.run(simtimer, state, well_state); } else { // With a deck, we may have more epochs etc. WellState 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); - Opm::WellsManager wells(*deck, *grid->c_grid(), props->permeability()); - Opm::SimulatorTwophase simulator(param, - *grid->c_grid(), - *props, - rock_comp->isActive() ? rock_comp.get() : 0, - wells.c_wells(), - src, - bcs.c_bcs(), - linsolver, - grav); - // @@@ 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); - } - Opm::SimulatorTimer simtimer; + + // Update the timer. if (deck->hasField("TSTEP")) { simtimer.init(*deck); } else { @@ -237,8 +213,41 @@ main(int argc, char** argv) simtimer.init(param); } simtimer.setCurrentStepNum(step); - simulator.run(simtimer, state, well_state); + 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_satate + 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); + } + + // Create and run simulator. + SimulatorTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells.c_wells(), + src, + bcs.c_bcs(), + linsolver, + grav); + SimulatorTwophase::SimulatorReport epoch_rep + = simulator.run(simtimer, state, well_state); + + // 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); }