From 8acd441149ec5528b135cedd4c349649a98db0c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 14 Jun 2012 14:26:54 +0200 Subject: [PATCH] Added simulator program handling scheduling. --- examples/Makefile.am | 6 +- examples/sim_poly2p_incomp_reorder.cpp | 306 +++++++++++++++++++++++++ 2 files changed, 311 insertions(+), 1 deletion(-) create mode 100644 examples/sim_poly2p_incomp_reorder.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index 1c0864156..3cf56a2bb 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -10,8 +10,12 @@ $(BOOST_SYSTEM_LIB) noinst_PROGRAMS = \ -polymer_reorder +polymer_reorder \ +sim_poly2p_incomp_reorder polymer_reorder_SOURCES = \ polymer_reorder.cpp +sim_poly2p_incomp_reorder_SOURCES = \ +sim_poly2p_incomp_reorder.cpp + diff --git a/examples/sim_poly2p_incomp_reorder.cpp b/examples/sim_poly2p_incomp_reorder.cpp new file mode 100644 index 000000000..295d74e91 --- /dev/null +++ b/examples/sim_poly2p_incomp_reorder.cpp @@ -0,0 +1,306 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + std::cout << "\n================ Test program for incompressible two-phase flow with polymer ===============\n\n"; + 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; + PolymerState state; + Opm::PolymerProperties poly_props; + // 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 EclipseGridParser(deck_filename)); + // Grid init + 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 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 RockCompressibility(*deck)); + // 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); + } else { + initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); + } + // Init polymer properties. + poly_props.readFromDeck(*deck); + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + 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 GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + // Init Polymer state + if (param.has("poly_init")) { + double poly_init = param.getDefault("poly_init", 0.0); + for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) { + double smin[2], smax[2]; + props->satRange(1, &cell, smin, smax); + if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) { + state.concentration()[cell] = poly_init; + state.maxconcentration()[cell] = poly_init; + } else { + state.saturation()[2*cell + 0] = 0.; + state.saturation()[2*cell + 1] = 1.; + state.concentration()[cell] = 0.; + state.maxconcentration()[cell] = 0.; + } + } + } + // Init polymer properties. + // Setting defaults to provide a simple example case. + double c_max = param.getDefault("c_max_limit", 5.0); + double mix_param = param.getDefault("mix_param", 1.0); + double rock_density = param.getDefault("rock_density", 1000.0); + double dead_pore_vol = param.getDefault("dead_pore_vol", 0.15); + double res_factor = param.getDefault("res_factor", 1.) ; // res_factor = 1 gives no change in permeability + double c_max_ads = param.getDefault("c_max_ads", 1.); + int ads_index = param.getDefault("ads_index", Opm::PolymerProperties::NoDesorption); + std::vector c_vals_visc(2, -1e100); + c_vals_visc[0] = 0.0; + c_vals_visc[1] = 7.0; + std::vector visc_mult_vals(2, -1e100); + visc_mult_vals[0] = 1.0; + // poly_props.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0); + visc_mult_vals[1] = 20.0; + std::vector c_vals_ads(3, -1e100); + c_vals_ads[0] = 0.0; + c_vals_ads[1] = 2.0; + c_vals_ads[2] = 8.0; + std::vector ads_vals(3, -1e100); + ads_vals[0] = 0.0; + ads_vals[1] = 0.0015; + ads_vals[2] = 0.0025; + // ads_vals[1] = 0.0; + // ads_vals[2] = 0.0; + poly_props.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads, + static_cast(ads_index), + c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals); + } + + // Warn if gravity but no density difference. + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + if (use_gravity) { + if (props->density()[0] == props->density()[1]) { + std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; + } + } + 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/unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + 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(), FlowBCManager::Side(pside), pside_pressure); + } + + // Linear solver. + LinearSolverFactory linsolver(param); + + // Warn if any parameters are unused. + // if (param.anyUnused()) { + // std::cout << "-------------------- Unused parameters: --------------------\n"; + // param.displayUsage(); + // std::cout << "----------------------------------------------------------------" << std::endl; + // } + + // Write parameters used for later reference. + // if (output) { + // param.writeParam(output_dir + "/spu_2p.param"); + // } + + + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << " (number of epochs: " + << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush; + + SimulatorReport rep; + if (!use_deck) { + // Simple simulation without a deck. + SimulatorPolymer simulator(param, + *grid->c_grid(), + *props, + poly_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); + 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); + + // Update the timer. + if (deck->hasField("TSTEP")) { + simtimer.init(*deck); + } else { + if (epoch != 0) { + THROW("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_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. + SimulatorPolymer simulator(param, + *grid->c_grid(), + *props, + poly_props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells.c_wells(), + src, + bcs.c_bcs(), + linsolver, + grav); + 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); +}