diff --git a/examples/flow_polymer.cpp b/examples/flow_polymer.cpp new file mode 100644 index 000000000..7973fe934 --- /dev/null +++ b/examples/flow_polymer.cpp @@ -0,0 +1,284 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2014 STATOIL ASA. + + 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 . +*/ +#include "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 +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + + +namespace +{ + void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) + { + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + } +} // anon namespace + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +try +{ + using namespace Opm; + + std::cout << "\n================ Test program for fully implicit three-phase black-oil-polymer flow ===============\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"); + if (!use_deck) { + OPM_THROW(std::runtime_error, "This program must be run with an input deck. " + "Specify the deck with deck_filename=deckname.data (for example)."); + } + std::shared_ptr grid; + std::shared_ptr props; + std::shared_ptr new_props; + std::shared_ptr rock_comp; + PolymerBlackoilState state; + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + std::string deck_filename = param.get("deck_filename"); + + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::string output_dir; + if (output) { + // Create output directory if needed. + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + // Write simulation parameters. + param.writeParam(output_dir + "/simulation.param"); + } + + std::string logFile = output_dir + "/LOGFILE.txt"; + Opm::ParserPtr parser(new Opm::Parser()); + { + std::shared_ptr streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); + std::shared_ptr counterLog = std::make_shared(Opm::Log::DefaultMessageTypes); + + Opm::OpmLog::addBackend( "STREAM" , streamLog ); + Opm::OpmLog::addBackend( "COUNTER" , counterLog ); + } + + + Opm::DeckConstPtr deck; + std::shared_ptr eclipseState; + try { + deck = parser->parseFile(deck_filename); + Opm::checkDeck(deck); + eclipseState.reset(new Opm::EclipseState(deck)); + } + catch (const std::invalid_argument& e) { + std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl; + std::cerr << "Exception caught: " << e.what() << std::endl; + return EXIT_FAILURE; + } + + // Grid init + std::vector porv = eclipseState->getDoubleGridProperty("PORV")->getData(); + grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv)); + auto &cGrid = *grid->c_grid(); + const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); + Opm::BlackoilOutputWriter outputWriter(cGrid, + param, + eclipseState, + pu ); + + // Rock and fluid init + props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param)); + new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid())); + const bool polymer = deck->hasKeyword("POLYMER"); + const bool use_wpolymer = deck->hasKeyword("WPOLYMER"); + PolymerProperties polymer_props(deck, eclipseState); + PolymerPropsAd polymer_props_ad(polymer_props); + // 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, eclipseState)); + + // Gravity. + gravity[2] = deck->hasKeyword("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.blackoilState()); + initBlackoilSurfvol(*grid->c_grid(), *props, state.blackoilState()); + enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; + if (pu.phase_used[Oil] && pu.phase_used[Gas]) { + const int np = props->numPhases(); + const int nc = grid->c_grid()->number_of_cells; + for (int c = 0; c < nc; ++c) { + state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]] + / state.surfacevol()[c*np + pu.phase_pos[Oil]]; + } + } + } else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) { + state.init(*grid->c_grid(), props->numPhases()); + const double grav = param.getDefault("gravity", unit::gravity); + initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state.blackoilState()); + state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0); + } else { + initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState()); + } + + // The capillary pressure is scaled in new_props to match the scaled capillary pressure in props. + if (deck->hasKeyword("SWATINIT")) { + const int nc = grid->c_grid()->number_of_cells; + std::vector cells(nc); + for (int c = 0; c < nc; ++c) { cells[c] = c; } + std::vector pc = state.saturation(); + props->capPress(nc, state.saturation().data(), cells.data(), pc.data(),NULL); + new_props->setSwatInitScaling(state.saturation(),pc); + } + + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + const double *grav = use_gravity ? &gravity[0] : 0; + + // Solver for Newton iterations. + std::unique_ptr fis_solver; + if (param.getDefault("use_cpr", true)) { + fis_solver.reset(new NewtonIterationBlackoilCPR(param)); + } else { + fis_solver.reset(new NewtonIterationBlackoilSimple(param)); + } + + Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap()); + SimulatorTimer simtimer; + + // initialize variables + simtimer.init(timeMap); + if (polymer){ + if (!use_wpolymer) { + OPM_MESSAGE("Warning: simulate polymer injection without WPOLYMER."); + } else { + if (param.has("polymer_start_days")) { + OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck." + "You seem to be trying to control it via parameter poly_start_days (etc.) as well."); + } + } + } else { + if (use_wpolymer) { + OPM_MESSAGE("Warning: use WPOLYMER in a non-polymer scenario."); + } + } + + bool use_local_perm = param.getDefault("use_local_perm", true); + Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, use_local_perm, grav); + + std::vector threshold_pressures = thresholdPressures(eclipseState, *grid->c_grid()); + + SimulatorFullyImplicitBlackoilPolymer simulator(param, + *grid->c_grid(), + geology, + *new_props, + polymer_props_ad, + rock_comp->isActive() ? rock_comp.get() : 0, + *fis_solver, + grav, + deck->hasKeyword("DISGAS"), + deck->hasKeyword("VAPOIL"), + polymer, + eclipseState, + outputWriter, + deck, + threshold_pressures); + + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << std::flush; + + SimulatorReport fullReport = simulator.run(simtimer, state); + + std::cout << "\n\n================ End of simulation ===============\n\n"; + fullReport.report(std::cout); + + if (output) { + std::string filename = output_dir + "/walltime.txt"; + std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); + fullReport.reportParam(tot_os); + warnIfUnusedParams(param); + } +} +catch (const std::exception &e) { + std::cerr << "Program threw an exception: " << e.what() << "\n"; + throw; +} + diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp index af74ffe21..f285c4650 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp @@ -1310,15 +1310,19 @@ namespace detail { // Update primary variables, if necessary. if (bhp_changed) { ADB::V new_bhp = Eigen::Map(xw.bhp().data(), nw); - bhp = ADB::function(new_bhp, bhp.derivative()); + // Avoiding the copy below would require a value setter method + // in AutoDiffBlock. + std::vector old_derivs = bhp.derivative(); + bhp = ADB::function(std::move(new_bhp), std::move(old_derivs)); } if (rates_changed) { // Need to reshuffle well rates, from phase running fastest // to wells running fastest. // The transpose() below switches the ordering. const DataBlock wrates = Eigen::Map(xw.wellRates().data(), nw, np).transpose(); - const ADB::V new_qs = Eigen::Map(wrates.data(), nw*np); - well_phase_flow_rate = ADB::function(new_qs, well_phase_flow_rate.derivative()); + ADB::V new_qs = Eigen::Map(wrates.data(), nw*np); + std::vector old_derivs = well_phase_flow_rate.derivative(); + well_phase_flow_rate = ADB::function(std::move(new_qs), std::move(old_derivs)); } } @@ -2396,8 +2400,8 @@ namespace detail { for (int block = 0; block < num_blocks; ++block) { fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]); } - return ADB::function(pm, jacs); - } else { + return ADB::function(std::move(pm), std::move(jacs)); + } else { return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); } } @@ -2424,7 +2428,7 @@ namespace detail { for (int block = 0; block < num_blocks; ++block) { fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]); } - return ADB::function(tm, jacs); + return ADB::function(std::move(tm), std::move(jacs)); } else { return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); } diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 485f9ac46..7ceec001e 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -994,7 +994,7 @@ namespace { for (int block = 0; block < num_blocks; ++block) { jacs[block] = dpm_diag * p.derivative()[block]; } - return ADB::function(pm, jacs); + return ADB::function(std::move(pm), std::move(jacs)); } else { return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); } @@ -1021,7 +1021,7 @@ namespace { for (int block = 0; block < num_blocks; ++block) { jacs[block] = dtm_diag * p.derivative()[block]; } - return ADB::function(tm, jacs); + return ADB::function(std::move(tm), std::move(jacs)); } else { return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); } diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp index ed0123a85..184426a31 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp @@ -117,7 +117,7 @@ namespace Opm { for (int block = 0; block < num_blocks; ++block) { jacs[block] = dim_diag * c.derivative()[block]; } - return ADB::function(inv_mu_w_eff, jacs); + return ADB::function(std::move(inv_mu_w_eff), std::move(jacs)); } @@ -165,7 +165,7 @@ namespace Opm { jacs[block] = dmc_diag * c.derivative()[block]; } - return ADB::function(mc, jacs); + return ADB::function(std::move(mc), std::move(jacs)); } @@ -213,7 +213,7 @@ namespace Opm { jacs[block] = dads_diag * c.derivative()[block]; } - return ADB::function(ads, jacs); + return ADB::function(std::move(ads), std::move(jacs)); }