diff --git a/.gitignore b/.gitignore index 9c1b8581d..18532844f 100644 --- a/.gitignore +++ b/.gitignore @@ -68,4 +68,10 @@ examples/upscale_* examples/grdecldips data steadystate_test_implicit -steadystate_test_explicit \ No newline at end of file +steadystate_test_explicit +======= +*.iml +.idea +*/opm-simulation + + diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp new file mode 100644 index 000000000..53220be7f --- /dev/null +++ b/examples/sim_poly_fi2p_comp_ad.cpp @@ -0,0 +1,225 @@ +/* + Copyright 2014 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 . +*/ + + +#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 +#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 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"); + + Opm::ParserPtr newParser(new Opm::Parser()); + bool strict_parsing = param.getDefault("strict_parsing", true); + Opm::DeckConstPtr deck = newParser->parseFile(deck_filename, strict_parsing); + std::shared_ptr eclipseState(new EclipseState(deck)); + // Grid init + std::vector porv; + if (eclipseState->hasDoubleGridProperty("PORV")) { + porv = eclipseState->getDoubleGridProperty("PORV")->getData(); + } + grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv)); + auto &cGrid = *grid->c_grid(); + const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); + Opm::EclipseWriter outputWriter(param, + eclipseState, + pu, + cGrid.number_of_cells, + cGrid.global_cell); + // Rock and fluid init + props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param)); + new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid())); + PolymerProperties polymer_props(deck, eclipseState); + PolymerPropsAd polymer_props_ad(polymer_props); + // 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); + initBlackoilSurfvol(*grid->c_grid(), *props, state); + } else { + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + } + + 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)); + } + + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::string output_dir; + if (output) { + 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); + } + param.writeParam(output_dir + "/simulation.param"); + } + + Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap()); + SimulatorTimer simtimer; + simtimer.init(timeMap); + + + SimulatorReport rep; + // With a deck, we may have more epochs etc. + WellState well_state; + // Check for WPOLYMER presence in last epoch to decide + // polymer injection control type. + const bool use_wpolymer = deck->hasKeyword("WPOLYMER"); + if (use_wpolymer) { + if (param.has("poly_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."); + } + } + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << std::flush; + SimulatorReport fullReport; + // Create and run simulator. + Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav); + SimulatorFullyImplicitCompressiblePolymer simulator(param, + *grid->c_grid(), + geology, + *new_props, + polymer_props_ad, + rock_comp->isActive() ? rock_comp.get() : 0, + eclipseState, + outputWriter, + deck, + *fis_solver, + grav); + 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.param"; + 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/examples/sim_poly_fibo_ad.cpp b/examples/sim_poly_fibo_ad.cpp new file mode 100644 index 000000000..3f1536a16 --- /dev/null +++ b/examples/sim_poly_fibo_ad.cpp @@ -0,0 +1,253 @@ +/* + Copyright 2014 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 + + +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 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"); + + Opm::ParserPtr newParser(new Opm::Parser() ); + bool strict_parsing = param.getDefault("strict_parsing", true); + Opm::DeckConstPtr deck = newParser->parseFile(deck_filename, strict_parsing); + std::shared_ptr eclipseState(new EclipseState(deck)); + + // Grid init + std::vector porv; + if (eclipseState->hasDoubleGridProperty("PORV")) { + porv = eclipseState->getDoubleGridProperty("PORV")->getData(); + } + grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv)); + auto &cGrid = *grid->c_grid(); + const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); + Opm::EclipseWriter outputWriter(param, + eclipseState, + pu, + cGrid.number_of_cells, + cGrid.global_cell); + + // 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); + initBlackoilSurfvol(*grid->c_grid(), *props, state); + 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 { + state.init(*grid->c_grid(), props->numPhases()); + initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState()); + } + + 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)); + } + + // 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"); + } + + 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."); + } + } + std::cout << "\n\n================ Starting main simulation loop ===============\n" + << std::flush; + SimulatorReport fullReport; + Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav); + + std::vector threshold_pressures = thresholdPressures(deck, 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); + + + 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/PolymerBlackoilState.hpp b/opm/polymer/PolymerBlackoilState.hpp index 8116ec2c8..078ac27b2 100644 --- a/opm/polymer/PolymerBlackoilState.hpp +++ b/opm/polymer/PolymerBlackoilState.hpp @@ -44,6 +44,10 @@ namespace Opm concentration_.resize(number_of_cells, 0.0); cmax_.resize(number_of_cells, 0.0); } + int numPhases() const + { + return state_blackoil_.numPhases(); + } enum ExtremalSat { MinSat = BlackoilState::MinSat, MaxSat = BlackoilState::MaxSat }; @@ -61,6 +65,8 @@ namespace Opm std::vector& facepressure() { return state_blackoil_.facepressure(); } std::vector& faceflux () { return state_blackoil_.faceflux(); } std::vector& saturation () { return state_blackoil_.saturation(); } + std::vector& gasoilratio () { return state_blackoil_.gasoilratio(); } + std::vector& rv () { return state_blackoil_.rv(); } std::vector& concentration() { return concentration_; } std::vector& maxconcentration() { return cmax_; } @@ -69,6 +75,8 @@ namespace Opm const std::vector& facepressure() const { return state_blackoil_.facepressure(); } const std::vector& faceflux () const { return state_blackoil_.faceflux(); } const std::vector& saturation () const { return state_blackoil_.saturation(); } + const std::vector& gasoilratio() const { return state_blackoil_.gasoilratio(); } + const std::vector& rv () const { return state_blackoil_.rv(); } const std::vector& concentration() const { return concentration_; } const std::vector& maxconcentration() const { return cmax_; } diff --git a/opm/polymer/PolymerState.hpp b/opm/polymer/PolymerState.hpp index d5928edf9..d2374fb26 100644 --- a/opm/polymer/PolymerState.hpp +++ b/opm/polymer/PolymerState.hpp @@ -55,6 +55,10 @@ namespace Opm state2p_.setFirstSat(cells, props, static_cast(es)); } + inline int numPhases() const + { + return state2p_.numPhases(); + } std::vector& pressure () { return state2p_.pressure(); } std::vector& facepressure() { return state2p_.facepressure(); } std::vector& faceflux () { return state2p_.faceflux(); } @@ -71,7 +75,7 @@ namespace Opm TwophaseState& twophaseState() { return state2p_; } const TwophaseState& twophaseState() const { return state2p_; } - + private: TwophaseState state2p_; std::vector concentration_; diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp index 544543d78..dbb2efe81 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.cpp +++ b/opm/polymer/SimulatorCompressiblePolymer.cpp @@ -292,207 +292,206 @@ namespace Opm wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(), state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); } - for (; !timer.done(); ++timer) { - // Report timestep and (optionally) write state to disk. - timer.report(std::cout); - if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + // Report timestep and (optionally) write state to disk. + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + } + + initial_pressure = state.pressure(); + + // Solve pressure equation. + if (check_well_controls_) { + computeFractionalFlow(props_, poly_props_, allcells_, + state.pressure(), state.surfacevol(), state.saturation(), + state.concentration(), state.maxconcentration(), + fractional_flows); + wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); + } + bool well_control_passed = !check_well_controls_; + int well_control_iteration = 0; + do { + // Run solver + pressure_timer.start(); + psolver_.solve(timer.currentStepLength(), state, well_state); + + // Renormalize pressure if both fluids and rock are + // incompressible, and there are no pressure + // conditions (bcs or wells). It is deemed sufficient + // for now to renormalize using geometric volume + // instead of pore volume. + if (psolver_.singularPressure()) { + // Compute average pressures of previous and last + // step, and total volume. + double av_prev_press = 0.0; + double av_press = 0.0; + double tot_vol = 0.0; + const int num_cells = grid_.number_of_cells; + for (int cell = 0; cell < num_cells; ++cell) { + av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell]; + av_press += state.pressure()[cell]*grid_.cell_volumes[cell]; + tot_vol += grid_.cell_volumes[cell]; + } + // Renormalization constant + const double ren_const = (av_prev_press - av_press)/tot_vol; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += ren_const; + } + const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells; + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += ren_const; } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); } - initial_pressure = state.pressure(); + // Stop timer and report + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; - // Solve pressure equation. + // Optionally, check if well controls are satisfied. if (check_well_controls_) { - computeFractionalFlow(props_, poly_props_, allcells_, - state.pressure(), state.surfacevol(), state.saturation(), - state.concentration(), state.maxconcentration(), - fractional_flows); - wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); - } - bool well_control_passed = !check_well_controls_; - int well_control_iteration = 0; - do { - // Run solver - pressure_timer.start(); - psolver_.solve(timer.currentStepLength(), state, well_state); - - // Renormalize pressure if both fluids and rock are - // incompressible, and there are no pressure - // conditions (bcs or wells). It is deemed sufficient - // for now to renormalize using geometric volume - // instead of pore volume. - if (psolver_.singularPressure()) { - // Compute average pressures of previous and last - // step, and total volume. - double av_prev_press = 0.0; - double av_press = 0.0; - double tot_vol = 0.0; - const int num_cells = grid_.number_of_cells; - for (int cell = 0; cell < num_cells; ++cell) { - av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell]; - av_press += state.pressure()[cell]*grid_.cell_volumes[cell]; - tot_vol += grid_.cell_volumes[cell]; - } - // Renormalization constant - const double ren_const = (av_prev_press - av_press)/tot_vol; - for (int cell = 0; cell < num_cells; ++cell) { - state.pressure()[cell] += ren_const; - } - const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells; - for (int well = 0; well < num_wells; ++well) { - well_state.bhp()[well] += ren_const; - } + Opm::computePhaseFlowRatesPerWell(*wells_, + well_state.perfRates(), + fractional_flows, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { + OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); } - - // Stop timer and report - pressure_timer.stop(); - double pt = pressure_timer.secsSinceStart(); - std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; - ptime += pt; - - // Optionally, check if well controls are satisfied. - if (check_well_controls_) { - Opm::computePhaseFlowRatesPerWell(*wells_, - well_state.perfRates(), - fractional_flows, - well_resflows_phase); - std::cout << "Checking well conditions." << std::endl; - // For testing we set surface := reservoir - well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); - ++well_control_iteration; - if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { - OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); - } - if (!well_control_passed) { - std::cout << "Well controls not passed, solving again." << std::endl; - } else { - std::cout << "Well conditions met." << std::endl; - } - } - } while (!well_control_passed); - - // Update pore volumes if rock is compressible. - if (rock_comp_props_ && rock_comp_props_->isActive()) { - initial_porevol = porevol; - computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); - } - - // Process transport sources (to include bdy terms and well flows). - Opm::computeTransportSource(props_, wells_, well_state, transport_src); - - // Find inflow rate. - const double current_time = timer.simulationTimeElapsed(); - double stepsize = timer.currentStepLength(); - polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c); - - // Solve transport. - transport_timer.start(); - if (num_transport_substeps_ != 1) { - stepsize /= double(num_transport_substeps_); - std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; - } - double injected[2] = { 0.0 }; - double produced[2] = { 0.0 }; - double polyinj = 0.0; - double polyprod = 0.0; - for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { - tsolver_.solve(&state.faceflux()[0], initial_pressure, - state.pressure(), &initial_porevol[0], &porevol[0], - &transport_src[0], &polymer_inflow_c[0], stepsize, - state.saturation(), state.surfacevol(), - state.concentration(), state.maxconcentration()); - double substep_injected[2] = { 0.0 }; - double substep_produced[2] = { 0.0 }; - double substep_polyinj = 0.0; - double substep_polyprod = 0.0; - Opm::computeInjectedProduced(props_, poly_props_, - state, - transport_src, polymer_inflow_c, stepsize, - substep_injected, substep_produced, - substep_polyinj, substep_polyprod); - injected[0] += substep_injected[0]; - injected[1] += substep_injected[1]; - produced[0] += substep_produced[0]; - produced[1] += substep_produced[1]; - polyinj += substep_polyinj; - polyprod += substep_polyprod; - if (gravity_ != 0 && use_segregation_split_) { - tsolver_.solveGravity(columns_, stepsize, - state.saturation(), state.surfacevol(), - state.concentration(), state.maxconcentration()); + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; } } - transport_timer.stop(); - double tt = transport_timer.secsSinceStart(); - std::cout << "Transport solver took: " << tt << " seconds." << std::endl; - ttime += tt; + } while (!well_control_passed); - // Report volume balances. - Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol); - polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); - polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_, - state, rock_comp_props_); - tot_injected[0] += injected[0]; - tot_injected[1] += injected[1]; - tot_produced[0] += produced[0]; - tot_produced[1] += produced[1]; - tot_polyinj += polyinj; - tot_polyprod += polyprod; - std::cout.precision(5); - const int width = 18; - std::cout << "\nMass balance: " - " water(surfvol) oil(surfvol) polymer(kg)\n"; - std::cout << " In-place: " - << std::setw(width) << inplace_surfvol[0] - << std::setw(width) << inplace_surfvol[1] - << std::setw(width) << polymass << std::endl; - std::cout << " Adsorbed: " - << std::setw(width) << 0.0 - << std::setw(width) << 0.0 - << std::setw(width) << polymass_adsorbed << std::endl; - std::cout << " Injected: " - << std::setw(width) << injected[0] - << std::setw(width) << injected[1] - << std::setw(width) << polyinj << std::endl; - std::cout << " Produced: " - << std::setw(width) << produced[0] - << std::setw(width) << produced[1] - << std::setw(width) << polyprod << std::endl; - std::cout << " Total inj: " - << std::setw(width) << tot_injected[0] - << std::setw(width) << tot_injected[1] - << std::setw(width) << tot_polyinj << std::endl; - std::cout << " Total prod: " - << std::setw(width) << tot_produced[0] - << std::setw(width) << tot_produced[1] - << std::setw(width) << tot_polyprod << std::endl; - const double balance[3] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0], - init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1], - init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed }; - std::cout << " Initial - inplace + inj - prod: " - << std::setw(width) << balance[0] - << std::setw(width) << balance[1] - << std::setw(width) << balance[2] - << std::endl; - std::cout << " Relative mass error: " - << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0]) - << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1]) - << std::setw(width) << balance[2]/(init_polymass + tot_polyinj) - << std::endl; - std::cout.precision(8); + // Update pore volumes if rock is compressible. + if (rock_comp_props_ && rock_comp_props_->isActive()) { + initial_porevol = porevol; + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } - watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(), - produced[0]/(produced[0] + produced[1]), - tot_produced[0]/tot_porevol_init); - if (wells_) { - wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(), - state.saturation(), timer.simulationTimeElapsed() + timer.currentStepLength(), - well_state.bhp(), well_state.perfRates()); + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(props_, wells_, well_state, transport_src); + + // Find inflow rate. + const double current_time = timer.simulationTimeElapsed(); + double stepsize = timer.currentStepLength(); + polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c); + + + // Solve transport. + transport_timer.start(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double polyinj = 0.0; + double polyprod = 0.0; + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], initial_pressure, + state.pressure(), &initial_porevol[0], &porevol[0], + &transport_src[0], &polymer_inflow_c[0], stepsize, + state.saturation(), state.surfacevol(), + state.concentration(), state.maxconcentration()); + double substep_injected[2] = { 0.0 }; + double substep_produced[2] = { 0.0 }; + double substep_polyinj = 0.0; + double substep_polyprod = 0.0; + Opm::computeInjectedProduced(props_, poly_props_, + state, + transport_src, polymer_inflow_c, stepsize, + substep_injected, substep_produced, + substep_polyinj, substep_polyprod); + injected[0] += substep_injected[0]; + injected[1] += substep_injected[1]; + produced[0] += substep_produced[0]; + produced[1] += substep_produced[1]; + polyinj += substep_polyinj; + polyprod += substep_polyprod; + if (gravity_ != 0 && use_segregation_split_) { + tsolver_.solveGravity(columns_, stepsize, + state.saturation(), state.surfacevol(), + state.concentration(), state.maxconcentration()); } } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol); + polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); + polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_, + state, rock_comp_props_); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + tot_polyinj += polyinj; + tot_polyprod += polyprod; + std::cout.precision(5); + const int width = 18; + std::cout << "\nMass balance: " + " water(surfvol) oil(surfvol) polymer(kg)\n"; + std::cout << " In-place: " + << std::setw(width) << inplace_surfvol[0] + << std::setw(width) << inplace_surfvol[1] + << std::setw(width) << polymass << std::endl; + std::cout << " Adsorbed: " + << std::setw(width) << 0.0 + << std::setw(width) << 0.0 + << std::setw(width) << polymass_adsorbed << std::endl; + std::cout << " Injected: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] + << std::setw(width) << polyinj << std::endl; + std::cout << " Produced: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] + << std::setw(width) << polyprod << std::endl; + std::cout << " Total inj: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] + << std::setw(width) << tot_polyinj << std::endl; + std::cout << " Total prod: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] + << std::setw(width) << tot_polyprod << std::endl; + const double balance[3] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0], + init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1], + init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed }; + std::cout << " Initial - inplace + inj - prod: " + << std::setw(width) << balance[0] + << std::setw(width) << balance[1] + << std::setw(width) << balance[2] + << std::endl; + std::cout << " Relative mass error: " + << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0]) + << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1]) + << std::setw(width) << balance[2]/(init_polymass + tot_polyinj) + << std::endl; + std::cout.precision(8); + + watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(), + state.saturation(), timer.simulationTimeElapsed() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } if (output_) { if (output_vtk_) { diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp index 827b80ddd..8d943d328 100644 --- a/opm/polymer/SimulatorPolymer.cpp +++ b/opm/polymer/SimulatorPolymer.cpp @@ -315,196 +315,194 @@ namespace Opm well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); } - for (; !timer.done(); ++timer) { - // Report timestep and (optionally) write state to disk. - timer.report(std::cout); - if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + // Report timestep and (optionally) write state to disk. + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + if (output_binary_) { + outputStateBinary(grid_, state, timer, output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + } + + // Solve pressure. + if (check_well_controls_) { + computeFractionalFlow(props_, poly_props_, allcells_, + state.saturation(), state.concentration(), state.maxconcentration(), + fractional_flows); + wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); + } + bool well_control_passed = !check_well_controls_; + int well_control_iteration = 0; + do { + // Run solver. + pressure_timer.start(); + std::vector initial_pressure = state.pressure(); + psolver_.solve(timer.currentStepLength(), state, well_state); + + // Renormalize pressure if rock is incompressible, and + // there are no pressure conditions (bcs or wells). + // It is deemed sufficient for now to renormalize + // using geometric volume instead of pore volume. + if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive()) + && allNeumannBCs(bcs_) && allRateWells(wells_)) { + // Compute average pressures of previous and last + // step, and total volume. + double av_prev_press = 0.0; + double av_press = 0.0; + double tot_vol = 0.0; + const int num_cells = grid_.number_of_cells; + for (int cell = 0; cell < num_cells; ++cell) { + av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell]; + av_press += state.pressure()[cell]*grid_.cell_volumes[cell]; + tot_vol += grid_.cell_volumes[cell]; } - if (output_binary_) { - outputStateBinary(grid_, state, timer, output_dir_); + // Renormalization constant + const double ren_const = (av_prev_press - av_press)/tot_vol; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += ren_const; + } + const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells; + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += ren_const; } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); } - // Solve pressure. + // Stop timer and report. + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + + // Optionally, check if well controls are satisfied. if (check_well_controls_) { - computeFractionalFlow(props_, poly_props_, allcells_, - state.saturation(), state.concentration(), state.maxconcentration(), - fractional_flows); - wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); - } - bool well_control_passed = !check_well_controls_; - int well_control_iteration = 0; - do { - // Run solver. - pressure_timer.start(); - std::vector initial_pressure = state.pressure(); - psolver_.solve(timer.currentStepLength(), state, well_state); - - // Renormalize pressure if rock is incompressible, and - // there are no pressure conditions (bcs or wells). - // It is deemed sufficient for now to renormalize - // using geometric volume instead of pore volume. - if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive()) - && allNeumannBCs(bcs_) && allRateWells(wells_)) { - // Compute average pressures of previous and last - // step, and total volume. - double av_prev_press = 0.0; - double av_press = 0.0; - double tot_vol = 0.0; - const int num_cells = grid_.number_of_cells; - for (int cell = 0; cell < num_cells; ++cell) { - av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell]; - av_press += state.pressure()[cell]*grid_.cell_volumes[cell]; - tot_vol += grid_.cell_volumes[cell]; - } - // Renormalization constant - const double ren_const = (av_prev_press - av_press)/tot_vol; - for (int cell = 0; cell < num_cells; ++cell) { - state.pressure()[cell] += ren_const; - } - const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells; - for (int well = 0; well < num_wells; ++well) { - well_state.bhp()[well] += ren_const; - } + Opm::computePhaseFlowRatesPerWell(*wells_, + well_state.perfRates(), + fractional_flows, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { + OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); } - - // Stop timer and report. - pressure_timer.stop(); - double pt = pressure_timer.secsSinceStart(); - std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; - ptime += pt; - - // Optionally, check if well controls are satisfied. - if (check_well_controls_) { - Opm::computePhaseFlowRatesPerWell(*wells_, - well_state.perfRates(), - fractional_flows, - well_resflows_phase); - std::cout << "Checking well conditions." << std::endl; - // For testing we set surface := reservoir - well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); - ++well_control_iteration; - if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { - OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); - } - if (!well_control_passed) { - std::cout << "Well controls not passed, solving again." << std::endl; - } else { - std::cout << "Well conditions met." << std::endl; - } - } - } while (!well_control_passed); - - // Update pore volumes if rock is compressible. - if (rock_comp_props_ && rock_comp_props_->isActive()) { - initial_porevol = porevol; - computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); - } - - // Process transport sources (to include bdy terms and well flows). - Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, - wells_, well_state.perfRates(), transport_src); - - // Find inflow rate. - const double current_time = timer.simulationTimeElapsed(); - double stepsize = timer.currentStepLength(); - polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c); - - // Solve transport. - transport_timer.start(); - if (num_transport_substeps_ != 1) { - stepsize /= double(num_transport_substeps_); - std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; - } - double substep_injected[2] = { 0.0 }; - double substep_produced[2] = { 0.0 }; - double substep_polyinj = 0.0; - double substep_polyprod = 0.0; - injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0; - for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { - tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize, - state.saturation(), state.concentration(), state.maxconcentration()); - Opm::computeInjectedProduced(props_, poly_props_, - state, - transport_src, polymer_inflow_c, stepsize, - substep_injected, substep_produced, substep_polyinj, substep_polyprod); - injected[0] += substep_injected[0]; - injected[1] += substep_injected[1]; - produced[0] += substep_produced[0]; - produced[1] += substep_produced[1]; - polyinj += substep_polyinj; - polyprod += substep_polyprod; - if (use_segregation_split_) { - tsolver_.solveGravity(columns_, &porevol[0], stepsize, - state.saturation(), state.concentration(), state.maxconcentration()); + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; } } - transport_timer.stop(); - double tt = transport_timer.secsSinceStart(); - std::cout << "Transport solver took: " << tt << " seconds." << std::endl; - ttime += tt; + } while (!well_control_passed); - // Report volume balances. - Opm::computeSaturatedVol(porevol, state.saturation(), satvol); - polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); - polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration()); - tot_injected[0] += injected[0]; - tot_injected[1] += injected[1]; - tot_produced[0] += produced[0]; - tot_produced[1] += produced[1]; - tot_polyinj += polyinj; - tot_polyprod += polyprod; - std::cout.precision(5); - const int width = 18; - std::cout << "\nVolume and polymer mass balance: " - " water(pv) oil(pv) polymer(kg)\n"; - std::cout << " Saturated volumes: " - << std::setw(width) << satvol[0]/tot_porevol_init - << std::setw(width) << satvol[1]/tot_porevol_init - << std::setw(width) << polymass << std::endl; - std::cout << " Adsorbed volumes: " - << std::setw(width) << 0.0 - << std::setw(width) << 0.0 - << std::setw(width) << polymass_adsorbed << std::endl; - std::cout << " Injected volumes: " - << std::setw(width) << injected[0]/tot_porevol_init - << std::setw(width) << injected[1]/tot_porevol_init - << std::setw(width) << polyinj << std::endl; - std::cout << " Produced volumes: " - << std::setw(width) << produced[0]/tot_porevol_init - << std::setw(width) << produced[1]/tot_porevol_init - << std::setw(width) << polyprod << std::endl; - std::cout << " Total inj volumes: " - << std::setw(width) << tot_injected[0]/tot_porevol_init - << std::setw(width) << tot_injected[1]/tot_porevol_init - << std::setw(width) << tot_polyinj << std::endl; - std::cout << " Total prod volumes: " - << std::setw(width) << tot_produced[0]/tot_porevol_init - << std::setw(width) << tot_produced[1]/tot_porevol_init - << std::setw(width) << tot_polyprod << std::endl; - std::cout << " In-place + prod - inj: " - << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init - << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init - << std::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl; - std::cout << " Init - now - pr + inj: " - << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init - << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init - << std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed) - << std::endl; - std::cout.precision(8); + // Update pore volumes if rock is compressible. + if (rock_comp_props_ && rock_comp_props_->isActive()) { + initial_porevol = porevol; + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } - watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(), - produced[0]/(produced[0] + produced[1]), - tot_produced[0]/tot_porevol_init); - if (wells_) { - wellreport.push(props_, *wells_, state.saturation(), - timer.simulationTimeElapsed() + timer.currentStepLength(), - well_state.bhp(), well_state.perfRates()); + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, + wells_, well_state.perfRates(), transport_src); + + // Find inflow rate. + const double current_time = timer.simulationTimeElapsed(); + double stepsize = timer.currentStepLength(); + polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c); + + // Solve transport. + transport_timer.start(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + double substep_injected[2] = { 0.0 }; + double substep_produced[2] = { 0.0 }; + double substep_polyinj = 0.0; + double substep_polyprod = 0.0; + injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0; + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize, + state.saturation(), state.concentration(), state.maxconcentration()); + Opm::computeInjectedProduced(props_, poly_props_, + state, + transport_src, polymer_inflow_c, stepsize, + substep_injected, substep_produced, substep_polyinj, substep_polyprod); + injected[0] += substep_injected[0]; + injected[1] += substep_injected[1]; + produced[0] += substep_produced[0]; + produced[1] += substep_produced[1]; + polyinj += substep_polyinj; + polyprod += substep_polyprod; + if (use_segregation_split_) { + tsolver_.solveGravity(columns_, &porevol[0], stepsize, + state.saturation(), state.concentration(), state.maxconcentration()); } } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.saturation(), satvol); + polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); + polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration()); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + tot_polyinj += polyinj; + tot_polyprod += polyprod; + std::cout.precision(5); + const int width = 18; + std::cout << "\nVolume and polymer mass balance: " + " water(pv) oil(pv) polymer(kg)\n"; + std::cout << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init + << std::setw(width) << polymass << std::endl; + std::cout << " Adsorbed volumes: " + << std::setw(width) << 0.0 + << std::setw(width) << 0.0 + << std::setw(width) << polymass_adsorbed << std::endl; + std::cout << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init + << std::setw(width) << polyinj << std::endl; + std::cout << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init + << std::setw(width) << polyprod << std::endl; + std::cout << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init + << std::setw(width) << tot_polyinj << std::endl; + std::cout << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init + << std::setw(width) << tot_polyprod << std::endl; + std::cout << " In-place + prod - inj: " + << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init + << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init + << std::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl; + std::cout << " Init - now - pr + inj: " + << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init + << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init + << std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed) + << std::endl; + std::cout.precision(8); + + watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, state.saturation(), + timer.simulationTimeElapsed() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } if (output_) { if (output_vtk_) { diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp new file mode 100644 index 000000000..ab3e3f38c --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp @@ -0,0 +1,372 @@ +/* + 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 . +*/ + +#ifndef OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED +#define OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +struct UnstructuredGrid; +struct Wells; + +namespace Opm { + + namespace parameter { class ParameterGroup; } + class DerivedGeology; + class RockCompressibility; + class NewtonIterationBlackoilInterface; + class PolymerBlackoilState; + class WellStateFullyImplicitBlackoil; + + + /// A fully implicit solver for the black-oil-polymer problem. + /// + /// The simulator is capable of handling three-phase problems + /// where gas can be dissolved in oil (but not vice versa). It + /// uses an industry-standard TPFA discretization with per-phase + /// upwind weighting of mobilities. + /// + /// It uses automatic differentiation via the class AutoDiffBlock + /// to simplify assembly of the jacobian matrix. + template + class FullyImplicitBlackoilPolymerSolver + { + public: + /// \brief The type of the grid that we use. + typedef T Grid; + /// Construct a solver. It will retain references to the + /// arguments of this functions, and they are expected to + /// remain in scope for the lifetime of the solver. + /// \param[in] param parameters + /// \param[in] grid grid data structure + /// \param[in] fluid fluid properties + /// \param[in] geo rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] wells well structure + /// \param[in] linsolver linear solver + FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param, + const Grid& grid , + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const PolymerPropsAd& polymer_props_ad, + const Wells& wells, + const NewtonIterationBlackoilInterface& linsolver, + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer); + + /// \brief Set threshold pressures that prevent or reduce flow. + /// This prevents flow across faces if the potential + /// difference is less than the threshold. If the potential + /// difference is greater, the threshold value is subtracted + /// before calculating flow. This is treated symmetrically, so + /// flow is prevented or reduced in both directions equally. + /// \param[in] threshold_pressures_by_face array of size equal to the number of faces + /// of the grid passed in the constructor. + void setThresholdPressures(const std::vector& threshold_pressures_by_face); + + /// Take a single forward step, modifiying + /// state.pressure() + /// state.faceflux() + /// state.saturation() + /// state.gasoilratio() + /// wstate.bhp() + /// \param[in] dt time step size + /// \param[in] state reservoir state + /// \param[in] wstate well state + void + step(const double dt , + PolymerBlackoilState& state , + WellStateFullyImplicitBlackoil& wstate, + const std::vector& polymer_inflow); + + private: + // Types and enums + typedef AutoDiffBlock ADB; + typedef ADB::V V; + typedef ADB::M M; + typedef Eigen::Array DataBlock; + + struct ReservoirResidualQuant { + ReservoirResidualQuant(); + std::vector accum; // Accumulations + ADB mflux; // Mass flux (surface conditions) + ADB b; // Reciprocal FVF + ADB head; // Pressure drop across int. interfaces + ADB mob; // Phase mobility (per cell) + }; + + struct SolutionState { + SolutionState(const int np); + ADB pressure; + std::vector saturation; + ADB rs; + ADB rv; + ADB concentration; + ADB qs; + ADB bhp; + }; + + struct WellOps { + WellOps(const Wells& wells); + M w2p; // well -> perf (scatter) + M p2w; // perf -> well (gather) + }; + + enum { Water = BlackoilPropsAdInterface::Water, + Oil = BlackoilPropsAdInterface::Oil , + Gas = BlackoilPropsAdInterface::Gas }; + + // the Newton relaxation type + enum RelaxType { DAMPEN, SOR }; + enum PrimalVariables { Sg = 0, RS = 1, RV = 2 }; + + // Member data + const Grid& grid_; + const BlackoilPropsAdInterface& fluid_; + const DerivedGeology& geo_; + const RockCompressibility* rock_comp_props_; + const PolymerPropsAd& polymer_props_ad_; + const Wells& wells_; + const NewtonIterationBlackoilInterface& linsolver_; + // For each canonical phase -> true if active + const std::vector active_; + // Size = # active phases. Maps active -> canonical phase indices. + const std::vector canph_; + const std::vector cells_; // All grid cells + HelperOps ops_; + const WellOps wops_; + V cmax_; + const bool has_disgas_; + const bool has_vapoil_; + const bool has_polymer_; + const int poly_pos_; + double dp_max_rel_; + double ds_max_; + double drs_max_rel_; + enum RelaxType relax_type_; + double relax_max_; + double relax_increment_; + double relax_rel_tol_; + int max_iter_; + bool use_threshold_pressure_; + V threshold_pressures_by_interior_face_; + + std::vector rq_; + std::vector phaseCondition_; + V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. + + LinearisedBlackoilResidual residual_; + + std::vector primalVariable_; + + // Private methods. + SolutionState + constantState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw); + + SolutionState + variableState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw); + + void + computeAccum(const SolutionState& state, + const int aix ); + + void computeWellConnectionPressures(const SolutionState& state, + const WellStateFullyImplicitBlackoil& xw); + + void + addWellControlEq(const SolutionState& state, + const WellStateFullyImplicitBlackoil& xw, + const V& aliveWells); + + void + addWellEq(const SolutionState& state, + WellStateFullyImplicitBlackoil& xw, + V& aliveWells, + const std::vector& polymer_inflow); + + void updateWellControls(ADB& bhp, + ADB& well_phase_flow_rate, + WellStateFullyImplicitBlackoil& xw) const; + + void + assemble(const V& dtpv, + const PolymerBlackoilState& x, + WellStateFullyImplicitBlackoil& xw, + const std::vector& polymer_inflow); + + V solveJacobianSystem() const; + + void updateState(const V& dx, + PolymerBlackoilState& state, + WellStateFullyImplicitBlackoil& well_state); + + std::vector + computePressures(const SolutionState& state) const; + + std::vector + computeRelPerm(const SolutionState& state) const; + + std::vector + computeRelPermWells(const SolutionState& state, + const DataBlock& well_s, + const std::vector& well_cells) const; + + void + computeMassFlux(const int actph , + const V& transi, + const ADB& kr , + const ADB& p , + const SolutionState& state ); + + void + computeMassFlux(const V& trans, + const std::vector& kr, + const std::vector& phasePressure, + const SolutionState& state); + + void + computeCmax(PolymerBlackoilState& state, + const ADB& c); + + ADB + computeMc(const SolutionState& state) const; + + ADB + rockPorosity(const ADB& p) const; + + ADB + rockPermeability(const ADB& p) const; + void applyThresholdPressures(ADB& dp); + + double + residualNorm() const; + + std::vector residuals() const; + + ADB + fluidViscosity(const int phase, + const ADB& p , + const ADB& rs , + const ADB& rv , + const std::vector& cond, + const std::vector& cells) const; + + ADB + fluidReciprocFVF(const int phase, + const ADB& p , + const ADB& rs , + const ADB& rv , + const std::vector& cond, + const std::vector& cells) const; + + ADB + fluidDensity(const int phase, + const ADB& p , + const ADB& rs , + const ADB& rv , + const std::vector& cond, + const std::vector& cells) const; + + V + fluidRsSat(const V& p, + const V& so, + const std::vector& cells) const; + + ADB + fluidRsSat(const ADB& p, + const ADB& so, + const std::vector& cells) const; + + V + fluidRvSat(const V& p, + const V& so, + const std::vector& cells) const; + + ADB + fluidRvSat(const ADB& p, + const ADB& so, + const std::vector& cells) const; + + + ADB + poroMult(const ADB& p) const; + + ADB + transMult(const ADB& p) const; + + void + classifyCondition(const SolutionState& state, + std::vector& cond ) const; + + const std::vector + phaseCondition() const {return phaseCondition_;} + + void + classifyCondition(const PolymerBlackoilState& state); + + + /// update the primal variable for Sg, Rv or Rs. The Gas phase must + /// be active to call this method. + void + updatePrimalVariableFromState(const PolymerBlackoilState& state); + + /// Update the phaseCondition_ member based on the primalVariable_ member. + void + updatePhaseCondFromPrimalVariable(); + + /// Compute convergence based on total mass balance (tol_mb) and maximum + /// residual mass balance (tol_cnv). + bool getConvergence(const double dt); + + void detectNewtonOscillations(const std::vector>& residual_history, + const int it, const double relaxRelTol, + bool& oscillate, bool& stagnate) const; + + void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const; + + double dpMaxRel() const { return dp_max_rel_; } + double dsMax() const { return ds_max_; } + double drsMaxRel() const { return drs_max_rel_; } + enum RelaxType relaxType() const { return relax_type_; } + double relaxMax() const { return relax_max_; }; + double relaxIncrement() const { return relax_increment_; }; + double relaxRelTol() const { return relax_rel_tol_; }; + double maxIter() const { return max_iter_; } + + }; +} // namespace Opm + +#include "FullyImplicitBlackoilPolymerSolver_impl.hpp" + +#endif // OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp new file mode 100644 index 000000000..100040572 --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp @@ -0,0 +1,2335 @@ +/* + 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 +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +//#include + +// A debugging utility. +#define DUMP(foo) \ + do { \ + std::cout << "==========================================\n" \ + << #foo ":\n" \ + << collapseJacs(foo) << std::endl; \ + } while (0) + +#define DUMPVAL(foo) \ + do { \ + std::cout << "==========================================\n" \ + << #foo ":\n" \ + << foo.value() << std::endl; \ + } while (0) + +#define DISKVAL(foo) \ + do { \ + std::ofstream os(#foo); \ + os.precision(16); \ + os << foo.value() << std::endl; \ + } while (0) + + +namespace Opm { + +typedef AutoDiffBlock ADB; +typedef ADB::V V; +typedef ADB::M M; +typedef Eigen::Array DataBlock; + + +namespace { + + + std::vector + buildAllCells(const int nc) + { + std::vector all_cells(nc); + + for (int c = 0; c < nc; ++c) { all_cells[c] = c; } + + return all_cells; + } + + + + template + V computePerfPress(const Grid& grid, const Wells& wells, const V& rho, const double grav) + { + using namespace Opm::AutoDiffGrid; + const int nw = wells.number_of_wells; + const int nperf = wells.well_connpos[nw]; + const int dim = dimensions(grid); + V wdp = V::Zero(nperf,1); + assert(wdp.size() == rho.size()); + + // Main loop, iterate over all perforations, + // using the following formula: + // wdp(perf) = g*(perf_z - well_ref_z)*rho(perf) + // where the total density rho(perf) is taken to be + // sum_p (rho_p*saturation_p) in the perforation cell. + // [although this is computed on the outside of this function]. + for (int w = 0; w < nw; ++w) { + const double ref_depth = wells.depth_ref[w]; + for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) { + const int cell = wells.well_cells[j]; + const double cell_depth = cellCentroid(grid, cell)[dim - 1]; + wdp[j] = rho[j]*grav*(cell_depth - ref_depth); + } + } + return wdp; + } + + + + template + std::vector + activePhases(const PU& pu) + { + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + std::vector active(maxnp, false); + + for (int p = 0; p < pu.MaxNumPhases; ++p) { + active[ p ] = pu.phase_used[ p ] != 0; + } + + return active; + } + + + + template + std::vector + active2Canonical(const PU& pu) + { + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + std::vector act2can(maxnp, -1); + + for (int phase = 0; phase < maxnp; ++phase) { + if (pu.phase_used[ phase ]) { + act2can[ pu.phase_pos[ phase ] ] = phase; + } + } + + return act2can; + } + + + template + int polymerPos(const PU& pu) + { + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + int pos = 0; + for (int phase = 0; phase < maxnp; ++phase) { + if (pu.phase_used[phase]) { + pos++; + } + } + + return pos; + } +} // Anonymous namespace + + + + template + FullyImplicitBlackoilPolymerSolver:: + FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param, + const Grid& grid , + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const PolymerPropsAd& polymer_props_ad, + const Wells& wells, + const NewtonIterationBlackoilInterface& linsolver, + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer) + : grid_ (grid) + , fluid_ (fluid) + , geo_ (geo) + , rock_comp_props_(rock_comp_props) + , polymer_props_ad_(polymer_props_ad) + , wells_ (wells) + , linsolver_ (linsolver) + , active_(activePhases(fluid.phaseUsage())) + , canph_ (active2Canonical(fluid.phaseUsage())) + , cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid))) + , ops_ (grid) + , wops_ (wells) + , cmax_(V::Zero(Opm::AutoDiffGrid::numCells(grid))) + , has_disgas_(has_disgas) + , has_vapoil_(has_vapoil) + , has_polymer_(has_polymer) + , poly_pos_(polymerPos(fluid.phaseUsage())) + , dp_max_rel_ (1.0e9) + , ds_max_ (0.2) + , drs_max_rel_ (1.0e9) + , relax_type_ (DAMPEN) + , relax_max_ (0.5) + , relax_increment_ (0.1) + , relax_rel_tol_ (0.2) + , max_iter_ (15) + , use_threshold_pressure_(false) + , rq_ (fluid.numPhases()) + , phaseCondition_(AutoDiffGrid::numCells(grid)) + , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), + ADB::null(), + ADB::null() } ) + { + if (has_polymer_) { + if (!active_[Water]) { + OPM_THROW(std::logic_error, "Polymer must solved in water!\n"); + } + // If deck has polymer, residual_ should contain polymer equation. + rq_.resize(fluid_.numPhases()+1); + residual_.material_balance_eq.resize(fluid_.numPhases()+1, ADB::null()); + assert(poly_pos_ == fluid_.numPhases()); + } + dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_); + ds_max_ = param.getDefault("ds_max", ds_max_); + drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_); + relax_max_ = param.getDefault("relax_max", relax_max_); + max_iter_ = param.getDefault("max_iter", max_iter_); + + std::string relaxation_type = param.getDefault("relax_type", std::string("dampen")); + if (relaxation_type == "dampen") { + relax_type_ = DAMPEN; + } else if (relaxation_type == "sor") { + relax_type_ = SOR; + } else { + OPM_THROW(std::runtime_error, "Unknown Relaxtion Type " << relaxation_type); + } + } + + + + template + void + FullyImplicitBlackoilPolymerSolver:: + setThresholdPressures(const std::vector& threshold_pressures_by_face) + { + const int num_faces = AutoDiffGrid::numFaces(grid_); + if (int(threshold_pressures_by_face.size()) != num_faces) { + OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures_by_face input, must be equal to number of faces."); + } + use_threshold_pressure_ = true; + // Map to interior faces. + const int num_ifaces = ops_.internal_faces.size(); + threshold_pressures_by_interior_face_.resize(num_ifaces); + for (int ii = 0; ii < num_ifaces; ++ii) { + threshold_pressures_by_interior_face_[ii] = threshold_pressures_by_face[ops_.internal_faces[ii]]; + } + } + + + + + template + void + FullyImplicitBlackoilPolymerSolver:: + step(const double dt, + PolymerBlackoilState& x , + WellStateFullyImplicitBlackoil& xw, + const std::vector& polymer_inflow) + { + const V pvdt = geo_.poreVolume() / dt; + + if (active_[Gas]) { updatePrimalVariableFromState(x); } + + { + const SolutionState state = constantState(x, xw); + computeCmax(x, state.concentration); + computeAccum(state, 0); + computeWellConnectionPressures(state, xw); + } + + std::vector> residual_history; + + assemble(pvdt, x, xw, polymer_inflow); + + + bool converged = false; + double omega = 1.; + const double r0 = residualNorm(); + + residual_history.push_back(residuals()); + + converged = getConvergence(dt); + + int it = 0; + std::cout << "\nIteration Residual\n" + << std::setw(9) << it << std::setprecision(9) + << std::setw(18) << r0 << std::endl; + + const int sizeNonLinear = residual_.sizeNonLinear(); + + V dxOld = V::Zero(sizeNonLinear); + + bool isOscillate = false; + bool isStagnate = false; + const enum RelaxType relaxtype = relaxType(); + + while ((!converged) && (it < maxIter())) { + V dx = solveJacobianSystem(); + + detectNewtonOscillations(residual_history, it, relaxRelTol(), isOscillate, isStagnate); + + if (isOscillate) { + omega -= relaxIncrement(); + omega = std::max(omega, relaxMax()); + std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl; + } + + stablizeNewton(dx, dxOld, omega, relaxtype); + + updateState(dx, x, xw); + + assemble(pvdt, x, xw, polymer_inflow); + + const double r = residualNorm(); + + residual_history.push_back(residuals()); + + converged = getConvergence(dt); + + it += 1; + std::cout << std::setw(9) << it << std::setprecision(9) + << std::setw(18) << r << std::endl; + } + + if (!converged) { + std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n"; + // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); + } + } + + + + + + template + FullyImplicitBlackoilPolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant() + : accum(2, ADB::null()) + , mflux( ADB::null()) + , b ( ADB::null()) + , head ( ADB::null()) + , mob ( ADB::null()) + { + } + + + + + + template + FullyImplicitBlackoilPolymerSolver::SolutionState::SolutionState(const int np) + : pressure ( ADB::null()) + , saturation(np, ADB::null()) + , rs ( ADB::null()) + , rv ( ADB::null()) + , concentration( ADB::null()) + , qs ( ADB::null()) + , bhp ( ADB::null()) + { + } + + + + + + template + FullyImplicitBlackoilPolymerSolver:: + WellOps::WellOps(const Wells& wells) + : w2p(wells.well_connpos[ wells.number_of_wells ], + wells.number_of_wells) + , p2w(wells.number_of_wells, + wells.well_connpos[ wells.number_of_wells ]) + { + const int nw = wells.number_of_wells; + const int* const wpos = wells.well_connpos; + + typedef Eigen::Triplet Tri; + + std::vector scatter, gather; + scatter.reserve(wpos[nw]); + gather .reserve(wpos[nw]); + + for (int w = 0, i = 0; w < nw; ++w) { + for (; i < wpos[ w + 1 ]; ++i) { + scatter.push_back(Tri(i, w, 1.0)); + gather .push_back(Tri(w, i, 1.0)); + } + } + + w2p.setFromTriplets(scatter.begin(), scatter.end()); + p2w.setFromTriplets(gather .begin(), gather .end()); + } + + + + + + template + typename FullyImplicitBlackoilPolymerSolver::SolutionState + FullyImplicitBlackoilPolymerSolver::constantState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw) + { + auto state = variableState(x, xw); + + // 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()); + state.concentration = ADB::constant(state.concentration.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; + } + + + + + + template + typename FullyImplicitBlackoilPolymerSolver::SolutionState + FullyImplicitBlackoilPolymerSolver::variableState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw) + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + const int np = x.numPhases(); + + std::vector vars0; + // p, Sw and Rs, Rv or Sg, concentration are used as primary depending on solution conditions + vars0.reserve(np + 2); + // Initial pressure. + assert (not x.pressure().empty()); + const V p = Eigen::Map(& x.pressure()[0], nc, 1); + vars0.push_back(p); + + // Initial saturation. + assert (not x.saturation().empty()); + const DataBlock s = Eigen::Map(& x.saturation()[0], nc, np); + const Opm::PhaseUsage pu = fluid_.phaseUsage(); + // We do not handle a Water/Gas situation correctly, guard against it. + assert (active_[ Oil]); + if (active_[ Water ]) { + const V sw = s.col(pu.phase_pos[ Water ]); + vars0.push_back(sw); + } + + // store cell status in vectors + V isRs = V::Zero(nc,1); + V isRv = V::Zero(nc,1); + V isSg = V::Zero(nc,1); + + if (active_[ Gas ]){ + for (int c = 0; c < nc ; c++ ) { + switch (primalVariable_[c]) { + case PrimalVariables::RS: + isRs[c] = 1; + break; + + case PrimalVariables::RV: + isRv[c] = 1; + break; + + default: + isSg[c] = 1; + break; + } + } + + + // define new primary variable xvar depending on solution condition + V xvar(nc); + const V sg = s.col(pu.phase_pos[ Gas ]); + const V rs = Eigen::Map(& x.gasoilratio()[0], x.gasoilratio().size()); + const V rv = Eigen::Map(& x.rv()[0], x.rv().size()); + xvar = isRs*rs + isRv*rv + isSg*sg; + vars0.push_back(xvar); + } + + // Initial polymer concentration. + if (has_polymer_) { + assert (not x.concentration().empty()); + const V c = Eigen::Map(& x.concentration()[0], nc, 1); + vars0.push_back(c); + } + + // Initial well rates. + assert (not xw.wellRates().empty()); + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + 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); + vars0.push_back(qs); + + // Initial well bottom-hole pressure. + assert (not xw.bhp().empty()); + const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + vars0.push_back(bhp); + + std::vector vars = ADB::variables(vars0); + + SolutionState state(np); + + // Pressure. + int nextvar = 0; + state.pressure = vars[ nextvar++ ]; + + // Saturations + const std::vector& bpat = vars[0].blockPattern(); + { + ADB so = ADB::constant(V::Ones(nc, 1), bpat); + + if (active_[ Water ]) { + ADB& sw = vars[ nextvar++ ]; + state.saturation[pu.phase_pos[ Water ]] = sw; + so = so - sw; + } + + if (active_[ Gas]) { + // Define Sg Rs and Rv in terms of xvar. + const ADB& xvar = vars[ nextvar++ ]; + const ADB& sg = isSg*xvar + isRv* so; + state.saturation[ pu.phase_pos[ Gas ] ] = sg; + so = so - sg; + const ADB rsSat = fluidRsSat(state.pressure, so, cells_); + const ADB rvSat = fluidRvSat(state.pressure, so, cells_); + + if (has_disgas_) { + state.rs = (1-isRs) * rsSat + isRs*xvar; + } else { + state.rs = rsSat; + } + if (has_vapoil_) { + state.rv = (1-isRv) * rvSat + isRv*xvar; + } else { + state.rv = rvSat; + } + } + + if (active_[ Oil ]) { + // Note that so is never a primary variable. + state.saturation[ pu.phase_pos[ Oil ] ] = so; + } + } + + // Concentration. + if (has_polymer_) { + state.concentration = vars[nextvar++]; + } + // Qs. + state.qs = vars[ nextvar++ ]; + + // Bhp. + state.bhp = vars[ nextvar++ ]; + + assert(nextvar == int(vars.size())); + + return state; + } + + + + + + template + void + FullyImplicitBlackoilPolymerSolver::computeAccum(const SolutionState& state, + const int aix ) + { + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + + const ADB& press = state.pressure; + const std::vector& sat = state.saturation; + const ADB& rs = state.rs; + const ADB& rv = state.rv; + const ADB& c = state.concentration; + + const std::vector cond = phaseCondition(); + std::vector pressure = computePressures(state); + + const ADB pv_mult = poroMult(press); + + const int maxnp = Opm::BlackoilPhases::MaxNumPhases; + for (int phase = 0; phase < maxnp; ++phase) { + if (active_[ phase ]) { + const int pos = pu.phase_pos[ phase ]; + rq_[pos].b = fluidReciprocFVF(phase, pressure[phase], rs, rv, cond, cells_); + rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos]; + // DUMP(rq_[pos].b); + // DUMP(rq_[pos].accum[aix]); + } + } + + if (active_[ Oil ] && active_[ Gas ]) { + // Account for gas dissolved in oil and vaporized oil + const int po = pu.phase_pos[ Oil ]; + const int pg = pu.phase_pos[ Gas ]; + + rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix]; + rq_[po].accum[aix] += state.rv * rq_[pg].accum[aix]; + //DUMP(rq_[pg].accum[aix]); + } + if (has_polymer_) { + // compute polymer properties. + const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); + const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); + const double rho_rock = polymer_props_ad_.rockDensity(); + const V phi = Eigen::Map(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1); + const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); + // compute total phases and determin polymer position. + rq_[poly_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol) + + pv_mult * rho_rock * (1. - phi) / phi * ads; + } + } + + + + + + template + void FullyImplicitBlackoilPolymerSolver::computeCmax(PolymerBlackoilState& state, + const ADB& c) + { + const int nc = AutoDiffGrid::numCells(grid_); + for (int i = 0; i < nc; ++i) { + cmax_(i) = std::max(cmax_(i), c.value()(i)); + } + std::copy(&cmax_[0], &cmax_[0] + nc, state.maxconcentration().begin()); + } + + + + + + template + void FullyImplicitBlackoilPolymerSolver::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. + 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(); + } + assert(active_[Oil]); + const ADB perf_so = subset(state.saturation[pu.phase_pos[Oil]], well_cells); + 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(), perf_so.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(), perf_so.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 = cellCentroidsZToEigen(grid_); + 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 = dimensions(grid_); + 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); + } + + + + + + template + void + FullyImplicitBlackoilPolymerSolver:: + assemble(const V& pvdt, + const PolymerBlackoilState& x , + WellStateFullyImplicitBlackoil& xw, + const std::vector& polymer_inflow) + { + using namespace Opm::AutoDiffGrid; + // Create the primary variables. + SolutionState state = variableState(x, xw); + + // DISKVAL(state.pressure); + // DISKVAL(state.saturation[0]); + // DISKVAL(state.saturation[1]); + // DISKVAL(state.saturation[2]); + // DISKVAL(state.rs); + // DISKVAL(state.rv); + // DISKVAL(state.qs); + // DISKVAL(state.bhp); + + // -------- Mass balance equations -------- + + // Compute b_p and the accumulation term b_p*s_p for each phase, + // except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o. + // These quantities are stored in rq_[phase].accum[1]. + // The corresponding accumulation terms from the start of + // the timestep (b^0_p*s^0_p etc.) were already computed + // in step() and stored in rq_[phase].accum[0]. + computeAccum(state, 1); + + // Set up the common parts of the mass balance equations + // for each active phase. + const V transi = subset(geo_.transmissibility(), ops_.internal_faces); + const std::vector kr = computeRelPerm(state); + const std::vector pressures = computePressures(state); + computeMassFlux(transi, kr, pressures, state); + for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { + // std::cout << "===== kr[" << phase << "] = \n" << std::endl; + // std::cout << kr[phase]; + // std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl; + // std::cout << rq_[phase].mflux; + residual_.material_balance_eq[ phaseIdx ] = + pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) + + ops_.div*rq_[phaseIdx].mflux; + + + // DUMP(ops_.div*rq_[phase].mflux); + // DUMP(residual_.material_balance_eq[phase]); + } + + // -------- Extra (optional) rs and rv contributions to the mass balance equations -------- + + // Add the extra (flux) terms to the mass balance equations + // From gas dissolved in the oil phase (rs) and oil vaporized in the gas phase (rv) + // The extra terms in the accumulation part of the equation are already handled. + if (active_[ Oil ] && active_[ Gas ]) { + const int po = fluid_.phaseUsage().phase_pos[ Oil ]; + const int pg = fluid_.phaseUsage().phase_pos[ Gas ]; + + const UpwindSelector upwindOil(grid_, ops_, + rq_[po].head.value()); + const ADB rs_face = upwindOil.select(state.rs); + + const UpwindSelector upwindGas(grid_, ops_, + rq_[pg].head.value()); + const ADB rv_face = upwindGas.select(state.rv); + + residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux); + residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux); + + // DUMP(residual_.material_balance_eq[ Gas ]); + + } + + // Add polymer equation. + if (has_polymer_) { + residual_.material_balance_eq[poly_pos_] = pvdt*(rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0]) + + ops_.div*rq_[poly_pos_].mflux; + } + + // Note: updateWellControls() can change all its arguments if + // a well control is switched. + updateWellControls(state.bhp, state.qs, xw); + V aliveWells; + addWellEq(state, xw, aliveWells, polymer_inflow); + addWellControlEq(state, xw, aliveWells); + } + + + + + + template + void FullyImplicitBlackoilPolymerSolver::addWellEq(const SolutionState& state, + WellStateFullyImplicitBlackoil& xw, + V& aliveWells, + const std::vector& polymer_inflow) + { + const int nc = Opm::AutoDiffGrid::numCells(grid_); + 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 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. + const ADB p_perfcell = subset(state.pressure, well_cells); + + // DUMPVAL(p_perfcell); + // DUMPVAL(state.bhp); + // DUMPVAL(ADB::constant(cdp)); + + // 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_.p2w * 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(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]; + wbqt += wbq[phase]; + } + // DUMPVAL(wbqt); + + // check for dead wells + aliveWells = V::Constant(nw, 1.0); + for (int w = 0; w < nw; ++w) { + if (wbqt.value()[w] == 0) { + aliveWells[w] = 0.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) { + const int pos = pu.phase_pos[phase]; + mix_s[phase] = notDeadWells_selector.select(ADB::constant(compi.col(pos)), 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); + } + + // DUMPVAL(ADB::constant(isInjInx)); + // DUMPVAL(ADB::constant(Tw)); + // DUMPVAL(mt); + // DUMPVAL(drawdown); + + // 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); + } + + // DUMPVAL(cqt_i); + // DUMPVAL(volRat); + + // 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; + } + + // DUMPVAL(mix_s[2]); + // DUMPVAL(cq_ps[2]); + + // Add well contributions to mass balance equations + for (int phase = 0; phase < np; ++phase) { + residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc); + } + + // Add well contributions to polymer mass balance equation + if (has_polymer_) { + const ADB mc = computeMc(state); + const V polyin = Eigen::Map(&polymer_inflow[0], nc); + const V poly_in_perf = subset(polyin, well_cells); + const V poly_mc_perf = subset(mc, well_cells).value(); + residual_.material_balance_eq[poly_pos_] -= superset(cq_ps[pu.phase_pos[Water]]*poly_mc_perf + + (wops_.w2p * mix_s[pu.phase_pos[Water]])*cqt_is*poly_in_perf, + well_cells, nc); + } + + // Add WELL EQUATIONS + ADB qs = state.qs; + for (int phase = 0; phase < np; ++phase) { + qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + + } + + + V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); + } + + std::vector cq_d(cq.data(), cq.data() + nperf*np); + xw.perfPhaseRates() = cq_d; + + residual_.well_flux_eq = qs; + } + + + + + + namespace + { + double rateToCompare(const ADB& well_phase_flow_rate, + const int well, + 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) { + // 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; + } + + 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); + + bool broken = false; + + switch (well_type) { + case INJECTOR: + { + switch (ctrl_type) { + case BHP: + broken = bhp.value()[well] > target; + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + broken = rateToCompare(well_phase_flow_rate, + well, num_phases, distr) > target; + break; + } + } + break; + + case PRODUCER: + { + switch (ctrl_type) { + case BHP: + broken = bhp.value()[well] < target; + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + // Note that the rates compared below are negative, + // so breaking the constraints means: too high flow rate + // (as for injection). + broken = rateToCompare(well_phase_flow_rate, + well, num_phases, distr) < target; + break; + } + } + break; + + default: + OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells."); + } + + return broken; + } + } // anonymous namespace + + + + + template + void FullyImplicitBlackoilPolymerSolver::updateWellControls(ADB& bhp, + 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; + bool bhp_changed = false; + bool rates_changed = false; + 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 (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)] << std::endl; + xw.currentControls()[w] = ctrl_index; + // Also updating well state and primary variables. + // We can only be switching to BHP and SURFACE_RATE + // controls since we do not support RESERVOIR_RATE. + const double target = well_controls_iget_target(wc, ctrl_index); + const double* distr = well_controls_iget_distr(wc, ctrl_index); + switch (well_controls_iget_type(wc, ctrl_index)) { + case BHP: + xw.bhp()[w] = target; + bhp_changed = true; + break; + + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + // + // Just record the fact that we need to recompute + // the 'well_phase_flow_rate'. + rates_changed = true; + break; + + case SURFACE_RATE: + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + rates_changed = true; + break; + } + } + } + + // 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()); + } + 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()); + } + } + + + + + template + void FullyImplicitBlackoilPolymerSolver::addWellControlEq(const SolutionState& state, + const WellStateFullyImplicitBlackoil& xw, + const V& aliveWells) + { + const int np = wells_.number_of_phases; + const int nw = wells_.number_of_wells; + + V bhp_targets = V::Zero(nw); + V rate_targets = V::Zero(nw); + M rate_distr(nw, np*nw); + 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]; + + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + bhp_targets (w) = well_controls_iget_target(wc, current); + rate_targets(w) = -1e100; + } + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + // RESERVOIR and SURFACE rates look the same, from a + // high-level point of view, in the system of + // simultaneous linear equations. + + const double* const distr = + well_controls_iget_distr(wc, current); + + for (int p = 0; p < np; ++p) { + rate_distr.insert(w, p*nw + w) = distr[p]; + } + + bhp_targets (w) = -1.0e100; + rate_targets(w) = well_controls_iget_target(wc, current); + } + break; + } + } + 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); + // For wells that are dead (not flowing), and therefore not communicating + // with the reservoir, we set the equation to be equal to the well's total + // flow. This will be a solution only if the target rate is also zero. + M rate_summer(nw, np*nw); + for (int w = 0; w < nw; ++w) { + for (int phase = 0; phase < np; ++phase) { + rate_summer.insert(w, phase*nw + w) = 1.0; + } + } + Selector alive_selector(aliveWells, Selector::NotEqualZero); + residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs); + // DUMP(residual_.well_eq); + } + + + + + + template + V FullyImplicitBlackoilPolymerSolver::solveJacobianSystem() const + { + return linsolver_.computeNewtonIncrement(residual_); + } + + + + + + namespace { + struct Chop01 { + double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); } + }; + } + + + + + + template + void FullyImplicitBlackoilPolymerSolver::updateState(const V& dx, + PolymerBlackoilState& state, + WellStateFullyImplicitBlackoil& well_state) + { + using namespace Opm::AutoDiffGrid; + const int np = fluid_.numPhases(); + const int nc = numCells(grid_); + const int nw = wells_.number_of_wells; + const V null; + assert(null.size() == 0); + const V zero = V::Zero(nc); + const V one = V::Constant(nc, 1.0); + + // store cell status in vectors + V isRs = V::Zero(nc,1); + V isRv = V::Zero(nc,1); + V isSg = V::Zero(nc,1); + if (active_[Gas]) { + for (int c = 0; c < nc; ++c) { + switch (primalVariable_[c]) { + case PrimalVariables::RS: + isRs[c] = 1; + break; + + case PrimalVariables::RV: + isRv[c] = 1; + break; + + default: + isSg[c] = 1; + break; + } + } + } + + // Extract parts of dx corresponding to each part. + const V dp = subset(dx, Span(nc)); + int varstart = nc; + const V dsw = active_[Water] ? subset(dx, Span(nc, 1, varstart)) : null; + varstart += dsw.size(); + + const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null; + varstart += dxvar.size(); + + const V dc = (has_polymer_) ? subset(dx, Span(nc, 1, varstart)) : null; + varstart += dc.size(); + + const V dqs = subset(dx, Span(np*nw, 1, varstart)); + varstart += dqs.size(); + const V dbhp = subset(dx, Span(nw, 1, varstart)); + varstart += dbhp.size(); + assert(varstart == dx.size()); + + // Pressure update. + const double dpmaxrel = dpMaxRel(); + const V p_old = Eigen::Map(&state.pressure()[0], nc, 1); + const V absdpmax = dpmaxrel*p_old.abs(); + const V dp_limited = sign(dp) * dp.abs().min(absdpmax); + const V p = (p_old - dp_limited).max(zero); + std::copy(&p[0], &p[0] + nc, state.pressure().begin()); + + + // Saturation updates. + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + const DataBlock s_old = Eigen::Map(& state.saturation()[0], nc, np); + const double dsmax = dsMax(); + V so; + V sw; + V sg; + + { + V maxVal = zero; + V dso = zero; + if (active_[Water]){ + maxVal = dsw.abs().max(maxVal); + dso = dso - dsw; + } + + V dsg; + if (active_[Gas]){ + dsg = isSg * dxvar - isRv * dsw; + maxVal = dsg.abs().max(maxVal); + dso = dso - dsg; + } + + maxVal = dso.abs().max(maxVal); + + V step = dsmax/maxVal; + step = step.min(1.); + + if (active_[Water]) { + const int pos = pu.phase_pos[ Water ]; + const V sw_old = s_old.col(pos); + sw = sw_old - step * dsw; + } + + if (active_[Gas]) { + const int pos = pu.phase_pos[ Gas ]; + const V sg_old = s_old.col(pos); + sg = sg_old - step * dsg; + } + + const int pos = pu.phase_pos[ Oil ]; + const V so_old = s_old.col(pos); + so = so_old - step * dso; + } + + // Appleyard chop process. + auto ixg = sg < 0; + for (int c = 0; c < nc; ++c) { + if (ixg[c]) { + sw[c] = sw[c] / (1-sg[c]); + so[c] = so[c] / (1-sg[c]); + sg[c] = 0; + } + } + + + auto ixo = so < 0; + for (int c = 0; c < nc; ++c) { + if (ixo[c]) { + sw[c] = sw[c] / (1-so[c]); + sg[c] = sg[c] / (1-so[c]); + so[c] = 0; + } + } + + auto ixw = sw < 0; + for (int c = 0; c < nc; ++c) { + if (ixw[c]) { + so[c] = so[c] / (1-sw[c]); + sg[c] = sg[c] / (1-sw[c]); + sw[c] = 0; + } + } + + const V sumSat = sw + so + sg; + sw = sw / sumSat; + so = so / sumSat; + sg = sg / sumSat; + + // Update the state + for (int c = 0; c < nc; ++c) { + state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c]; + } + + for (int c = 0; c < nc; ++c) { + state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c]; + } + + if (active_[ Oil ]) { + const int pos = pu.phase_pos[ Oil ]; + for (int c = 0; c < nc; ++c) { + state.saturation()[c*np + pos] = so[c]; + } + } + + // Update rs and rv + const double drsmaxrel = drsMaxRel(); + const double drvmax = 1e9;//% same as in Mrst + V rs; + if (has_disgas_) { + const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); + const V drs = isRs * dxvar; + const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel); + rs = rs_old - drs_limited; + } + V rv; + if (has_vapoil_) { + const V rv_old = Eigen::Map(&state.rv()[0], nc); + const V drv = isRv * dxvar; + const V drv_limited = sign(drv) * drv.abs().min(drvmax); + rv = rv_old - drv_limited; + } + + // Update the state + if (has_disgas_) + std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin()); + + if (has_vapoil_) + std::copy(&rv[0], &rv[0] + nc, state.rv().begin()); + + // Sg is used as primal variable for water only cells. + const double epsilon = std::sqrt(std::numeric_limits::epsilon()); + auto watOnly = sw > (1 - epsilon); + + // phase translation sg <-> rs + const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_); + const V rsSat = fluidRsSat(p, so, cells_); + + std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg); + + if (has_disgas_) { + // The obvious case + auto hasGas = (sg > 0 && isRs == 0); + + // keep oil saturated if previous sg is sufficient large: + const int pos = pu.phase_pos[ Gas ]; + auto hadGas = (sg <= 0 && s_old.col(pos) > epsilon); + // Set oil saturated if previous rs is sufficiently large + const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); + auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) ); + + auto useSg = watOnly || hasGas || hadGas || gasVaporized; + for (int c = 0; c < nc; ++c) { + if (useSg[c]) { rs[c] = rsSat[c];} + else { primalVariable_[c] = PrimalVariables::RS; } + + } + } + + // phase transitions so <-> rv + const V rvSat0 = fluidRvSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_); + const V rvSat = fluidRvSat(p, so, cells_); + + if (has_vapoil_) { + // The obvious case + auto hasOil = (so > 0 && isRv == 0); + + // keep oil saturated if previous so is sufficient large: + const int pos = pu.phase_pos[ Oil ]; + auto hadOil = (so <= 0 && s_old.col(pos) > epsilon ); + // Set oil saturated if previous rv is sufficiently large + const V rv_old = Eigen::Map(&state.rv()[0], nc); + auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) ); + auto useSg = watOnly || hasOil || hadOil || oilCondensed; + for (int c = 0; c < nc; ++c) { + if (useSg[c]) { rv[c] = rvSat[c]; } + else {primalVariable_[c] = PrimalVariables::RV; } + + } + + } + + //Polymer concentration updates. + if (has_polymer_) { + const V c_old = Eigen::Map(&state.concentration()[0], nc, 1); + const V c = (c_old - dc).max(zero); + std::copy(&c[0], &c[0] + nc, state.concentration().begin()); + } + + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const V dwr = Eigen::Map(wwr.data(), nw*np); + const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const V wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + + // Bhp update. + const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); + const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); + const V bhp = bhp_old - dbhp_limited; + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + + // Update phase conditions used for property calculations. + updatePhaseCondFromPrimalVariable(); + } + + + + + + template + std::vector + FullyImplicitBlackoilPolymerSolver::computeRelPerm(const SolutionState& state) const + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + const std::vector& bpat = state.pressure.blockPattern(); + + const ADB null = ADB::constant(V::Zero(nc, 1), bpat); + + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + const ADB sw = (active_[ Water ] + ? state.saturation[ pu.phase_pos[ Water ] ] + : null); + + const ADB so = (active_[ Oil ] + ? state.saturation[ pu.phase_pos[ Oil ] ] + : null); + + const ADB sg = (active_[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : null); + + return fluid_.relperm(sw, so, sg, cells_); + } + + + template + std::vector + FullyImplicitBlackoilPolymerSolver::computePressures(const SolutionState& state) const + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + const std::vector& bpat = state.pressure.blockPattern(); + + const ADB null = ADB::constant(V::Zero(nc, 1), bpat); + + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + const ADB sw = (active_[ Water ] + ? state.saturation[ pu.phase_pos[ Water ] ] + : null); + + const ADB so = (active_[ Oil ] + ? state.saturation[ pu.phase_pos[ Oil ] ] + : null); + + const ADB sg = (active_[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : null); + + // convert the pressure offsets to the capillary pressures + std::vector pressure = fluid_.capPress(sw, so, sg, cells_); + for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) { + // The reference pressure is always the liquid phase (oil) pressure. + if (phaseIdx == BlackoilPhases::Liquid) + continue; + pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid]; + } + + // Since pcow = po - pw, but pcog = pg - po, + // we have + // pw = po - pcow + // pg = po + pcgo + // This is an unfortunate inconsistency, but a convention we must handle. + for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) { + if (phaseIdx == BlackoilPhases::Aqua) { + pressure[phaseIdx] = state.pressure - pressure[phaseIdx]; + } else { + pressure[phaseIdx] += state.pressure; + } + } + + return pressure; + } + + + + + template + std::vector + FullyImplicitBlackoilPolymerSolver::computeRelPermWells(const SolutionState& state, + const DataBlock& well_s, + const std::vector& well_cells) const + { + const int nw = wells_.number_of_wells; + const int nperf = wells_.well_connpos[nw]; + const std::vector& bpat = state.pressure.blockPattern(); + + const ADB null = ADB::constant(V::Zero(nperf), bpat); + + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + const ADB sw = (active_[ Water ] + ? ADB::constant(well_s.col(pu.phase_pos[ Water ]), bpat) + : null); + + const ADB so = (active_[ Oil ] + ? ADB::constant(well_s.col(pu.phase_pos[ Oil ]), bpat) + : null); + + const ADB sg = (active_[ Gas ] + ? ADB::constant(well_s.col(pu.phase_pos[ Gas ]), bpat) + : null); + + return fluid_.relperm(sw, so, sg, well_cells); + } + + + + + + template + void + FullyImplicitBlackoilPolymerSolver::computeMassFlux(const V& transi, + const std::vector& kr , + const std::vector& phasePressure, + const SolutionState& state) + { + for (int phase = 0; phase < fluid_.numPhases(); ++phase) { + const int canonicalPhaseIdx = canph_[phase]; + const std::vector cond = phaseCondition(); + + const ADB tr_mult = transMult(state.pressure); + const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_); + rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu; + + const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.rs, state.rv, cond, cells_); + + ADB& head = rq_[phase].head; + + // compute gravity potensial using the face average as in eclipse and MRST + const ADB rhoavg = ops_.caver * rho; + + ADB dp = ops_.ngrad * phasePressure[canonicalPhaseIdx] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix())); + + if (use_threshold_pressure_) { + applyThresholdPressures(dp); + } + + head = transi*dp; + if (canonicalPhaseIdx == Water) { + if(has_polymer_) { + const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); + const ADB mc = computeMc(state); + ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, + cmax, + kr[canonicalPhaseIdx], + state.saturation[canonicalPhaseIdx]); + ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value().data()); + rq_[phase].mob = tr_mult * krw_eff * inv_wat_eff_visc; + rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_wat_eff_visc; + rq_[poly_pos_].b = rq_[phase].b; + rq_[poly_pos_].head = rq_[phase].head; + UpwindSelector upwind(grid_, ops_, rq_[poly_pos_].head.value()); + rq_[poly_pos_].mflux = upwind.select(rq_[poly_pos_].b * rq_[poly_pos_].mob) * rq_[poly_pos_].head; + } + } + //head = transi*(ops_.ngrad * phasePressure) + gflux; + + UpwindSelector upwind(grid_, ops_, head.value()); + + const ADB& b = rq_[phase].b; + const ADB& mob = rq_[phase].mob; + rq_[phase].mflux = upwind.select(b * mob) * head; + } + } + + + + template + void + FullyImplicitBlackoilPolymerSolver::applyThresholdPressures(ADB& dp) + { + // We support reversible threshold pressures only. + // Method: if the potential difference is lower (in absolute + // value) than the threshold for any face, then the potential + // (and derivatives) is set to zero. If it is above the + // threshold, the threshold pressure is subtracted from the + // absolute potential (the potential is moved towards zero). + + // Identify the set of faces where the potential is under the + // threshold, that shall have zero flow. Storing the bool + // Array as a V (a double Array) with 1 and 0 elements, a + // 1 where flow is allowed, a 0 where it is not. + const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast(); + + // Create a sparse vector that nullifies the low potential elements. + const M keep_high_potential = spdiag(high_potential); + + // Find the current sign for the threshold modification + const V sign_dp = sign(dp.value()); + const V threshold_modification = sign_dp * threshold_pressures_by_interior_face_; + + // Modify potential and nullify where appropriate. + dp = keep_high_potential * (dp - threshold_modification); + } + + + + + + template + double + FullyImplicitBlackoilPolymerSolver::residualNorm() const + { + double globalNorm = 0; + std::vector::const_iterator quantityIt = residual_.material_balance_eq.begin(); + const std::vector::const_iterator endQuantityIt = residual_.material_balance_eq.end(); + for (; quantityIt != endQuantityIt; ++quantityIt) { + const double quantityResid = (*quantityIt).value().matrix().norm(); + if (!std::isfinite(quantityResid)) { + const int trouble_phase = quantityIt - residual_.material_balance_eq.begin(); + OPM_THROW(Opm::NumericalProblem, + "Encountered a non-finite residual in material balance equation " + << trouble_phase); + } + globalNorm = std::max(globalNorm, quantityResid); + } + globalNorm = std::max(globalNorm, residual_.well_flux_eq.value().matrix().norm()); + globalNorm = std::max(globalNorm, residual_.well_eq.value().matrix().norm()); + + return globalNorm; + } + + + template + std::vector + FullyImplicitBlackoilPolymerSolver::residuals() const + { + std::vector residual; + + std::vector::const_iterator massBalanceIt = residual_.material_balance_eq.begin(); + const std::vector::const_iterator endMassBalanceIt = residual_.material_balance_eq.end(); + + for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) { + const double massBalanceResid = (*massBalanceIt).value().matrix().template lpNorm(); + if (!std::isfinite(massBalanceResid)) { + OPM_THROW(Opm::NumericalProblem, + "Encountered a non-finite residual"); + } + residual.push_back(massBalanceResid); + } + + // the following residuals are not used in the oscillation detection now + const double wellFluxResid = residual_.well_flux_eq.value().matrix().template lpNorm(); + if (!std::isfinite(wellFluxResid)) { + OPM_THROW(Opm::NumericalProblem, + "Encountered a non-finite residual"); + } + residual.push_back(wellFluxResid); + + const double wellResid = residual_.well_eq.value().matrix().template lpNorm(); + if (!std::isfinite(wellResid)) { + OPM_THROW(Opm::NumericalProblem, + "Encountered a non-finite residual"); + } + residual.push_back(wellResid); + + return residual; + } + + template + void + FullyImplicitBlackoilPolymerSolver::detectNewtonOscillations(const std::vector>& residual_history, + const int it, const double relaxRelTol, + bool& oscillate, bool& stagnate) const + { + // The detection of oscillation in two primary variable results in the report of the detection + // of oscillation for the solver. + // Only the saturations are used for oscillation detection for the black oil model. + // Stagnate is not used for any treatment here. + + if ( it < 2 ) { + oscillate = false; + stagnate = false; + return; + } + + int oscillatePhase = 0; + const std::vector& F0 = residual_history[it]; + const std::vector& F1 = residual_history[it - 1]; + const std::vector& F2 = residual_history[it - 2]; + for (int p= 0; p < fluid_.numPhases(); ++p){ + const double d1 = std::abs((F0[p] - F2[p]) / F0[p]); + const double d2 = std::abs((F0[p] - F1[p]) / F0[p]); + + oscillatePhase += (d1 < relaxRelTol) && (relaxRelTol < d2); + + stagnate = ! (std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3); + } + + oscillate = (oscillatePhase > 1); + } + + + template + void + FullyImplicitBlackoilPolymerSolver::stablizeNewton(V& dx, V& dxOld, const double omega, + const RelaxType relax_type) const + { + // The dxOld is updated with dx. + // If omega is equal to 1., no relaxtion will be appiled. + + const V tempDxOld = dxOld; + dxOld = dx; + + switch (relax_type) { + case DAMPEN: + if (omega == 1.) { + return; + } + dx = dx*omega; + return; + case SOR: + if (omega == 1.) { + return; + } + dx = dx*omega + (1.-omega)*tempDxOld; + return; + default: + OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type."); + } + + return; + } + + template + bool + FullyImplicitBlackoilPolymerSolver::getConvergence(const double dt) + { + const double tol_mb = 1.0e-7; + const double tol_cnv = 1.0e-3; + + const int nc = Opm::AutoDiffGrid::numCells(grid_); + const Opm::PhaseUsage& pu = fluid_.phaseUsage(); + + const V pv = geo_.poreVolume(); + const double pvSum = pv.sum(); + + const std::vector cond = phaseCondition(); + + double CNVW = 0.; + double CNVO = 0.; + double CNVG = 0.; + + double RW_sum = 0.; + double RO_sum = 0.; + double RG_sum = 0.; + + double BW_avg = 0.; + double BO_avg = 0.; + double BG_avg = 0.; + + if (active_[Water]) { + const int pos = pu.phase_pos[Water]; + const ADB& tempBW = rq_[pos].b; + V BW = 1./tempBW.value(); + V RW = residual_.material_balance_eq[Water].value(); + BW_avg = BW.sum()/nc; + const V tempV = RW.abs()/pv; + + CNVW = BW_avg * dt * tempV.maxCoeff(); + RW_sum = RW.sum(); + } + + if (active_[Oil]) { + // Omit the disgas here. We should add it. + const int pos = pu.phase_pos[Oil]; + const ADB& tempBO = rq_[pos].b; + V BO = 1./tempBO.value(); + V RO = residual_.material_balance_eq[Oil].value(); + BO_avg = BO.sum()/nc; + const V tempV = RO.abs()/pv; + + CNVO = BO_avg * dt * tempV.maxCoeff(); + RO_sum = RO.sum(); + } + + if (active_[Gas]) { + // Omit the vapoil here. We should add it. + const int pos = pu.phase_pos[Gas]; + const ADB& tempBG = rq_[pos].b; + V BG = 1./tempBG.value(); + V RG = residual_.material_balance_eq[Gas].value(); + BG_avg = BG.sum()/nc; + const V tempV = RG.abs()/pv; + + CNVG = BG_avg * dt * tempV.maxCoeff(); + RG_sum = RG.sum(); + } + + double tempValue = tol_mb * pvSum /dt; + + bool converged_MB = (fabs(BW_avg*RW_sum) < tempValue) + && (fabs(BO_avg*RO_sum) < tempValue) + && (fabs(BG_avg*RG_sum) < tempValue); + + bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv); + + double residualWellFlux = residual_.well_flux_eq.value().matrix().template lpNorm(); + double residualWell = residual_.well_eq.value().matrix().template lpNorm(); + + bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa); + + bool converged = converged_MB && converged_CNV && converged_Well; + +#ifdef OPM_VERBOSE + std::cout << " CNVW " << CNVW << " CNVO " << CNVO << " CNVG " << CNVG << std::endl; + std::cout << " converged_MB " << converged_MB << " converged_CNV " << converged_CNV + << " converged_Well " << converged_Well << " converged " << converged << std::endl; +#endif + return converged; + } + + + template + ADB + FullyImplicitBlackoilPolymerSolver::fluidViscosity(const int phase, + const ADB& p , + const ADB& rs , + const ADB& rv , + const std::vector& cond, + const std::vector& cells) const + { + switch (phase) { + case Water: + return fluid_.muWat(p, cells); + case Oil: { + return fluid_.muOil(p, rs, cond, cells); + } + case Gas: + return fluid_.muGas(p, rv, cond, cells); + default: + OPM_THROW(std::runtime_error, "Unknown phase index " << phase); + } + } + + + + + + template + ADB + FullyImplicitBlackoilPolymerSolver::fluidReciprocFVF(const int phase, + const ADB& p , + const ADB& rs , + const ADB& rv , + const std::vector& cond, + const std::vector& cells) const + { + switch (phase) { + case Water: + return fluid_.bWat(p, cells); + case Oil: { + return fluid_.bOil(p, rs, cond, cells); + } + case Gas: + return fluid_.bGas(p, rv, cond, cells); + default: + OPM_THROW(std::runtime_error, "Unknown phase index " << phase); + } + } + + + + + + template + ADB + FullyImplicitBlackoilPolymerSolver::fluidDensity(const int phase, + const ADB& p , + const ADB& rs , + const ADB& rv , + const std::vector& cond, + const std::vector& cells) const + { + const double* rhos = fluid_.surfaceDensity(); + ADB b = fluidReciprocFVF(phase, p, rs, rv, cond, cells); + ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b; + if (phase == Oil && active_[Gas]) { + // It is correct to index into rhos with canonical phase indices. + rho += V::Constant(p.size(), 1, rhos[Gas]) * rs * b; + } + if (phase == Gas && active_[Oil]) { + // It is correct to index into rhos with canonical phase indices. + rho += V::Constant(p.size(), 1, rhos[Oil]) * rv * b; + } + return rho; + } + + + + + + template + V + FullyImplicitBlackoilPolymerSolver::fluidRsSat(const V& p, + const V& satOil, + const std::vector& cells) const + { + return fluid_.rsSat(p, satOil, cells); + } + + + + + + template + ADB + FullyImplicitBlackoilPolymerSolver::fluidRsSat(const ADB& p, + const ADB& satOil, + const std::vector& cells) const + { + return fluid_.rsSat(p, satOil, cells); + } + + template + V + FullyImplicitBlackoilPolymerSolver::fluidRvSat(const V& p, + const V& satOil, + const std::vector& cells) const + { + return fluid_.rvSat(p, satOil, cells); + } + + + + + + template + ADB + FullyImplicitBlackoilPolymerSolver::fluidRvSat(const ADB& p, + const ADB& satOil, + const std::vector& cells) const + { + return fluid_.rvSat(p, satOil, cells); + } + + + + + + template + ADB + FullyImplicitBlackoilPolymerSolver::computeMc(const SolutionState& state) const + { + return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration); + } + + + + + template + ADB + FullyImplicitBlackoilPolymerSolver::poroMult(const ADB& p) const + { + const int n = p.size(); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + V pm(n); + V dpm(n); + for (int i = 0; i < n; ++i) { + pm[i] = rock_comp_props_->poroMult(p.value()[i]); + dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]); + } + ADB::M dpm_diag = spdiag(dpm); + const int num_blocks = p.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dpm_diag * p.derivative()[block]; + } + return ADB::function(pm, jacs); + } else { + return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); + } + } + + + + + + template + ADB + FullyImplicitBlackoilPolymerSolver::transMult(const ADB& p) const + { + const int n = p.size(); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + V tm(n); + V dtm(n); + for (int i = 0; i < n; ++i) { + tm[i] = rock_comp_props_->transMult(p.value()[i]); + dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]); + } + ADB::M dtm_diag = spdiag(dtm); + const int num_blocks = p.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dtm_diag * p.derivative()[block]; + } + return ADB::function(tm, jacs); + } else { + return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); + } + } + + + /* + template + void + FullyImplicitBlackoilPolymerSolver:: + classifyCondition(const SolutionState& state, + std::vector& cond ) const + { + const PhaseUsage& pu = fluid_.phaseUsage(); + + if (active_[ Gas ]) { + // Oil/Gas or Water/Oil/Gas system + const int po = pu.phase_pos[ Oil ]; + const int pg = pu.phase_pos[ Gas ]; + + const V& so = state.saturation[ po ].value(); + const V& sg = state.saturation[ pg ].value(); + + cond.resize(sg.size()); + + for (V::Index c = 0, e = sg.size(); c != e; ++c) { + if (so[c] > 0) { cond[c].setFreeOil (); } + if (sg[c] > 0) { cond[c].setFreeGas (); } + if (active_[ Water ]) { cond[c].setFreeWater(); } + } + } + else { + // Water/Oil system + assert (active_[ Water ]); + + const int po = pu.phase_pos[ Oil ]; + const V& so = state.saturation[ po ].value(); + + cond.resize(so.size()); + + for (V::Index c = 0, e = so.size(); c != e; ++c) { + cond[c].setFreeWater(); + + if (so[c] > 0) { cond[c].setFreeOil(); } + } + } + } */ + + + template + void + FullyImplicitBlackoilPolymerSolver::classifyCondition(const PolymerBlackoilState& state) + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + const int np = state.numPhases(); + + const PhaseUsage& pu = fluid_.phaseUsage(); + const DataBlock s = Eigen::Map(& state.saturation()[0], nc, np); + if (active_[ Gas ]) { + // Oil/Gas or Water/Oil/Gas system + const V so = s.col(pu.phase_pos[ Oil ]); + const V sg = s.col(pu.phase_pos[ Gas ]); + + for (V::Index c = 0, e = sg.size(); c != e; ++c) { + if (so[c] > 0) { phaseCondition_[c].setFreeOil (); } + if (sg[c] > 0) { phaseCondition_[c].setFreeGas (); } + if (active_[ Water ]) { phaseCondition_[c].setFreeWater(); } + } + } + else { + // Water/Oil system + assert (active_[ Water ]); + + const V so = s.col(pu.phase_pos[ Oil ]); + + + for (V::Index c = 0, e = so.size(); c != e; ++c) { + phaseCondition_[c].setFreeWater(); + + if (so[c] > 0) { phaseCondition_[c].setFreeOil(); } + } + } + + + } + + template + void + FullyImplicitBlackoilPolymerSolver::updatePrimalVariableFromState(const PolymerBlackoilState& state) + { + using namespace Opm::AutoDiffGrid; + const int nc = numCells(grid_); + const int np = state.numPhases(); + + const PhaseUsage& pu = fluid_.phaseUsage(); + const DataBlock s = Eigen::Map(& state.saturation()[0], nc, np); + + // Water/Oil/Gas system + assert (active_[ Gas ]); + + // reset the primary variables if RV and RS is not set Sg is used as primary variable. + primalVariable_.resize(nc); + std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg); + + const V sg = s.col(pu.phase_pos[ Gas ]); + const V so = s.col(pu.phase_pos[ Oil ]); + const V sw = s.col(pu.phase_pos[ Water ]); + + const double epsilon = std::sqrt(std::numeric_limits::epsilon()); + auto watOnly = sw > (1 - epsilon); + auto hasOil = so > 0; + auto hasGas = sg > 0; + + // For oil only cells Rs is used as primal variable. For cells almost full of water + // the default primal variable (Sg) is used. + if (has_disgas_) { + for (V::Index c = 0, e = sg.size(); c != e; ++c) { + if ( !watOnly[c] && hasOil[c] && !hasGas[c] ) {primalVariable_[c] = PrimalVariables::RS; } + } + } + + // For gas only cells Rv is used as primal variable. For cells almost full of water + // the default primal variable (Sg) is used. + if (has_vapoil_) { + for (V::Index c = 0, e = so.size(); c != e; ++c) { + if ( !watOnly[c] && hasGas[c] && !hasOil[c] ) {primalVariable_[c] = PrimalVariables::RV; } + } + } + updatePhaseCondFromPrimalVariable(); + } + + + + + + /// Update the phaseCondition_ member based on the primalVariable_ member. + template + void + FullyImplicitBlackoilPolymerSolver::updatePhaseCondFromPrimalVariable() + { + if (! active_[Gas]) { + OPM_THROW(std::logic_error, "updatePhaseCondFromPrimarVariable() logic requires active gas phase."); + } + const int nc = primalVariable_.size(); + for (int c = 0; c < nc; ++c) { + phaseCondition_[c] = PhasePresence(); // No free phases. + phaseCondition_[c].setFreeWater(); // Not necessary for property calculation usage. + switch (primalVariable_[c]) { + case PrimalVariables::Sg: + phaseCondition_[c].setFreeOil(); + phaseCondition_[c].setFreeGas(); + break; + case PrimalVariables::RS: + phaseCondition_[c].setFreeOil(); + break; + case PrimalVariables::RV: + phaseCondition_[c].setFreeGas(); + break; + default: + OPM_THROW(std::logic_error, "Unknown primary variable enum value in cell " << c << ": " << primalVariable_[c]); + } + } + } + + + + +} // namespace Opm diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp new file mode 100644 index 000000000..3cfbedb01 --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -0,0 +1,1029 @@ +/* + Copyright 2014 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 + + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// A debugging utility. +#define DUMP(foo) \ + do { \ + std::cout << "==========================================\n" \ + << #foo ":\n" \ + << collapseJacs(foo) << std::endl; \ + } while (0) + + + +namespace Opm { + +typedef AutoDiffBlock ADB; +typedef ADB::V V; +typedef ADB::M M; +typedef Eigen::Array DataBlock; + + + + + +namespace { + + std::vector + buildAllCells(const int nc) + { + std::vector all_cells(nc); + + for (int c = 0; c < nc; ++c) { all_cells[c] = c; } + + return all_cells; + } + + + + + + template + AutoDiffBlock::M + gravityOperator(const UnstructuredGrid& grid, + const HelperOps& ops , + const GeoProps& geo ) + { + const int nc = grid.number_of_cells; + + std::vector f2hf(2 * grid.number_of_faces, -1); + for (int c = 0, i = 0; c < nc; ++c) { + for (; i < grid.cell_facepos[c + 1]; ++i) { + const int f = grid.cell_faces[ i ]; + const int p = 0 + (grid.face_cells[2*f + 0] != c); + + f2hf[2*f + p] = i; + } + } + + typedef AutoDiffBlock::V V; + typedef AutoDiffBlock::M M; + + const V& gpot = geo.gravityPotential(); + const V& trans = geo.transmissibility(); + + const HelperOps::IFaces::Index ni = ops.internal_faces.size(); + + typedef Eigen::Triplet Tri; + std::vector grav; grav.reserve(2 * ni); + for (HelperOps::IFaces::Index i = 0; i < ni; ++i) { + const int f = ops.internal_faces[ i ]; + const int c1 = grid.face_cells[2*f + 0]; + const int c2 = grid.face_cells[2*f + 1]; + + assert ((c1 >= 0) && (c2 >= 0)); + + const double dG1 = gpot[ f2hf[2*f + 0] ]; + const double dG2 = gpot[ f2hf[2*f + 1] ]; + const double t = trans[ f ]; + + grav.push_back(Tri(i, c1, t * dG1)); + grav.push_back(Tri(i, c2, - t * dG2)); + } + + M G(ni, nc); G.setFromTriplets(grav.begin(), grav.end()); + + return G; + } + + + + + + V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav) + { + const int nw = wells.number_of_wells; + const int nperf = wells.well_connpos[nw]; + const int dim = grid.dimensions; + V wdp = V::Zero(nperf,1); + assert(wdp.size() == rho.size()); + + // Main loop, iterate over all perforations, + // using the following formula: + // wdp(perf) = g*(perf_z - well_ref_z)*rho(perf) + // where the total density rho(perf) is taken to be + // sum_p (rho_p*saturation_p) in the perforation cell. + // [although this is computed on the outside of this function]. + for (int w = 0; w < nw; ++w) { + const double ref_depth = wells.depth_ref[w]; + for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) { + const int cell = wells.well_cells[j]; + const double cell_depth = grid.cell_centroids[dim * cell + dim - 1]; + wdp[j] = rho[j]*grav*(cell_depth - ref_depth); + } + } + return wdp; + } + +} // Anonymous namespace + + + + + + FullyImplicitCompressiblePolymerSolver:: + FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid& grid, + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const PolymerPropsAd& polymer_props_ad, + const Wells& wells, + const NewtonIterationBlackoilInterface& linsolver) + : grid_ (grid) + , fluid_ (fluid) + , geo_ (geo) + , rock_comp_props_(rock_comp_props) + , polymer_props_ad_(polymer_props_ad) + , wells_ (wells) + , linsolver_ (linsolver) + , cells_ (buildAllCells(grid.number_of_cells)) + , ops_ (grid) + , wops_ (wells) + , grav_ (gravityOperator(grid_, ops_, geo_)) + , cmax_(V::Zero(grid.number_of_cells)) + , phaseCondition_ (grid.number_of_cells) + , rq_ (fluid.numPhases() + 1) + , residual_ ( { std::vector(fluid.numPhases() + 1, ADB::null()), + ADB::null(), + ADB::null() } ) + { + } + + + + + + void + FullyImplicitCompressiblePolymerSolver:: + step(const double dt, + PolymerBlackoilState& x , + WellStateFullyImplicitBlackoil& xw, + const std::vector& polymer_inflow) + { + + const SolutionState state = constantState(x, xw); + computeCmax(x, state.concentration); + computeAccum(state, 0); + + const double atol = 1.0e-12; + const double rtol = 5.0e-8; + const int maxit = 15; + assemble(dt, x, xw, polymer_inflow); + + const double r0 = residualNorm(); + const double r_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm(); + int it = 0; + std::cout << "\nIteration Residual Polymer Res\n" + << std::setw(9) << it << std::setprecision(9) + << std::setw(18) << r0 << std::setprecision(9) + << std::setw(18) << r_polymer << std::endl; + bool resTooLarge = r0 > atol; + while (resTooLarge && (it < maxit)) { + const V dx = solveJacobianSystem(); + + updateState(dx, x, xw); + assemble(dt, x, xw, polymer_inflow); + + const double r = residualNorm(); + + const double rr_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm(); + resTooLarge = (r > atol) && (r > rtol*r0); + + it += 1; + std::cout << std::setw(9) << it << std::setprecision(9) + << std::setw(18) << r << std::setprecision(9) + << std::setw(18) << rr_polymer << std::endl; + } + + if (resTooLarge) { + std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n"; + // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); + } + } + + + + + + FullyImplicitCompressiblePolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant() + : accum(2, ADB::null()) + , mflux( ADB::null()) + , b ( ADB::null()) + , head ( ADB::null()) + , mob ( ADB::null()) + , ads (2, ADB::null()) + { + } + + + + + + FullyImplicitCompressiblePolymerSolver::SolutionState::SolutionState(const int np) + : pressure ( ADB::null()) + , saturation(np, ADB::null()) + , concentration( ADB::null()) + , qs ( ADB::null()) + , bhp ( ADB::null()) + { + } + + + + + + FullyImplicitCompressiblePolymerSolver:: + WellOps::WellOps(const Wells& wells) + : w2p(wells.well_connpos[ wells.number_of_wells ], + wells.number_of_wells) + , p2w(wells.number_of_wells, + wells.well_connpos[ wells.number_of_wells ]) + { + const int nw = wells.number_of_wells; + const int* const wpos = wells.well_connpos; + + typedef Eigen::Triplet Tri; + + std::vector scatter, gather; + scatter.reserve(wpos[nw]); + gather .reserve(wpos[nw]); + + for (int w = 0, i = 0; w < nw; ++w) { + for (; i < wpos[ w + 1 ]; ++i) { + scatter.push_back(Tri(i, w, 1.0)); + gather .push_back(Tri(w, i, 1.0)); + } + } + + w2p.setFromTriplets(scatter.begin(), scatter.end()); + p2w.setFromTriplets(gather .begin(), gather .end()); + } + + + + + + FullyImplicitCompressiblePolymerSolver::SolutionState + FullyImplicitCompressiblePolymerSolver::constantState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw) + { + const int nc = grid_.number_of_cells; + const int np = x.numPhases(); + + // The block pattern assumes the following primary variables: + // pressure + // water saturation (if water present) + // polymer concentration + // 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. + std::vector bpat(np + 1, 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); + V so = V::Ones(nc, 1); + const V sw = s.col(0); + so -= sw; + state.saturation[0] = ADB::constant(sw, bpat); + state.saturation[1] = ADB::constant(so, bpat); + + // Concentration + assert(not x.concentration().empty()); + const V c = Eigen::Map(&x.concentration()[0], nc); + state.concentration = ADB::constant(c); + + // 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); + + return state; + } + + + + + + FullyImplicitCompressiblePolymerSolver::SolutionState + FullyImplicitCompressiblePolymerSolver::variableState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw) + { + const int nc = grid_.number_of_cells; + const int np = x.numPhases(); + + std::vector vars0; + + // Initial pressure. + assert (not x.pressure().empty()); + const V p = Eigen::Map(& x.pressure()[0], nc, 1); + vars0.push_back(p); + + // Initial saturation. + assert (not x.saturation().empty()); + const DataBlock s = Eigen::Map(& x.saturation()[0], nc, np); + const V sw = s.col(0); + vars0.push_back(sw); + + // Initial concentration. + assert (not x.concentration().empty()); + const V c = Eigen::Map(&x.concentration()[0], nc); + vars0.push_back(c); + + // Initial 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); + vars0.push_back(qs); + + // Initial well bottom-hole pressure. + assert (not xw.bhp().empty()); + const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + vars0.push_back(bhp); + + std::vector vars = ADB::variables(vars0); + + SolutionState state(np); + + // Pressure. + int nextvar = 0; + state.pressure = vars[ nextvar++ ]; + + // Saturation. + const std::vector& bpat = vars[0].blockPattern(); + { + ADB so = ADB::constant(V::Ones(nc, 1), bpat); + ADB& sw = vars[ nextvar++ ]; + state.saturation[0] = sw; + so = so - sw; + state.saturation[1] = so; + } + + // Concentration. + state.concentration = vars[nextvar++]; + + // Qs. + state.qs = vars[ nextvar++ ]; + + // Bhp. + state.bhp = vars[ nextvar++ ]; + + assert(nextvar == int(vars.size())); + + return state; + } + + + + + + void + FullyImplicitCompressiblePolymerSolver::computeAccum(const SolutionState& state, + const int aix ) + { + + const ADB& press = state.pressure; + const std::vector& sat = state.saturation; + const ADB& c = state.concentration; + + const std::vector cond = phaseCondition(); + std::vector pressure = computePressures(state); + + const ADB pv_mult = poroMult(press); + + for (int phase = 0; phase < 2; ++phase) { + rq_[phase].b = fluidReciprocFVF(phase, pressure[phase], cond, cells_); + } + rq_[0].accum[aix] = pv_mult * rq_[0].b * sat[0]; + rq_[1].accum[aix] = pv_mult * rq_[1].b * sat[1]; + const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); + const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); + const double rho_rock = polymer_props_ad_.rockDensity(); + const V phi = Eigen::Map(&fluid_.porosity()[0], grid_.number_of_cells, 1); + + const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); + rq_[2].accum[aix] = pv_mult * rq_[0].b * sat[0] * c * (1. - dead_pore_vol) + pv_mult * rho_rock * (1. - phi) / phi * ads; + } + + + + + + void + FullyImplicitCompressiblePolymerSolver:: + computeCmax(PolymerBlackoilState& state, + const ADB& c) + { + const int nc = grid_.number_of_cells; + for (int i = 0; i < nc; ++i) { + cmax_(i) = std::max(cmax_(i), c.value()(i)); + } + std::copy(&cmax_[0], &cmax_[0] + nc, state.maxconcentration().begin()); + + } + + + + + + void + FullyImplicitCompressiblePolymerSolver:: + assemble(const double dt, + const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw, + const std::vector& polymer_inflow) + { + // Create the primary variables. + // + const SolutionState state = variableState(x, xw); + const V pvdt = geo_.poreVolume() / dt; + // -------- Mass balance equations -------- + + // Compute b_p and the accumulation term b_p*s_p for each phase, + // except gas. For gas, we compute b_g*s_g + Rs*b_o*s_o. + // These quantities are stored in rq_[phase].accum[1]. + // The corresponding accumulation terms from the start of + // the timestep (b^0_p*s^0_p etc.) were already computed + // in step() and stored in rq_[phase].accum[0]. + computeAccum(state, 1); + + // Set up the common parts of the mass balance equations + // for each active phase. + const V trans = subset(geo_.transmissibility(), ops_.internal_faces); + const std::vector kr = computeRelPerm(state); + const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); + const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); + const ADB mc = computeMc(state); + computeMassFlux(trans, mc, kr[1], krw_eff, state); + residual_.material_balance_eq[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0]) + + ops_.div*rq_[0].mflux; + residual_.material_balance_eq[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0]) + + ops_.div*rq_[1].mflux; + residual_.material_balance_eq[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0]) //+ cell / dt * (rq_[2].ads[1] - rq_[2].ads[0]) + + ops_.div*rq_[2].mflux; + + // -------- Extra (optional) sg or rs equation, and rs contributions to the mass balance equations -------- + + // Add the extra (flux) terms to the gas mass balance equations + // from gas dissolved in the oil phase. + // The extra terms in the accumulation part of the equation are already handled. + // -------- Well equation, and well contributions to the mass balance equations -------- + + // Contribution to mass balance will have to wait. + + 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 std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + const V transw = Eigen::Map(wells_.WI, nperf); + + const ADB& bhp = state.bhp; + + const DataBlock well_s = wops_.w2p * Eigen::Map(wells_.comp_frac, nw, np).matrix(); + + // Extract variables for perforation cell pressures + // and corresponding perforation well pressures. + const ADB p_perfcell = subset(state.pressure, well_cells); + // Finally construct well perforation pressures and well flows. + + // Compute well pressure differentials. + // Construct pressure difference vector for wells. + const int dim = grid_.dimensions; + const double* g = geo_.gravity(); + if (g) { + // Guard against gravity in anything but last dimension. + for (int dd = 0; dd < dim - 1; ++dd) { + assert(g[dd] == 0.0); + } + } + ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern()); + std::vector press = computePressures(state); + const std::vector cond = phaseCondition(); + for (int phase = 0; phase < 2; ++phase) { + const ADB cell_rho = fluidDensity(phase, press[phase], cond, cells_); + cell_rho_total += state.saturation[phase] * cell_rho; + } + ADB inj_rho_total = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); + assert(np == wells_.number_of_phases); + const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); + for (int phase = 0; phase < 2; ++phase) { + const ADB cell_rho = fluidDensity(phase, press[phase], cond, cells_); + const V fraction = compi.col(phase); + inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells); + } + const V rho_perf_cell = subset(cell_rho_total, well_cells).value(); + const V rho_perf_well = inj_rho_total.value(); + V prodperfs = V::Constant(nperf, -1.0); + for (int w = 0; w < nw; ++w) { + if (wells_.type[w] == PRODUCER) { + std::fill(prodperfs.data() + wells_.well_connpos[w], + prodperfs.data() + wells_.well_connpos[w+1], 1.0); + } + } + const Selector producer(prodperfs); + const V rho_perf = producer.select(rho_perf_cell, rho_perf_well); + const V well_perf_dp = computePerfPress(grid_, wells_, rho_perf, g ? g[dim-1] : 0.0); + + const ADB p_perfwell = wops_.w2p * bhp + well_perf_dp; + const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); + // DUMP(nkgradp_well); + const Selector cell_to_well_selector(nkgradp_well.value()); + ADB well_rates_all = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern()); + ADB perf_total_mob = subset(rq_[0].mob, well_cells) + subset(rq_[1].mob, well_cells); + std::vector well_contribs(np, ADB::null()); + std::vector well_perf_rates(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const ADB& cell_b = rq_[phase].b; + const ADB perf_b = subset(cell_b, well_cells); + const ADB& cell_mob = rq_[phase].mob; + const V well_fraction = compi.col(phase); + // Using total mobilities for all phases for injection. + const ADB perf_mob_injector = (wops_.w2p * well_fraction.matrix()).array() * perf_total_mob; + const ADB perf_mob = producer.select(subset(cell_mob, well_cells), + perf_mob_injector); + const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations. + well_perf_rates[phase] = (perf_flux * perf_b); + const ADB well_rates = wops_.p2w * well_perf_rates[phase]; + well_rates_all += superset(well_rates, Span(nw, 1, phase*nw), nw*np); + + // const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc); + well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc); + // DUMP(well_contribs[phase]); + residual_.material_balance_eq[phase] += well_contribs[phase]; + } + + // well rates contribs to polymer mass balance eqn. + // for injection wells. + const V polyin = Eigen::Map(& polymer_inflow[0], nc); + const V poly_in_perf = subset(polyin, well_cells); + const V poly_mc_cell = subset(mc, well_cells).value(); + const V poly_in_c = poly_in_perf;// * poly_mc_cell; + const V poly_mc = producer.select(poly_mc_cell, poly_in_c); + + residual_.material_balance_eq[2] += superset(well_perf_rates[0] * poly_mc, well_cells, nc); + // Set the well flux equation + 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 FullyImplicitCompressiblePolymerSolver::solveJacobianSystem() const + { + return linsolver_.computeNewtonIncrement(residual_); + } + + + + + + namespace { + struct Chop01 { + double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); } + }; + } + + + + + + void FullyImplicitCompressiblePolymerSolver:: + updateState(const V& dx, + PolymerBlackoilState& state, + WellStateFullyImplicitBlackoil& well_state) const + { + const int np = fluid_.numPhases(); + const int nc = grid_.number_of_cells; + const int nw = wells_.number_of_wells; + const V one = V::Constant(nc, 1.0); + const V zero = V::Zero(nc); + // Extract parts of dx corresponding to each part. + const V dp = subset(dx, Span(nc)); + int varstart = nc; + const V dsw =subset(dx, Span(nc, 1, varstart)); + varstart += dsw.size(); + const V dc = subset(dx, Span(nc, 1, varstart)); + varstart += dc.size(); + const V dqs = subset(dx, Span(np*nw, 1, varstart)); + varstart += dqs.size(); + const V dbhp = subset(dx, Span(nw, 1, varstart)); + varstart += dbhp.size(); + assert(varstart == dx.size()); + + // Pressure update. + const double dpmaxrel = 0.8; + const V p_old = Eigen::Map(&state.pressure()[0], nc, 1); + const V absdpmax = dpmaxrel*p_old.abs(); + const V dp_limited = sign(dp) * dp.abs().min(absdpmax); + const V p = (p_old - dp_limited).max(zero); + std::copy(&p[0], &p[0] + nc, state.pressure().begin()); + + // Saturation updates. + const double dsmax = 0.3; + const DataBlock s_old = Eigen::Map(&state.saturation()[0], nc, np); + V so = one; + const V sw_old = s_old.col(0); + const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax); + const V sw = (sw_old - dsw_limited).unaryExpr(Chop01()); + so -= sw; + for (int c = 0; c < nc; ++c) { + state.saturation()[c*np] = sw[c]; + state.saturation()[c*np + 1] = so[c]; + } + + // Concentration updates. +// const double dcmax = 0.3 * polymer_props_ad_.cMax(); +// std::cout << "\n the max concentration: " << dcmax / 0.3 << std::endl; + const V c_old = Eigen::Map(&state.concentration()[0], nc, 1); +// const V dc_limited = sign(dc) * dc.abs().min(dcmax); +// const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02()); + const V c = (c_old - dc).max(zero); + std::copy(&c[0], &c[0] + nc, state.concentration().begin()); + + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const V dwr = Eigen::Map(wwr.data(), nw*np); + const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const V wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + + // Bhp update. + const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); + const V bhp = bhp_old - dbhp; + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + + } + + + + + + std::vector + FullyImplicitCompressiblePolymerSolver::computeRelPerm(const SolutionState& state) const + { + const int nc = grid_.number_of_cells; + const std::vector& bpat = state.pressure.blockPattern(); + + const ADB null = ADB::constant(V::Zero(nc, 1), bpat); + + const ADB sw = state.saturation[0]; + const ADB so = state.saturation[1]; + const ADB sg = null; + return fluid_.relperm(sw, so, sg, cells_); + } + + + + + + std::vector + FullyImplicitCompressiblePolymerSolver::computeRelPermWells(const SolutionState& state, + const DataBlock& well_s, + const std::vector& well_cells) const + { + const int nw = wells_.number_of_wells; + const int nperf = wells_.well_connpos[nw]; + const std::vector& bpat = state.pressure.blockPattern(); + + const ADB null = ADB::constant(V::Zero(nperf), bpat); + + const ADB sw = state.saturation[0]; + const ADB so = state.saturation[1]; + const ADB sg = null; + return fluid_.relperm(sw, so, sg, well_cells); + } + + + + + + std::vector + FullyImplicitCompressiblePolymerSolver:: + computePressures(const SolutionState& state) const + { + const int nc = grid_.number_of_cells; + const std::vector& bpat = state.pressure.blockPattern(); + + const ADB null = ADB::constant(V::Zero(nc, 1), bpat); + + const ADB sw = state.saturation[0]; + + const ADB so = state.saturation[1]; + + const ADB sg = null; + + // convert the pressure offsets to the capillary pressures + std::vector pressure = fluid_.capPress(sw, so, sg, cells_); + pressure[0] = pressure[0] - pressure[1]; + + // add the total pressure to the capillary pressures + pressure[0] = state.pressure - pressure[0]; + pressure[1] = state.pressure + pressure[1]; + + return pressure; + } + + + + + + void + FullyImplicitCompressiblePolymerSolver::computeMassFlux( + const V& transi, + const ADB& mc, + const ADB& kro, + const ADB& krw_eff, + const SolutionState& state ) + { + const ADB tr_mult = transMult(state.pressure); + const std::vector cond = phaseCondition(); + std::vector press = computePressures(state); + + const ADB mu_w = fluidViscosity(0, press[0], cond, cells_); + ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu_w.value().data()); + rq_[0].mob = tr_mult * krw_eff * inv_wat_eff_vis; + rq_[2].mob = tr_mult * mc * krw_eff * inv_wat_eff_vis; + const ADB mu_o = fluidViscosity(1, press[1], cond, cells_); + rq_[1].mob = tr_mult * kro / mu_o; + for (int phase = 0; phase < 2; ++phase) { + const ADB rho = fluidDensity(phase, press[phase], cond, cells_); + ADB& head = rq_[ phase ].head; + // compute gravity potensial using the face average as in eclipse and MRST + const ADB rhoavg = ops_.caver * rho; + const ADB dp = ops_.ngrad * press[phase] - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix())); + head = transi*dp; + UpwindSelector upwind(grid_, ops_, head.value()); + const ADB& b = rq_[phase].b; + const ADB& mob = rq_[phase].mob; + rq_[phase].mflux = upwind.select(b * mob) * head; + } + rq_[2].b = rq_[0].b; + rq_[2].head = rq_[0].head; + UpwindSelector upwind(grid_, ops_, rq_[2].head.value()); + rq_[2].mflux = upwind.select(rq_[2].b * rq_[2].mob) * rq_[2].head; + } + + + + + + double + FullyImplicitCompressiblePolymerSolver::residualNorm() const + { + double r = 0; + for (std::vector::const_iterator + b = residual_.material_balance_eq.begin(), + e = residual_.material_balance_eq.end(); + b != e; ++b) + { + r = std::max(r, (*b).value().matrix().lpNorm()); + } + r = std::max(r, residual_.well_flux_eq.value().matrix().lpNorm()); + r = std::max(r, residual_.well_eq.value().matrix().lpNorm()); + + return r; + } + + + + + + ADB + FullyImplicitCompressiblePolymerSolver::fluidViscosity(const int phase, + const ADB& p , + const std::vector& cond, + const std::vector& cells) const + { + const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern()); + switch (phase) { + case Water: + return fluid_.muWat(p, cells); + case Oil: { + return fluid_.muOil(p, null, cond, cells); + } + default: + OPM_THROW(std::runtime_error, "Unknown phase index " << phase); + } + } + + + + + + ADB + FullyImplicitCompressiblePolymerSolver::fluidReciprocFVF(const int phase, + const ADB& p , + const std::vector& cond, + const std::vector& cells) const + { + const ADB null = ADB::constant(V::Zero(grid_.number_of_cells, 1), p.blockPattern()); + switch (phase) { + case Water: + return fluid_.bWat(p, cells); + case Oil: { + return fluid_.bOil(p, null, cond, cells); + } + default: + OPM_THROW(std::runtime_error, "Unknown phase index " << phase); + } + } + + + + + + ADB + FullyImplicitCompressiblePolymerSolver::fluidDensity(const int phase, + const ADB& p , + const std::vector& cond, + const std::vector& cells) const + { + const double* rhos = fluid_.surfaceDensity(); + ADB b = fluidReciprocFVF(phase, p, cond, cells); + ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b; + return rho; + } + + + + + + // here mc means m(c) * c. + ADB + FullyImplicitCompressiblePolymerSolver::computeMc(const SolutionState& state) const + { + ADB c = state.concentration; + return polymer_props_ad_.polymerWaterVelocityRatio(c); + } + + + + + + ADB + FullyImplicitCompressiblePolymerSolver::poroMult(const ADB& p) const + { + const int n = p.size(); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + V pm(n); + V dpm(n); + for (int i = 0; i < n; ++i) { + pm[i] = rock_comp_props_->poroMult(p.value()[i]); + dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]); + } + ADB::M dpm_diag = spdiag(dpm); + const int num_blocks = p.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dpm_diag * p.derivative()[block]; + } + return ADB::function(pm, jacs); + } else { + return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); + } + } + + + + + + ADB + FullyImplicitCompressiblePolymerSolver::transMult(const ADB& p) const + { + const int n = p.size(); + if (rock_comp_props_ && rock_comp_props_->isActive()) { + V tm(n); + V dtm(n); + for (int i = 0; i < n; ++i) { + tm[i] = rock_comp_props_->transMult(p.value()[i]); + dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]); + } + ADB::M dtm_diag = spdiag(dtm); + const int num_blocks = p.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dtm_diag * p.derivative()[block]; + } + return ADB::function(tm, jacs); + } else { + return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); + } + } + + + void + FullyImplicitCompressiblePolymerSolver::classifyCondition(const PolymerBlackoilState& state) + { + const int nc = grid_.number_of_cells; + const DataBlock s = Eigen::Map(& state.saturation()[0], nc, 2); + + const V so = s.col(1); + for (V::Index c = 0, e = so.size(); c != e; ++c) { + phaseCondition_[c].setFreeWater(); + if (so[c] > 0) { + phaseCondition_[c].setFreeOil(); + } + } + } + +} //namespace Opm diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp new file mode 100644 index 000000000..282241511 --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -0,0 +1,246 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED +#define OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include + +struct UnstructuredGrid; +struct Wells; + +namespace Opm { + + class DerivedGeology; + class RockCompressibility; + class NewtonIterationBlackoilInterface; + class PolymerBlackoilState; + class WellStateFullyImplicitBlackoil; + + /// A fully implicit solver for the oil-water with polymer problem. + /// + /// The simulator is capable of handling oil-water-polymer problems + /// It uses an industry-standard TPFA discretization with per-phase + /// upwind weighting of mobilities. + /// + /// It uses automatic differentiation via the class AutoDiffBlock + /// to simplify assembly of the jacobian matrix. + class FullyImplicitCompressiblePolymerSolver + { + public: + /// Construct a solver. It will retain references to the + /// arguments of this functions, and they are expected to + /// remain in scope for the lifetime of the solver. + /// \param[in] grid grid data structure + /// \param[in] fluid fluid properties + /// \param[in] geo rock properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] polymer_props_ad polymer properties + /// \param[in] wells well structure + /// \param[in] linsolver linear solver + FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid& grid , + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const PolymerPropsAd& polymer_props_ad, + const Wells& wells, + const NewtonIterationBlackoilInterface& linsolver); + + /// Take a single forward step, modifiying + /// state.pressure() + /// state.faceflux() + /// state.saturation() + /// state.concentration() + /// wstate.bhp() + /// \param[in] dt time step size + /// \param[in] state reservoir state + /// \param[in] wstate well state + /// \param[in] polymer_inflow polymer influx + void + step(const double dt, + PolymerBlackoilState& state , + WellStateFullyImplicitBlackoil& wstate, + const std::vector& polymer_inflow); + + private: + typedef AutoDiffBlock ADB; + typedef ADB::V V; + typedef ADB::M M; + typedef Eigen::Array DataBlock; + + struct ReservoirResidualQuant { + ReservoirResidualQuant(); + std::vector accum; // Accumulations + ADB mflux; // Mass flux (surface conditions) + ADB b; // Reciprocal FVF + ADB head; // Pressure drop across int. interfaces + ADB mob; // Phase mobility (per cell) + std::vector ads; // Adsorption term. + }; + + struct SolutionState { + SolutionState(const int np); + ADB pressure; + std::vector saturation; + ADB concentration; + ADB qs; + ADB bhp; + }; + + struct WellOps { + WellOps(const Wells& wells); + M w2p; // well -> perf (scatter) + M p2w; // perf -> well (gather) + }; + + enum { Water = BlackoilPropsAdInterface::Water, + Oil = BlackoilPropsAdInterface::Oil }; + + // Member data + const UnstructuredGrid& grid_; + const BlackoilPropsAdInterface& fluid_; + const DerivedGeology& geo_; + const RockCompressibility* rock_comp_props_; + const PolymerPropsAd& polymer_props_ad_; + const Wells& wells_; + const NewtonIterationBlackoilInterface& linsolver_; + const std::vector cells_; // All grid cells + HelperOps ops_; + const WellOps wops_; + const M grav_; + V cmax_; + std::vector phaseCondition_; + std::vector rq_; + // The mass_balance vector has one element for each active phase, + // each of which has size equal to the number of cells. + // The well_eq has size equal to the number of wells. + LinearisedBlackoilResidual residual_; + // Private methods. + SolutionState + constantState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw); + + SolutionState + variableState(const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw); + + void + computeAccum(const SolutionState& state, + const int aix ); + + void + assemble(const double dt, + const PolymerBlackoilState& x, + const WellStateFullyImplicitBlackoil& xw, + const std::vector& polymer_inflow); + + V solveJacobianSystem() const; + + void updateState(const V& dx, + PolymerBlackoilState& state, + WellStateFullyImplicitBlackoil& well_state) const; + + std::vector + computeRelPerm(const SolutionState& state) const; + + std::vector + computeRelPermWells(const SolutionState& state, + const DataBlock& well_s, + const std::vector& well_cells) const; + std::vector + computePressures(const SolutionState& state) const; + + void + computeMassFlux(const int actph , + const V& transi, + const std::vector& kr , + const SolutionState& state ); + + void + computeMassFlux(const V& trans, + const ADB& mc, + const ADB& kro, + const ADB& krw_eff, + const SolutionState& state); + + std::vector + computeFracFlow(const ADB& kro, + const ADB& krw_eff, + const ADB& c) const; + + void + computeCmax(PolymerBlackoilState& state, + const ADB& c); + + ADB + computeMc(const SolutionState& state) const; + + ADB + rockPorosity(const ADB& p) const; + + ADB + rockPermeability(const ADB& p) const; + + double + residualNorm() const; + + ADB + fluidViscosity(const int phase, + const ADB& p , + const std::vector& cond, + const std::vector& cells) const; + + ADB + fluidReciprocFVF(const int phase, + const ADB& p , + const std::vector& cond, + const std::vector& cells) const; + + ADB + fluidDensity(const int phase, + const ADB& p , + const std::vector& cond, + const std::vector& cells) const; + + ADB + poroMult(const ADB& p) const; + + ADB + transMult(const ADB& p) const; + + const std::vector + phaseCondition() const { return phaseCondition_; } + + void + classifyCondition(const PolymerBlackoilState& state); + }; +} // namespace Opm + + +#endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp new file mode 100644 index 000000000..6a83b826d --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -0,0 +1,896 @@ +/* + Copyright 2014 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 + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace Opm { + + + + + +typedef AutoDiffBlock ADB; +typedef ADB::V V; +typedef ADB::M M; +typedef Eigen::Array DataBlock; + + + + + +namespace { + + std::vector + buildAllCells(const int nc) + { + std::vector all_cells(nc); + for (int c = 0; c < nc; ++c) { all_cells[c] = c; } + + return all_cells; + } + struct Chop01 { + double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); } + }; + + + + + + V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav) + { + const int nw = wells.number_of_wells; + const int nperf = wells.well_connpos[nw]; + const int dim = grid.dimensions; + V wdp = V::Zero(nperf,1); + assert(wdp.size() == rho.size()); + + // Main loop, iterate over all perforations, + // using the following formula: + // wdp(perf) = g*(perf_z - well_ref_z)*rho(perf) + // where the total density rho(perf) is taken to be + // sum_p (rho_p*saturation_p) in the perforation cell. + // [although this is computed on the outside of this function]. + for (int w = 0; w < nw; ++w) { + const double ref_depth = wells.depth_ref[w]; + for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) { + const int cell = wells.well_cells[j]; + const double cell_depth = grid.cell_centroids[dim * cell + dim - 1]; + wdp[j] = rho[j]*grav*(cell_depth - ref_depth); + } + } + return wdp; + } + +} //anonymous namespace + + + + + + FullyImplicitTwophasePolymerSolver:: + FullyImplicitTwophasePolymerSolver(const UnstructuredGrid& grid, + const IncompPropsAdInterface& fluid, + const PolymerPropsAd& polymer_props_ad, + const LinearSolverInterface& linsolver, + const Wells& wells, + const double* gravity) + : grid_ (grid) + , fluid_(fluid) + , polymer_props_ad_ (polymer_props_ad) + , linsolver_(linsolver) + , wells_(wells) + , gravity_(gravity) + , cells_ (buildAllCells(grid.number_of_cells)) + , ops_(grid) + , wops_(wells) +// , mob_(std::vector(fluid.numPhases() + 1, ADB::null())) + , cmax_(V::Zero(grid.number_of_cells)) + , rq_(fluid.numPhases() + 1) + , residual_( { std::vector(fluid.numPhases() + 1, ADB::null()), ADB::null(), ADB::null()}) + { + } + + + + + + FullyImplicitTwophasePolymerSolver:: + WellOps::WellOps(const Wells& wells) + : w2p(wells.well_connpos[ wells.number_of_wells ], + wells.number_of_wells) + , p2w(wells.number_of_wells, + wells.well_connpos[ wells.number_of_wells ]) + { + const int nw = wells.number_of_wells; + const int* const wpos = wells.well_connpos; + + typedef Eigen::Triplet Tri; + + std::vector scatter, gather; + scatter.reserve(wpos[nw]); + gather .reserve(wpos[nw]); + + for (int w = 0, i = 0; w < nw; ++w) { + for (; i < wpos[ w + 1 ]; ++i) { + scatter.push_back(Tri(i, w, 1.0)); + gather .push_back(Tri(w, i, 1.0)); + } + } + + w2p.setFromTriplets(scatter.begin(), scatter.end()); + p2w.setFromTriplets(gather .begin(), gather .end()); + } + + + + + + void + FullyImplicitTwophasePolymerSolver:: + step(const double dt, + PolymerState& x, + WellState& xw, + const std::vector& polymer_inflow, + std::vector& src) + { + + V pvol(grid_.number_of_cells); + // Pore volume + const V::Index nc = grid_.number_of_cells; + V rho = V::Constant(pvol.size(), 1, *fluid_.porosity()); + std::transform(grid_.cell_volumes, grid_.cell_volumes + nc, + rho.data(), pvol.data(), + std::multiplies()); + + const V pvdt = pvol / dt; + + const SolutionState old_state = constantState(x, xw); + computeCmax(x, old_state.concentration); + computeAccum(old_state, 0); + const double atol = 1.0e-12; + const double rtol = 5.0e-8; + const int maxit = 40; + assemble(pvdt, x, xw, polymer_inflow, src); + + const double r0 = residualNorm(); + int it = 0; + std::cout << "\nIteration Residual\n" + << std::setw(9) << it << std::setprecision(9) + << std::setw(18) << r0 << std::endl; + bool resTooLarge = r0 > atol; + while (resTooLarge && (it < maxit)) { + const V dx = solveJacobianSystem(); + updateState(dx, x, xw); + + assemble(pvdt, x, xw, polymer_inflow, src); + + const double r = residualNorm(); + + resTooLarge = (r > atol) && (r > rtol*r0); + + it += 1; + std::cout << std::setw(9) << it << std::setprecision(9) + << std::setw(18) << r << std::endl; + } + + if (resTooLarge) { + std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n"; + // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); + } + } + + + + + + FullyImplicitTwophasePolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant() + : accum(2, ADB::null()) + , mflux( ADB::null()) + , b ( ADB::null()) + , head ( ADB::null()) + , mob ( ADB::null()) + { + } + + + + + + FullyImplicitTwophasePolymerSolver::SolutionState::SolutionState(const int np) + : pressure ( ADB::null()) + , saturation (np, ADB::null()) + , concentration ( ADB::null()) + , qs ( ADB::null()) + , bhp ( ADB::null()) + { + } + + + + + + FullyImplicitTwophasePolymerSolver::SolutionState + FullyImplicitTwophasePolymerSolver::constantState(const PolymerState& x, + const WellState& xw) + { + const int nc = grid_.number_of_cells; + const int np = x.numPhases(); + // The block pattern assumes the following primary variables: + // pressure + // water saturation + // polymer concentration + // well surface rates + // well bottom-hole pressure + // Note that oil is assumed to always be present, but is never + // a primary variable. + std::vector bpat(np + 1, 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); + state.pressure = ADB::constant(p); + + // Saturation. + assert (not x.saturation().empty()); + const DataBlock s_all = Eigen::Map(& x.saturation()[0], nc, np); + for (int phase = 0; phase < np; ++phase) { + state.saturation[phase] = ADB::constant(s_all.col(phase)); + } + + // Concentration + assert(not x.concentration().empty()); + const V c = Eigen::Map(&x.concentration()[0], nc); + state.concentration = ADB::constant(c); + + // 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); + + // 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); + return state; + } + + + + + + FullyImplicitTwophasePolymerSolver::SolutionState + FullyImplicitTwophasePolymerSolver::variableState(const PolymerState& x, + const WellState& xw) + { + const int nc = grid_.number_of_cells; + const int np = x.numPhases(); + + std::vector vars0; + vars0.reserve(np + 3); + + // Initial pressure. + assert (not x.pressure().empty()); + const V p = Eigen::Map(& x.pressure()[0], nc); + vars0.push_back(p); + + // Initial saturation. + assert (not x.saturation().empty()); + const DataBlock s_all = Eigen::Map(& x.saturation()[0], nc, np); + const V sw = s_all.col(0); + vars0.push_back(sw); + + // Initial concentration. + assert (not x.concentration().empty()); + const V c = Eigen::Map(&x.concentration()[0], nc); + vars0.push_back(c); + + // Initial 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); + vars0.push_back(qs); + + // Initial well bottom hole pressure. + assert (not xw.bhp().empty()); + const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + vars0.push_back(bhp); + + std::vector vars = ADB::variables(vars0); + + SolutionState state(np); + + // Pressure. + int nextvar = 0; + state.pressure = vars[ nextvar++ ]; + + // Saturation. + const std::vector& bpat = vars[0].blockPattern(); + { + ADB so = ADB::constant(V::Ones(nc, 1), bpat); + ADB sw = vars[ nextvar++ ]; + state.saturation[0] = sw; + so = so - sw; + state.saturation[1] = so; + } + + // Concentration. + state.concentration = vars[nextvar++]; + + // Qs. + state.qs = vars[ nextvar++ ]; + + // BHP. + state.bhp = vars[ nextvar++ ]; + assert(nextvar == int(vars.size())); + + return state; + } + + + + + + void + FullyImplicitTwophasePolymerSolver:: + computeCmax(PolymerState& state, + const ADB& c) + { + const int nc = grid_.number_of_cells; + for (int i = 0; i < nc; ++i) { + cmax_(i) = std::max(cmax_(i), c.value()(i)); + } + std::copy(&cmax_[0], &cmax_[0] + nc, state.maxconcentration().begin()); + + } + + + + + + void + FullyImplicitTwophasePolymerSolver:: + computeAccum(const SolutionState& state, + const int aix ) + { + + const std::vector& sat = state.saturation; + const ADB& c = state.concentration; + + rq_[0].accum[aix] = sat[0]; + rq_[1].accum[aix] = sat[1]; + + const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); + const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); + const double rho_rock = polymer_props_ad_.rockDensity(); + const V phi = Eigen::Map(&fluid_.porosity()[0], grid_.number_of_cells, 1); + const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); + + rq_[2].accum[aix] = sat[0] * c * (1. - dead_pore_vol) + rho_rock * (1. - phi) / phi * ads; + } + + + + + + void + FullyImplicitTwophasePolymerSolver:: + assemble(const V& pvdt, + const PolymerState& x, + const WellState& xw, + const std::vector& polymer_inflow, + std::vector& src) + { + // Create the primary variables. + const SolutionState state = variableState(x, xw); + computeAccum(state, 1); + // -------- Mass balance equations for water and oil -------- + const V trans = subset(transmissibility(), ops_.internal_faces); + const std::vector kr = computeRelPerm(state); + const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); + const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); + const ADB mc = computeMc(state); + computeMassFlux(trans, mc, kr[1], krw_eff, state); + residual_.mass_balance[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0]) + + ops_.div*rq_[0].mflux; + residual_.mass_balance[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0]) + + ops_.div*rq_[1].mflux; + residual_.mass_balance[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0]) + + ops_.div*rq_[2].mflux; + + // -------- Well equation, and well contributions to the mass balance equations -------- + + // Contribution to mass balance will have to wait. + + 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 std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + const V transw = Eigen::Map(wells_.WI, nperf); + + const ADB& bhp = state.bhp; + + const DataBlock well_s = wops_.w2p * Eigen::Map(wells_.comp_frac, nw, np).matrix(); + + // Extract variables for perforation cell pressures + // and corresponding perforation well pressures. + const ADB p_perfcell = subset(state.pressure, well_cells); + // Finally construct well perforation pressures and well flows. + + // Compute well pressure differentials. + // Construct pressure difference vector for wells. + const int dim = grid_.dimensions; + if (gravity_) { + for (int dd = 0; dd < dim -1; ++dd) { + assert(g[dd] == 0.0); + } + } + ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern()); + for (int phase = 0; phase < 2; ++phase) { + // For incompressible flow cell rho is the same. + const ADB cell_rho = fluidDensity(phase, state.pressure); + cell_rho_total += state.saturation[phase] * cell_rho; + } + ADB inj_rho_total = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); + assert(np == wells_.number_of_phases); + const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); + for (int phase = 0; phase < 2; ++phase) { + const ADB cell_rho = fluidDensity(phase, state.pressure); + const V fraction = compi.col(phase); + inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells); + } + const V rho_perf_cell = subset(cell_rho_total, well_cells).value(); + const V rho_perf_well = inj_rho_total.value(); + V prodperfs = V::Constant(nperf, -1.0); + for (int w = 0; w < nw; ++w) { + if (wells_.type[w] == PRODUCER) { + std::fill(prodperfs.data() + wells_.well_connpos[w], + prodperfs.data() + wells_.well_connpos[w+1], 1.0); + } + } + const Selector producer(prodperfs); + const V rho_perf = producer.select(rho_perf_cell, rho_perf_well); + const V well_perf_dp = computePerfPress(grid_, wells_, rho_perf, gravity_ ? gravity_[dim - 1] : 0.0); + + const ADB p_perfwell = wops_.w2p * bhp + well_perf_dp; + const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); + // DUMP(nkgradp_well); + const Selector cell_to_well_selector(nkgradp_well.value()); + ADB well_rates_all = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern()); + + ADB perf_total_mob = subset(rq_[0].mob, well_cells) + subset(rq_[1].mob, well_cells); + + std::vector well_contribs(np, ADB::null()); + std::vector well_perf_rates(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { +// const ADB& cell_b = rq_[phase].b; + // const ADB perf_b = subset(cell_b, well_cells); + const ADB& cell_mob = rq_[phase].mob; + const V well_fraction = compi.col(phase); + // Using total mobilities for all phases for injection. + const ADB perf_mob_injector = (wops_.w2p * well_fraction.matrix()).array() * perf_total_mob; + const ADB perf_mob = producer.select(subset(cell_mob, well_cells), + perf_mob_injector); + const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations. + well_perf_rates[phase] = perf_flux; + const ADB well_rates = wops_.p2w * well_perf_rates[phase]; + well_rates_all += superset(well_rates, Span(nw, 1, phase*nw), nw*np); + + // const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc); + well_contribs[phase] = superset(perf_flux, well_cells, nc); + // DUMP(well_contribs[phase]); + residual_.mass_balance[phase] += well_contribs[phase]; + for (int i = 0; i < nc; ++i) { + src[i] += well_contribs[phase].value()[i]; + } + } + + // well rates contribs to polymer mass balance eqn. + // for injection wells. + const V polyin = Eigen::Map(& polymer_inflow[0], nc); + const V poly_in_perf = subset(polyin, well_cells); + const V poly_mc_cell = subset(mc, well_cells).value(); + const V poly_c = producer.select(poly_mc_cell, poly_in_perf); + residual_.mass_balance[2] += superset(well_perf_rates[0] * poly_c, well_cells, nc); + + // Set the well flux equation + 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 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); + // residual_.well_eq = bhp_residual; + + } + + + + + + std::vector + FullyImplicitTwophasePolymerSolver:: + computePressures(const SolutionState& state) const + { + const ADB sw = state.saturation[0]; + const ADB so = state.saturation[1]; + + // convert the pressure offsets to the capillary pressures + std::vector pressure = fluid_.capPress(sw, so, cells_); + pressure[0] = pressure[0] - pressure[1]; + + // add the total pressure to the capillary pressures + for (int phaseIdx = 0; phaseIdx < 2; ++phaseIdx) { + pressure[phaseIdx] += state.pressure; + } + + return pressure; + } + + + + + + void + FullyImplicitTwophasePolymerSolver::computeMassFlux(const V& trans, + const ADB& mc, + const ADB& kro, + const ADB& krw_eff, + const SolutionState& state ) + { + const double* mus = fluid_.viscosity(); + ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus); + rq_[0].mob = krw_eff * inv_wat_eff_vis; + rq_[1].mob = kro / V::Constant(kro.size(), 1, mus[1]); + rq_[2].mob = mc * krw_eff * inv_wat_eff_vis; + + const int nc = grid_.number_of_cells; + V z(nc); + // Compute z coordinates + for (int c = 0; c < nc; ++c){ + z[c] = grid_.cell_centroids[c * 3 + 2]; + } + std::vector press = computePressures(state); + for (int phase = 0; phase < 2; ++phase) { + const ADB rho = fluidDensity(phase, state.pressure); + ADB& head = rq_[phase].head; + const ADB rhoavg = ops_.caver * rho; + const ADB dp = ops_.ngrad * press[phase] + - gravity_[2] * (rhoavg * (ops_.ngrad * z.matrix())); + head = trans * dp; + UpwindSelector upwind(grid_, ops_, head.value()); + const ADB& mob = rq_[phase].mob; + rq_[phase].mflux = upwind.select(mob) * head; + } + rq_[2].head = rq_[0].head; + UpwindSelector upwind(grid_, ops_, rq_[2].head.value()); + rq_[2].mflux = upwind.select(rq_[2].mob) * rq_[2].head; + } + + + + + + std::vector + FullyImplicitTwophasePolymerSolver::accumSource(const ADB& kro, + const ADB& krw_eff, + const ADB& c, + const std::vector& src, + const std::vector& polymer_inflow_c) const + { + //extract the source to out and in source. + std::vector outsrc; + std::vector insrc; + std::vector::const_iterator it; + for (it = src.begin(); it != src.end(); ++it) { + if (*it < 0) { + outsrc.push_back(*it); + insrc.push_back(0.0); + } else if (*it > 0) { + insrc.push_back(*it); + outsrc.push_back(0.0); + } else { + outsrc.push_back(0); + insrc.push_back(0); + } + } + const V outSrc = Eigen::Map(& outsrc[0], grid_.number_of_cells); + const V inSrc = Eigen::Map(& insrc[0], grid_.number_of_cells); + const V polyin = Eigen::Map(& polymer_inflow_c[0], grid_.number_of_cells); + // compute the out-fracflow. + const std::vector f = computeFracFlow(); + // compute the in-fracflow. + V zero = V::Zero(grid_.number_of_cells); + V one = V::Ones(grid_.number_of_cells); + + std::vector source; + //water source + source.push_back(f[0] * outSrc + one * inSrc); + //oil source + source.push_back(f[1] * outSrc + zero * inSrc); + //polymer source + source.push_back(f[0] * outSrc * c + one * inSrc * polyin); + + return source; + } + + + + + + std::vector + FullyImplicitTwophasePolymerSolver::computeFracFlow() const + { + ADB total_mob = rq_[0].mob + rq_[1].mob; + + std::vector fracflow; + + fracflow.push_back(rq_[0].mob / total_mob); + fracflow.push_back(rq_[1].mob / total_mob); + + return fracflow; + } + + + + + + V + FullyImplicitTwophasePolymerSolver::solveJacobianSystem() const + { + const int np = fluid_.numPhases(); + if (np != 2) { + OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwophasePolymerSolver."); + } + ADB mass_res = vertcat(residual_.mass_balance[0], residual_.mass_balance[1]); + mass_res = vertcat(mass_res, residual_.mass_balance[2]); + ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq); + ADB total_res = collapseJacs(vertcat(mass_res, well_res)); + + const Eigen::SparseMatrix matr = total_res.derivative()[0]; + V dx(V::Zero(total_res.size())); + Opm::LinearSolverInterface::LinearSolverReport rep + = linsolver_.solve(matr.rows(), matr.nonZeros(), + matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), + total_res.value().data(), dx.data()); + if (!rep.converged) { + OPM_THROW(std::runtime_error, + "FullyImplicitBlackoilSolver::solveJacobianSystem(): " + "Linear solver convergence failure."); + } + return dx; + } + + + + + + void FullyImplicitTwophasePolymerSolver::updateState(const V& dx, + PolymerState& state, + WellState& well_state) const + { + const int np = fluid_.numPhases(); + const int nc = grid_.number_of_cells; + const int nw = wells_.number_of_wells; + const V one = V::Constant(nc, 1.0); + const V zero = V::Zero(nc); + + // Extract parts of dx corresponding to each part. + const V dp = subset(dx, Span(nc)); + int varstart = nc; + const V dsw = subset(dx, Span(nc, 1, varstart)); + varstart += dsw.size(); + const V dc = subset(dx, Span(nc, 1, varstart)); + varstart += dc.size(); + const V dqs = subset(dx, Span(np*nw, 1, varstart)); + varstart += dqs.size(); + const V dbhp = subset(dx, Span(nw, 1, varstart)); + varstart += dbhp.size(); + + assert(varstart == dx.size()); + + // Pressure update. + const V p_old = Eigen::Map(&state.pressure()[0], nc); + const V p = p_old - dp; + std::copy(&p[0], &p[0] + nc, state.pressure().begin()); + + // Saturation updates. + const double dsmax = 0.3; + const DataBlock s_old = Eigen::Map(& state.saturation()[0], nc, np); + V so = one; + const V sw_old = s_old.col(0); + const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax); + const V sw = (sw_old - dsw_limited).unaryExpr(Chop01()); + so -= sw; + for (int c = 0; c < nc; ++c) { + state.saturation()[c*np] = sw[c]; + state.saturation()[c*np + 1] = so[c]; + } + + // Concentration updates. + const V c_old = Eigen::Map(&state.concentration()[0], nc); + const V c = (c_old - dc).max(zero); + std::copy(&c[0], &c[0] + nc, state.concentration().begin()); + + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const V dwr = Eigen::Map(wwr.data(), nw*np); + const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const V wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + + // Bhp update. + const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); + const V bhp = bhp_old - dbhp; + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + } + + + + + + std::vector + FullyImplicitTwophasePolymerSolver::computeRelPerm(const SolutionState& state) const + { + + const ADB sw = state.saturation[0]; + const ADB so = state.saturation[1]; + + return fluid_.relperm(sw, so, cells_); + } + + + + + + double + FullyImplicitTwophasePolymerSolver::residualNorm() const + { + double r = 0; + for (std::vector::const_iterator + b = residual_.mass_balance.begin(), + e = residual_.mass_balance.end(); + b != e; ++b) + { + r = std::max(r, (*b).value().matrix().lpNorm()); + } + + r = std::max(r, residual_.well_flux_eq.value().matrix().lpNorm()); + r = std::max(r, residual_.well_eq.value().matrix().lpNorm()); + + return r; + } + + + + + + ADB + FullyImplicitTwophasePolymerSolver::fluidDensity(const int phase, + const ADB p) const + { + const double* rhos = fluid_.surfaceDensity(); + ADB rho = ADB::constant(V::Constant(grid_.number_of_cells, 1, rhos[phase]), + p.blockPattern()); + + return rho; + } + + + + + + V + FullyImplicitTwophasePolymerSolver::transmissibility() const + { + const V::Index nc = grid_.number_of_cells; + V htrans(grid_.cell_facepos[nc]); + V trans(grid_.cell_facepos[nc]); + UnstructuredGrid* ug = const_cast(& grid_); + tpfa_htrans_compute(ug, fluid_.permeability(), htrans.data()); + tpfa_trans_compute (ug, htrans.data(), trans.data()); + + return trans; + } + + + + + + // here mc means m(c) * c. + ADB + FullyImplicitTwophasePolymerSolver::computeMc(const SolutionState& state) const + { + ADB c = state.concentration; + return polymer_props_ad_.polymerWaterVelocityRatio(c); + } + + + + + +} //namespace Opm diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp new file mode 100644 index 000000000..cb7be065f --- /dev/null +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -0,0 +1,217 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED +#define OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include + +struct UnstructuredGrid; +struct Wells; + +namespace Opm { + + class LinearSolverInterface; + class PolymerState; + class WellState; + + /// A fully implicit solver for the incompressible oil-water flow wtih polymer problem. + /// + /// The simulator is capable of handling oil-water with polymer problems + /// It uses an industry-standard TPFA discretization with per-phase + /// upwind weighting of mobilities. + /// + /// It uses automatic differentiation via the class AutoDiffBlock + /// to simplify assembly of the jacobian matrix. + class FullyImplicitTwophasePolymerSolver + { + public: + /// Construct a solver. It will retain references to the + /// arguments of this functions, and they are expected to + /// remain in scope for the lifetime of the solver. + /// \param[in] grid grid data structure + /// \param[in] fluid fluid properties + /// \param[in] polymer_props_ad polymer properties + /// \param[in] wells well structure + /// \param[in] linsolver linear solver + /// \param[in] gravity gravity + FullyImplicitTwophasePolymerSolver(const UnstructuredGrid& grid, + const IncompPropsAdInterface& fluid, + const PolymerPropsAd& polymer_props_ad, + const LinearSolverInterface& linsolver, + const Wells& wells, + const double* gravity); + + /// Take a single forward step, modifiying + /// state.pressure() + /// state.faceflux() + /// state.saturation() + /// state.concentration() + /// wstate.bhp() + /// \param[in] dt time step size + /// \param[in] state reservoir state with polymer + /// \param[in] wstate well state + /// \param[in] polymer_inflow polymer influx + /// \param[in] src to caculate wc + void step(const double dt, + PolymerState& state, + WellState& well_state, + const std::vector& polymer_inflow, + std::vector& src); + private: + typedef AutoDiffBlock ADB; + typedef ADB::V V; + typedef ADB::M M; + typedef Eigen::Array DataBlock; + + struct ReservoirResidualQuant { + ReservoirResidualQuant(); + std::vector accum; // Accumulations + ADB mflux; // Mass flux (surface conditions) + ADB b; // Reciprocal FVF + ADB head; // Pressure drop across int. interfaces + ADB mob; // Phase mobility (per cell) + }; + + struct SolutionState { + SolutionState(const int np); + ADB pressure; + std::vector saturation; + ADB concentration; + ADB qs; + ADB bhp; + }; + + struct WellOps { + WellOps(const Wells& wells); + M w2p; // well -> perf (scatter) + M p2w; // perf -> well (gather) + }; + + const UnstructuredGrid& grid_; + const IncompPropsAdInterface& fluid_; + const PolymerPropsAd& polymer_props_ad_; + const LinearSolverInterface& linsolver_; + const Wells& wells_; + const double* gravity_; + const std::vector cells_; + HelperOps ops_; + const WellOps wops_; + V cmax_; + std::vector rq_; + + struct { + std::vector mass_balance; + ADB well_eq; + ADB well_flux_eq; + } residual_; + + SolutionState + constantState(const PolymerState& x, + const WellState& xw); + SolutionState + variableState(const PolymerState& x, + const WellState& xw); + void + assemble(const V& pvdt, + const PolymerState& x, + const WellState& xw, + const std::vector& polymer_inflow, + std::vector& src); + + V solveJacobianSystem() const; + + void updateState(const V& dx, + PolymerState& x, + WellState& xw) const; + std::vector + computeRelPerm(const SolutionState& state) const; + + V + transmissibility() const; + + std::vector + computePressures(const SolutionState& state) const; + + void + computeMassFlux(const V& trans, + const ADB& mc, + const ADB& kro, + const ADB& krw_eff, + const SolutionState& state ); + + std::vector + accumSource(const ADB& kro, + const ADB& krw_eff, + const ADB& c, + const std::vector& src, + const std::vector& polymer_inflow_c) const; + + + std::vector + computeFracFlow() const; + + double + residualNorm() const; + + ADB + polymerSource(const std::vector& kr, + const std::vector& src, + const std::vector& polymer_inflow_c, + const SolutionState& state) const; + + void + computeCmax(PolymerState& state, + const ADB& c); + void + computeAccum(const SolutionState& state, + const int aix ); + ADB + computeMc(const SolutionState& state) const; + + ADB + rockPorosity(const ADB& p) const; + + ADB + rockPermeability(const ADB& p) const; + + const double + fluidDensity(const int phase) const; + + ADB + fluidDensity(const int phase, + const ADB p) const; + + ADB + transMult(const ADB& p) const; + }; + + +} // namespace Opm + +#endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp new file mode 100644 index 000000000..897a8f3d9 --- /dev/null +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp @@ -0,0 +1,272 @@ +/* + Copyright 2014 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 +#include +#include +#include +#include + +namespace Opm { + + typedef PolymerPropsAd::ADB ADB; + typedef PolymerPropsAd::V V; + + + + + + double + PolymerPropsAd::rockDensity() const + { + return polymer_props_.rockDensity(); + } + + + + + + double + PolymerPropsAd::deadPoreVol() const + { + return polymer_props_.deadPoreVol(); + } + + + + + + double + PolymerPropsAd::cMax() const + { + return polymer_props_.cMax(); + } + + + + + + PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props) + : polymer_props_ (polymer_props) + { + } + + + + + + PolymerPropsAd::~PolymerPropsAd() + { + } + + + + + + V PolymerPropsAd::effectiveInvWaterVisc(const V& c, + const double* visc) const + { + const int nc = c.size(); + V inv_mu_w_eff(nc); + for (int i = 0; i < nc; ++i) { + double im = 0; + polymer_props_.effectiveInvVisc(c(i), visc, im); + inv_mu_w_eff(i) = im; + } + + return inv_mu_w_eff; + } + + + + + + ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c, + const double* visc) const + { + const int nc = c.size(); + V inv_mu_w_eff(nc); + V dinv_mu_w_eff(nc); + for (int i = 0; i < nc; ++i) { + double im = 0, dim = 0; + polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim); + inv_mu_w_eff(i) = im; + dinv_mu_w_eff(i) = dim; + } + ADB::M dim_diag = spdiag(dinv_mu_w_eff); + const int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dim_diag * c.derivative()[block]; + } + return ADB::function(inv_mu_w_eff, jacs); + } + + + + + + V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const + { + const int nc = c.size(); + V mc(nc); + + for (int i = 0; i < nc; ++i) { + double m = 0; + polymer_props_.computeMc(c(i), m); + mc(i) = m; + } + + return mc; + } + + + + + + ADB PolymerPropsAd::polymerWaterVelocityRatio(const ADB& c) const + { + + const int nc = c.size(); + V mc(nc); + V dmc(nc); + + for (int i = 0; i < nc; ++i) { + double m = 0; + double dm = 0; + polymer_props_.computeMcWithDer(c.value()(i), m, dm); + + mc(i) = m; + dmc(i) = dm; + } + + ADB::M dmc_diag = spdiag(dmc); + const int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dmc_diag * c.derivative()[block]; + } + + return ADB::function(mc, jacs); + } + + + + + + V PolymerPropsAd::adsorption(const V& c, const V& cmax_cells) const + { + const int nc = c.size(); + V ads(nc); + + for (int i = 0; i < nc; ++i) { + double c_ads = 0; + polymer_props_.adsorption(c(i), cmax_cells(i), c_ads); + + ads(i) = c_ads; + } + + return ads; + } + + + + + + ADB PolymerPropsAd::adsorption(const ADB& c, const ADB& cmax_cells) const + { + const int nc = c.value().size(); + + V ads(nc); + V dads(nc); + + for (int i = 0; i < nc; ++i) { + double c_ads = 0; + double dc_ads = 0; + polymer_props_.adsorptionWithDer(c.value()(i), cmax_cells.value()(i), c_ads, dc_ads); + ads(i) = c_ads; + dads(i) = dc_ads; + } + + ADB::M dads_diag = spdiag(dads); + int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dads_diag * c.derivative()[block]; + } + + return ADB::function(ads, jacs); + } + + + + + + V + PolymerPropsAd::effectiveRelPerm(const V& c, + const V& cmax_cells, + const V& krw) const + { + const int nc = c.size(); + + V one = V::Ones(nc); + V ads = adsorption(c, cmax_cells); + double max_ads = polymer_props_.cMaxAds(); + double res_factor = polymer_props_.resFactor(); + double factor = (res_factor -1.) / max_ads; + V rk = one + factor * ads; + V krw_eff = krw / rk; + + return krw_eff; + } + + + + + + ADB + PolymerPropsAd::effectiveRelPerm(const ADB& c, + const ADB& cmax_cells, + const ADB& krw, + const ADB& sw) const + { + const int nc = c.value().size(); + V one = V::Ones(nc); + ADB ads = adsorption(c, cmax_cells); + V krw_eff = effectiveRelPerm(c.value(), cmax_cells.value(), krw.value()); + + double max_ads = polymer_props_.cMaxAds(); + double res_factor = polymer_props_.resFactor(); + double factor = (res_factor - 1.) / max_ads; + ADB rk = one + ads * factor; + ADB dkrw_ds = krw / rk; + ADB dkrw_dc = -factor * krw / (rk * rk) * ads ; + + const int num_blocks = c.numBlocks(); + std::vector jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + jacs[block] = dkrw_ds.derivative()[block] * sw.derivative()[block] + + dkrw_dc.derivative()[block] * c.derivative()[block]; + } + + return ADB::function(krw_eff, jacs); + } + +}// namespace Opm diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp new file mode 100644 index 000000000..4b691144d --- /dev/null +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -0,0 +1,112 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_POLYMERPROPSAD_HEADED_INLCUDED +#define OPM_POLYMERPROPSAD_HEADED_INLCUDED + +#include +#include +#include +#include +#include + +namespace Opm { + + class PolymerPropsAd + { + public: + /// \return Reference rock density. + double rockDensity() const; + + /// \return The value of dead pore volume. + double deadPoreVol() const; + + /// \return The max concentration injected. + double cMax() const; + + typedef AutoDiffBlock ADB; + typedef ADB::V V; + + /// Constructor wrapping a polymer props. + PolymerPropsAd(const PolymerProperties& polymer_props); + + /// Destructor. + ~PolymerPropsAd(); + + /// \param[in] c Array of n polymer concentraion values. + /// \param[in] visc Array of 2 viscosity value. + /// \return value of inverse effective water viscosity. + V + effectiveInvWaterVisc(const V& c,const double* visc) const; + + /// \param[in] c Array of n polymer concentraion values. + /// \param[in] visc Array of 2 viscosity value + /// \return value of inverse effective water viscosity. + ADB + effectiveInvWaterVisc(const ADB& c,const double* visc) const; + + /// \param[in] c Array of n polymer concentraion values. + /// \return Array of n mc values, here mc means m(c) * c. + V + polymerWaterVelocityRatio(const V& c) const; + + /// \param[in] c Array of n polymer concentraion values. + /// \return Array of n mc values, here mc means m(c) * c. + ADB + polymerWaterVelocityRatio(const ADB& c) const; + + /// \param[in] c Array of n polymer concentraion values. + /// \param[in] cmax_cells Array of n polymer concentraion values + /// that the cell experienced. + /// \return Array of n adsorption values. + V + adsorption(const V& c, const V& cmax_cells) const; + + /// \param[in] c Array of n polymer concentraion values. + /// \param[in] cmax_cells Array of n polymer concentraion values + /// that the cell experienced. + /// \return Array of n adsorption values. + ADB + adsorption(const ADB& c, const ADB& cmax_cells) const; + + /// \param[in] c Array of n polymer concentraion values. + /// \param[in] cmax_cells Array of n polymer concentraion values + /// that the cell experienced. + /// \param[in] relperm Array of n relative water relperm values. + /// \return Array of n adsorption values. + V + effectiveRelPerm(const V& c, const V& cmax_cells, const V& relperm) const; + + + /// \param[in] c Array of n polymer concentraion values. + /// \param[in] cmax_cells Array of n polymer concentraion values + /// that the cell experienced. + /// \param[in] relperm Array of n relative water relperm values. + /// \return Array of n adsorption values. + ADB + effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw, const ADB& sw) const; + + private: + const PolymerProperties& polymer_props_; + }; + +} //namespace Opm + +#endif// OPM_POLYMERPROPSAD_HEADED_INLCUDED diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp new file mode 100644 index 000000000..8b165d8a2 --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp @@ -0,0 +1,117 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED +#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED + +#include +#include + +struct UnstructuredGrid; +struct Wells; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class BlackoilPropsAdInterface; + class RockCompressibility; + class DerivedGeology; + class NewtonIterationBlackoilInterface; + class SimulatorTimer; + class PolymerBlackoilState; + class WellStateFullyImplicitBlackoil; + class WellsManager; + class EclipseState; + class EclipseWriter; + class PolymerPropsAd; + class PolymerInflowInterface; + struct SimulatorReport; + + /// Class collecting all necessary components for a two-phase simulation. + template + class SimulatorFullyImplicitBlackoilPolymer + { + public: + /// \brief The type of the grid that we use. + typedef T Grid; + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] geo derived geological properties + /// \param[in] props fluid and rock properties + /// \param[in] polymer_props polymer properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + /// \param[in] disgas true for dissolved gas option + /// \param[in] vapoil true for vaporized oil option + /// \param[in] eclipse_state + /// \param[in] output_writer + /// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow + SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param, + const Grid& grid, + const DerivedGeology& geo, + BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity, + const bool disgas, + const bool vapoil, + const bool polymer, + std::shared_ptr eclipse_state, + EclipseWriter& output_writer, + Opm::DeckConstPtr& deck, + const std::vector& threshold_pressures_by_face); + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + /// \return simulation report, with timing data + SimulatorReport run(SimulatorTimer& timer, + PolymerBlackoilState& state); + + private: + class Impl; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + std::shared_ptr pimpl_; + }; + +} // namespace Opm + +#include "SimulatorFullyImplicitBlackoilPolymer_impl.hpp" +#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp new file mode 100644 index 000000000..cb3d710db --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp @@ -0,0 +1,211 @@ +#include "config.h" + +#include "SimulatorFullyImplicitBlackoilPolymerOutput.hpp" + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#ifdef HAVE_DUNE_CORNERPOINT +#include +#include +#include +#include +#endif +namespace Opm +{ + + + void outputStateVtk(const UnstructuredGrid& grid, + const PolymerBlackoilState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["concentration"] = &state.concentration(); + dm["maxconcentration"] = &state.maxconcentration(); + std::vector cell_velocity; + Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), + AutoDiffGrid::numFaces(grid), + AutoDiffGrid::beginFaceCentroids(grid), + AutoDiffGrid::faceCells(grid), + AutoDiffGrid::beginCellCentroids(grid), + AutoDiffGrid::beginCellVolumes(grid), + AutoDiffGrid::dimensions(grid), + state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + } + + + void outputStateMatlab(const UnstructuredGrid& grid, + const PolymerBlackoilState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["surfvolume"] = &state.surfacevol(); + dm["rs"] = &state.gasoilratio(); + dm["rv"] = &state.rv(); + dm["concentration"] = &state.concentration(); + dm["maxconcentration"] = &state.maxconcentration(); + std::vector cell_velocity; + Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), + AutoDiffGrid::numFaces(grid), + AutoDiffGrid::beginFaceCentroids(grid), + UgGridHelpers::faceCells(grid), + AutoDiffGrid::beginCellCentroids(grid), + AutoDiffGrid::beginCellVolumes(grid), + AutoDiffGrid::dimensions(grid), + state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + void outputWellStateMatlab(const Opm::WellState& well_state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["bhp"] = &well_state.bhp(); + dm["wellrates"] = &well_state.wellRates(); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error,"Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + +#if 0 + void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + watercut.write(os); + } + + void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + wellreport.write(os); + } +#endif + +#ifdef HAVE_DUNE_CORNERPOINT + void outputStateVtk(const Dune::CpGrid& grid, + const PolymerBlackoilState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + std::ostringstream vtkpath; + vtkpath << output_dir << "/vtk_files"; + vtkpath << "/output-" << std::setw(3) << std::setfill('0') << step; + boost::filesystem::path fpath(vtkpath.str()); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step; +#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3) + Dune::VTKWriter writer(grid.leafGridView(), Dune::VTK::nonconforming); +#else + Dune::VTKWriter writer(grid.leafView(), Dune::VTK::nonconforming); +#endif + writer.addCellData(state.saturation(), "saturation", state.numPhases()); + writer.addCellData(state.pressure(), "pressure", 1); + writer.addCellData(state.concentration(), "concentration", 1); + writer.addCellData(state.maxconcentration(), "maxconcentration", 1); + + std::vector cell_velocity; + Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), + AutoDiffGrid::numFaces(grid), + AutoDiffGrid::beginFaceCentroids(grid), + AutoDiffGrid::faceCells(grid), + AutoDiffGrid::beginCellCentroids(grid), + AutoDiffGrid::beginCellVolumes(grid), + AutoDiffGrid::dimensions(grid), + state.faceflux(), cell_velocity); + writer.addCellData(cell_velocity, "velocity", Dune::CpGrid::dimension); + writer.pwrite(vtkfilename.str(), vtkpath.str(), std::string("."), Dune::VTK::ascii); + } +#endif + +} diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp new file mode 100644 index 000000000..55a364324 --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp @@ -0,0 +1,96 @@ +#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED +#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#ifdef HAVE_DUNE_CORNERPOINT +#include +#endif +namespace Opm +{ + + void outputStateVtk(const UnstructuredGrid& grid, + const PolymerBlackoilState& state, + const int step, + const std::string& output_dir); + + + void outputStateMatlab(const UnstructuredGrid& grid, + const PolymerBlackoilState& state, + const int step, + const std::string& output_dir); + + void outputWellStateMatlab(const Opm::WellState& well_state, + const int step, + const std::string& output_dir); +#ifdef HAVE_DUNE_CORNERPOINT + void outputStateVtk(const Dune::CpGrid& grid, + const PolymerBlackoilState& state, + const int step, + const std::string& output_dir); +#endif + + template + void outputStateMatlab(const Grid& grid, + const PolymerBlackoilState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["surfvolume"] = &state.surfacevol(); + dm["rs"] = &state.gasoilratio(); + dm["rv"] = &state.rv(); + dm["concentration"] = &state.concentration(); + dm["maxconcentration"] = &state.maxconcentration(); + + std::vector cell_velocity; + Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), + AutoDiffGrid::numFaces(grid), + AutoDiffGrid::beginFaceCentroids(grid), + UgGridHelpers::faceCells(grid), + AutoDiffGrid::beginCellCentroids(grid), + AutoDiffGrid::beginCellVolumes(grid), + AutoDiffGrid::dimensions(grid), + state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + +} +#endif diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp new file mode 100644 index 000000000..7a745ef41 --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -0,0 +1,621 @@ +/* + Copyright 2014 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 +#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 +#include + +namespace Opm +{ + template + class SimulatorFullyImplicitBlackoilPolymer::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const Grid& grid, + const DerivedGeology& geo, + BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity, + bool has_disgas, + bool has_vapoil, + bool has_polymer, + std::shared_ptr eclipse_state, + EclipseWriter& output_writer, + Opm::DeckConstPtr& deck, + const std::vector& threshold_pressures_by_face); + + SimulatorReport run(SimulatorTimer& timer, + PolymerBlackoilState& state); + + private: + // Data. + typedef RateConverter:: + SurfaceToReservoirVoidage< BlackoilPropsAdInterface, + std::vector > RateConverterType; + + const parameter::ParameterGroup param_; + + // Parameters for output. + bool output_; + bool output_vtk_; + std::string output_dir_; + int output_interval_; + // Observed objects. + const Grid& grid_; + BlackoilPropsAdInterface& props_; + const PolymerPropsAd& polymer_props_; + const RockCompressibility* rock_comp_props_; + const double* gravity_; + // Solvers + const DerivedGeology& geo_; + NewtonIterationBlackoilInterface& solver_; + // Misc. data + std::vector allcells_; + const bool has_disgas_; + const bool has_vapoil_; + const bool has_polymer_; + // eclipse_state + std::shared_ptr eclipse_state_; + // output_writer + EclipseWriter& output_writer_; + Opm::DeckConstPtr& deck_; + RateConverterType rateConverter_; + // Threshold pressures. + std::vector threshold_pressures_by_face_; + + void + computeRESV(const std::size_t step, + const Wells* wells, + const PolymerBlackoilState& x, + WellStateFullyImplicitBlackoil& xw); + }; + + + + + template + SimulatorFullyImplicitBlackoilPolymer::SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param, + const Grid& grid, + const DerivedGeology& geo, + BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity, + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer, + std::shared_ptr eclipse_state, + EclipseWriter& output_writer, + Opm::DeckConstPtr& deck, + const std::vector& threshold_pressures_by_face) + + { + pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, + has_vapoil, has_polymer, eclipse_state, output_writer, deck, threshold_pressures_by_face)); + } + + + + + + template + SimulatorReport SimulatorFullyImplicitBlackoilPolymer::run(SimulatorTimer& timer, + PolymerBlackoilState& state) + { + return pimpl_->run(timer, state); + } + + + + static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["bhp"] = &well_state.bhp(); + dm["wellrates"] = &well_state.wellRates(); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error,"Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + +#if 0 + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + watercut.write(os); + } + + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + wellreport.write(os); + } +#endif + + + // \TODO: Treat bcs. + template + SimulatorFullyImplicitBlackoilPolymer::Impl::Impl(const parameter::ParameterGroup& param, + const Grid& grid, + const DerivedGeology& geo, + BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity, + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer, + std::shared_ptr eclipse_state, + EclipseWriter& output_writer, + Opm::DeckConstPtr& deck, + const std::vector& threshold_pressures_by_face) + : param_(param), + grid_(grid), + props_(props), + polymer_props_(polymer_props), + rock_comp_props_(rock_comp_props), + gravity_(gravity), + geo_(geo), + solver_(linsolver), + has_disgas_(has_disgas), + has_vapoil_(has_vapoil), + has_polymer_(has_polymer), + eclipse_state_(eclipse_state), + output_writer_(output_writer), + deck_(deck), + rateConverter_(props_, std::vector(AutoDiffGrid::numCells(grid_), 0)), + threshold_pressures_by_face_(threshold_pressures_by_face) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + + // Misc init. + const int num_cells = AutoDiffGrid::numCells(grid); + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + template + SimulatorReport SimulatorFullyImplicitBlackoilPolymer::Impl::run(SimulatorTimer& timer, + PolymerBlackoilState& state) + { + WellStateFullyImplicitBlackoil prev_well_state; + + // Create timers and file for writing timing info. + Opm::time::StopWatch solver_timer; + double stime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::string tstep_filename = output_dir_ + "/step_timing.txt"; + std::ofstream tstep_os(tstep_filename.c_str()); + + // Main simulation loop. + while (!timer.done()) { + // Report timestep. + step_timer.start(); + timer.report(std::cout); + + // Create wells and well state. + WellsManager wells_manager(eclipse_state_, + timer.currentStepNum(), + Opm::UgGridHelpers::numCells(grid_), + Opm::UgGridHelpers::globalCell(grid_), + Opm::UgGridHelpers::cartDims(grid_), + Opm::UgGridHelpers::dimensions(grid_), + Opm::UgGridHelpers::beginCellCentroids(grid_), + Opm::UgGridHelpers::cell2Faces(grid_), + Opm::UgGridHelpers::beginFaceCentroids(grid_), + props_.permeability()); + const Wells* wells = wells_manager.c_wells(); + WellStateFullyImplicitBlackoil well_state; + well_state.init(wells, state.blackoilState(), prev_well_state); + // compute polymer inflow + std::unique_ptr polymer_inflow_ptr; + if (deck_->hasKeyword("WPOLYMER")) { + if (wells_manager.c_wells() == 0) { + OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells."); + } + polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, *wells, Opm::UgGridHelpers::numCells(grid_))); + } else { + polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day, + 1.0*Opm::unit::day, + 0.0)); + } + std::vector polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_)); + polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), + timer.simulationTimeElapsed() + timer.currentStepLength(), + polymer_inflow_c); + // Output state at start of time step. + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); + } + if (output_) { + if (timer.currentStepNum() == 0) { + output_writer_.writeInit(timer); + } + output_writer_.writeTimeStep(timer, state.blackoilState(), well_state.basicWellState()); + } + + // Max oil saturation (for VPPARS), hysteresis update. + props_.updateSatOilMax(state.saturation()); + props_.updateSatHyst(state.saturation(), allcells_); + + // Compute reservoir volumes for RESV controls. + computeRESV(timer.currentStepNum(), wells, state, well_state); + + // Run a single step of the solver. + solver_timer.start(); + FullyImplicitBlackoilPolymerSolver solver(param_, grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells, solver_, has_disgas_, has_vapoil_, has_polymer_); + if (!threshold_pressures_by_face_.empty()) { + solver.setThresholdPressures(threshold_pressures_by_face_); + } + solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c); + solver_timer.stop(); + + // Report timing. + const double st = solver_timer.secsSinceStart(); + std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; + stime += st; + if (output_) { + SimulatorReport step_report; + step_report.pressure_time = st; + step_report.total_time = step_timer.secsSinceStart(); + step_report.reportParam(tstep_os); + } + + // Increment timer, remember well state. + ++timer; + prev_well_state = well_state; + } + + // Write final simulation state. + if (output_) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputWellStateMatlab(prev_well_state, timer.currentStepNum(), output_dir_); + output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state.basicWellState()); + } + + // Stop timer and create timing report + total_timer.stop(); + SimulatorReport report; + report.pressure_time = stime; + report.transport_time = 0.0; + report.total_time = total_timer.secsSinceStart(); + return report; + } + + namespace SimFIBODetails { + typedef std::unordered_map WellMap; + + inline WellMap + mapWells(const std::vector& wells) + { + WellMap wmap; + + for (std::vector::const_iterator + w = wells.begin(), e = wells.end(); + w != e; ++w) + { + wmap.insert(std::make_pair((*w)->name(), *w)); + } + + return wmap; + } + + inline int + resv_control(const WellControls* ctrl) + { + int i, n = well_controls_get_num(ctrl); + + bool match = false; + for (i = 0; (! match) && (i < n); ++i) { + match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE; + } + + if (! match) { i = 0; } + + return i - 1; // -1 if no match, undo final "++" otherwise + } + + inline bool + is_resv_prod(const Wells& wells, + const int w) + { + return ((wells.type[w] == PRODUCER) && + (0 <= resv_control(wells.ctrls[w]))); + } + + inline bool + is_resv_prod(const WellMap& wmap, + const std::string& name, + const std::size_t step) + { + bool match = false; + + WellMap::const_iterator i = wmap.find(name); + + if (i != wmap.end()) { + WellConstPtr wp = i->second; + + match = (wp->isProducer(step) && + wp->getProductionProperties(step) + .hasProductionControl(WellProducer::RESV)); + } + + return match; + } + + inline std::vector + resvProducers(const Wells& wells, + const std::size_t step, + const WellMap& wmap) + { + std::vector resv_prod; + + for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) { + if (is_resv_prod(wells, w) || + ((wells.name[w] != 0) && + is_resv_prod(wmap, wells.name[w], step))) + { + resv_prod.push_back(w); + } + } + + return resv_prod; + } + + inline void + historyRates(const PhaseUsage& pu, + const WellProductionProperties& p, + std::vector& rates) + { + assert (! p.predictionMode); + assert (rates.size() == + std::vector::size_type(pu.num_phases)); + + if (pu.phase_used[ BlackoilPhases::Aqua ]) { + const std::vector::size_type + i = pu.phase_pos[ BlackoilPhases::Aqua ]; + + rates[i] = p.WaterRate; + } + + if (pu.phase_used[ BlackoilPhases::Liquid ]) { + const std::vector::size_type + i = pu.phase_pos[ BlackoilPhases::Liquid ]; + + rates[i] = p.OilRate; + } + + if (pu.phase_used[ BlackoilPhases::Vapour ]) { + const std::vector::size_type + i = pu.phase_pos[ BlackoilPhases::Vapour ]; + + rates[i] = p.GasRate; + } + } + } // namespace SimFIBODetails + + template + void + SimulatorFullyImplicitBlackoilPolymer:: + Impl::computeRESV(const std::size_t step, + const Wells* wells, + const PolymerBlackoilState& x, + WellStateFullyImplicitBlackoil& xw) + { + typedef SimFIBODetails::WellMap WellMap; + + const std::vector& w_ecl = eclipse_state_->getSchedule()->getWells(step); + const WellMap& wmap = SimFIBODetails::mapWells(w_ecl); + + const std::vector& resv_prod = + SimFIBODetails::resvProducers(*wells, step, wmap); + + if (! resv_prod.empty()) { + const PhaseUsage& pu = props_.phaseUsage(); + const std::vector::size_type np = props_.numPhases(); + + rateConverter_.defineState(x.blackoilState()); + + std::vector distr (np); + std::vector hrates(np); + std::vector prates(np); + + for (std::vector::const_iterator + rp = resv_prod.begin(), e = resv_prod.end(); + rp != e; ++rp) + { + WellControls* ctrl = wells->ctrls[*rp]; + + // RESV control mode, all wells + { + const int rctrl = SimFIBODetails::resv_control(ctrl); + + if (0 <= rctrl) { + const std::vector::size_type off = (*rp) * np; + + // Convert to positive rates to avoid issues + // in coefficient calculations. + std::transform(xw.wellRates().begin() + (off + 0*np), + xw.wellRates().begin() + (off + 1*np), + prates.begin(), std::negate()); + + const int fipreg = 0; // Hack. Ignore FIP regions. + rateConverter_.calcCoeff(prates, fipreg, distr); + + well_controls_iset_distr(ctrl, rctrl, & distr[0]); + } + } + + // RESV control, WCONHIST wells. A bit of duplicate + // work, regrettably. + if (wells->name[*rp] != 0) { + WellMap::const_iterator i = wmap.find(wells->name[*rp]); + + if (i != wmap.end()) { + WellConstPtr wp = i->second; + + const WellProductionProperties& p = + wp->getProductionProperties(step); + + if (! p.predictionMode) { + // History matching (WCONHIST/RESV) + SimFIBODetails::historyRates(pu, p, hrates); + + const int fipreg = 0; // Hack. Ignore FIP regions. + rateConverter_.calcCoeff(hrates, fipreg, distr); + + // WCONHIST/RESV target is sum of all + // observed phase rates translated to + // reservoir conditions. Recall sign + // convention: Negative for producers. + const double target = + - std::inner_product(distr.begin(), distr.end(), + hrates.begin(), 0.0); + + well_controls_clear(ctrl); + well_controls_assert_number_of_phases(ctrl, int(np)); + + const int ok = + well_controls_add_new(RESERVOIR_RATE, target, + & distr[0], ctrl); + + if (ok != 0) { + xw.currentControls()[*rp] = 0; + well_controls_set_current(ctrl, 0); + } + } + } + } + } + } + } +} // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp new file mode 100644 index 000000000..f4d35ed51 --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp @@ -0,0 +1,474 @@ +/* + Copyright 2014 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 . +*/ + + +#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 + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +namespace Opm +{ + + + + namespace + { + static void outputStateVtk(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir); + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir); + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir); + } // anonymous namespace + + class SimulatorFullyImplicitCompressiblePolymer::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const DerivedGeology& geo, + const BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + std::shared_ptr eclipse_state, + EclipseWriter& output_writer, + Opm::DeckConstPtr& deck, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity); + + SimulatorReport run(SimulatorTimer& timer, + PolymerBlackoilState& state); + + private: + // Data. + + // Parameters for output. + bool output_; + bool output_vtk_; + std::string output_dir_; + int output_interval_; + // Parameters for well control + bool check_well_controls_; + int max_well_control_iterations_; + // Observed objects. + const UnstructuredGrid& grid_; + const BlackoilPropsAdInterface& props_; + const PolymerPropsAd& polymer_props_; + const RockCompressibility* rock_comp_props_; + std::shared_ptr eclipse_state_; + EclipseWriter& output_writer_; + Opm::DeckConstPtr& deck_; + NewtonIterationBlackoilInterface& linsolver_; + const double* gravity_; + // Solvers + DerivedGeology geo_; + // Misc. data + std::vector allcells_; + }; + + + + + SimulatorFullyImplicitCompressiblePolymer:: + SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const DerivedGeology& geo, + const BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + std::shared_ptr eclipse_state, + EclipseWriter& output_writer, + Opm::DeckConstPtr& deck, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity) + + { + pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, eclipse_state, output_writer, deck, linsolver, gravity)); + } + + + + + + SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer, + PolymerBlackoilState& state) + { + return pimpl_->run(timer, state); + } + + + + + + // \TODO: Treat bcs. + SimulatorFullyImplicitCompressiblePolymer::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const DerivedGeology& geo, + const BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + std::shared_ptr eclipse_state, + EclipseWriter& output_writer, + Opm::DeckConstPtr& deck, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + polymer_props_(polymer_props), + rock_comp_props_(rock_comp_props), + eclipse_state_(eclipse_state), + output_writer_(output_writer), + deck_(deck), + linsolver_(linsolver), + gravity_(gravity), + geo_(geo) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + + // Well control related init. + check_well_controls_ = param.getDefault("check_well_controls", false); + max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); + + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + SimulatorReport SimulatorFullyImplicitCompressiblePolymer::Impl::run(SimulatorTimer& timer, + PolymerBlackoilState& state) + { + WellStateFullyImplicitBlackoil prev_well_state; + // Initialisation. + std::vector porevol; + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } else { + computePorevolume(grid_, props_.porosity(), porevol); + } + std::vector initial_porevol = porevol; + + std::vector polymer_inflow_c(grid_.number_of_cells); + // Main simulation loop. + Opm::time::StopWatch solver_timer; + double stime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::string tstep_filename = output_dir_ + "/step_timing.txt"; + std::ofstream tstep_os(tstep_filename.c_str()); + + //Main simulation loop. + while (!timer.done()) { +#if 0 + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); + std::vector fractional_flows; + std::vector well_resflows_phase; + if (wells_) { + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + } + std::fstream tstep_os; + if (output_) { + std::string filename = output_dir_ + "/step_timing.param"; + tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); + } +#endif + // Report timestep and (optionally) write state to disk. + + step_timer.start(); + timer.report(std::cout); + + WellsManager wells_manager(eclipse_state_, + timer.currentStepNum(), + Opm::UgGridHelpers::numCells(grid_), + Opm::UgGridHelpers::globalCell(grid_), + Opm::UgGridHelpers::cartDims(grid_), + Opm::UgGridHelpers::dimensions(grid_), + Opm::UgGridHelpers::beginCellCentroids(grid_), + Opm::UgGridHelpers::cell2Faces(grid_), + Opm::UgGridHelpers::beginFaceCentroids(grid_), + props_.permeability()); + const Wells* wells = wells_manager.c_wells(); + WellStateFullyImplicitBlackoil well_state; + well_state.init(wells, state.blackoilState(), prev_well_state); + //Compute polymer inflow. + std::unique_ptr polymer_inflow_ptr; + if (deck_->hasKeyword("WPOLYMER")) { + if (wells_manager.c_wells() == 0) { + OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells."); + } + polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, *wells, Opm::UgGridHelpers::numCells(grid_))); + } else { + polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day, + 1.0*Opm::unit::day, + 0.0)); + } + std::vector polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_)); + polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), + timer.simulationTimeElapsed() + timer.currentStepLength(), + polymer_inflow_c); + + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + } + if (output_) { + if (timer.currentStepNum() == 0) { + output_writer_.writeInit(timer); + } + output_writer_.writeTimeStep(timer, state.blackoilState(), well_state.basicWellState()); + } + // Run solver. + solver_timer.start(); + FullyImplicitCompressiblePolymerSolver solver(grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells_manager.c_wells(), linsolver_); + solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c); + // Stop timer and report. + solver_timer.stop(); + const double st = solver_timer.secsSinceStart(); + std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; + + stime += st; + // Update pore volumes if rock is compressible. + if (rock_comp_props_ && rock_comp_props_->isActive()) { + initial_porevol = porevol; + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } +/* + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double polyinj = 0; + double polyprod = 0; + Opm::computeInjectedProduced(props_, polymer_props_, + state, + transport_src, polymer_inflow_c, timer.currentStepLength(), + injected, produced, + polyinj, polyprod); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + std::cout.precision(5); + const int width = 18; + std::cout << "\nMass balance report.\n"; + std::cout << " Injected reservoir volumes: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] << std::endl; + std::cout << " Produced reservoir volumes: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] << std::endl; + std::cout << " Total inj reservoir volumes: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] << std::endl; + std::cout << " Total prod reservoir volumes: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] << std::endl; +*/ + if (output_) { + SimulatorReport step_report; + step_report.pressure_time = st; + step_report.total_time = step_timer.secsSinceStart(); + step_report.reportParam(tstep_os); + } + ++timer; + prev_well_state = well_state; + } + // Write final simulation state. + if (output_) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state.basicWellState()); + } + + total_timer.stop(); + SimulatorReport report; + report.pressure_time = stime; + report.transport_time = 0.0; + report.total_time = total_timer.secsSinceStart(); + return report; + } + + + + namespace + { + + static void outputStateVtk(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.concentration(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + } + + + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::PolymerBlackoilState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.concentration(); + dm["surfvolume"] = &state.surfacevol(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + +#if 0 + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + watercut.write(os); + } +#endif + } + +} // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp new file mode 100644 index 000000000..df868c8af --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp @@ -0,0 +1,106 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED +#define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED + +#include +#include +#include + +struct UnstructuredGrid; +struct Wells; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class BlackoilPropsAdInterface; + class RockCompressibility; + class DerivedGeology; + class WellStateFullyImplicitBlackoil; + class WellsManager; + class EclipseWriter; + class EclipseState; + class NewtonIterationBlackoilInterface; + class SimulatorTimer; + class PolymerBlackoilState; + class PolymerPropsAd; + class PolymerInflowInterface; + struct SimulatorReport; + + /// Class collecting all necessary components for a two-phase simulation. + class SimulatorFullyImplicitCompressiblePolymer + { + public: + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] polymer_props polymer properties + /// \param[in] rock_comp_props if non-null, rock compressibility properties + /// \param[in] eclipse_state + /// \param[in] eclipse_writer + /// \param[in] deck + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const DerivedGeology& geo, + const BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + std::shared_ptr eclipse_state, + EclipseWriter& eclipse_writer, + Opm::DeckConstPtr& deck, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity); + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \return simulation report, with timing data + SimulatorReport run(SimulatorTimer& timer, + PolymerBlackoilState& state); + + private: + class Impl; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + std::shared_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp new file mode 100644 index 000000000..ba2fdfd06 --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp @@ -0,0 +1,513 @@ +/* + Copyright 2014 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 +#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 Opm +{ + + + + class SimulatorFullyImplicitTwophasePolymer::Impl + { + public: + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + LinearSolverInterface& linsolver, + WellsManager& wells_manager, + PolymerInflowInterface& polymer_inflow, + const double* gravity); + + SimulatorReport run(SimulatorTimer& timer, + PolymerState& state, + WellState& well_state); + + private: + + // Parameters for output. + bool output_; + bool output_vtk_; + std::string output_dir_; + int output_interval_; + // Parameters for well control + bool check_well_controls_; + int max_well_control_iterations_; + // Observed objects. + const UnstructuredGrid& grid_; + const IncompPropsAdInterface& props_; + const PolymerPropsAd& polymer_props_; + WellsManager& wells_manager_; + const Wells* wells_; + PolymerInflowInterface& polymer_inflow_; + // Solvers + FullyImplicitTwophasePolymerSolver solver_; + // Misc. data + std::vector allcells_; + }; + + + + + SimulatorFullyImplicitTwophasePolymer:: + SimulatorFullyImplicitTwophasePolymer(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + LinearSolverInterface& linsolver, + WellsManager& wells_manager, + PolymerInflowInterface& polymer_inflow, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, polymer_props, linsolver, wells_manager, polymer_inflow, gravity)); + } + + + + + + SimulatorReport SimulatorFullyImplicitTwophasePolymer::run(SimulatorTimer& timer, + PolymerState& state, + WellState& well_state) + { + return pimpl_->run(timer, state, well_state); + } + + + + static void outputStateVtk(const UnstructuredGrid& grid, + const Opm::PolymerState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/vtk_files"; + boost::filesystem::path fpath(vtkfilename.str()); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.concentration(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + } + + + static void outputStateMatlab(const UnstructuredGrid& grid, + const Opm::PolymerState& state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.concentration(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + watercut.write(os); + } +/* + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + OPM_THROW(std::runtime_error, "Failed to open " << fname); + } + wellreport.write(os); + } + */ + + static void outputWellStateMatlab(WellState& well_state, + const int step, + const std::string& output_dir) + { + Opm::DataMap dm; + dm["bhp"] = &well_state.bhp(); + dm["wellrates"] = &well_state.wellRates(); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first; + boost::filesystem::path fpath = fname.str(); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath); + } + fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; + std::ofstream file(fname.str().c_str()); + if (!file) { + OPM_THROW(std::runtime_error,"Failed to open " << fname.str()); + } + file.precision(15); + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + + + SimulatorFullyImplicitTwophasePolymer::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + LinearSolverInterface& linsolver, + WellsManager& wells_manager, + PolymerInflowInterface& polymer_inflow, + const double* gravity) + : grid_(grid), + props_(props), + polymer_props_(polymer_props), + wells_manager_(wells_manager), + wells_(wells_manager.c_wells()), + polymer_inflow_(polymer_inflow), + solver_(grid_, props_, polymer_props_, linsolver, *wells_manager.c_wells(), gravity) + + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_vtk_ = param.getDefault("output_vtk", true); + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + SimulatorReport SimulatorFullyImplicitTwophasePolymer::Impl::run(SimulatorTimer& timer, + PolymerState& state, + WellState& well_state) + { + + // Initialisation. + std::vector porevol; + Opm::computePorevolume(grid_, props_.porosity(), porevol); + + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + std::vector polymer_inflow_c(grid_.number_of_cells); + std::vector transport_src(grid_.number_of_cells); + + // Main simulation loop. + Opm::time::StopWatch solver_timer; + double stime = 0.0; + Opm::time::StopWatch step_timer; + Opm::time::StopWatch total_timer; + total_timer.start(); + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); +#if 0 + // These must be changed for three-phase. + double init_surfvol[2] = { 0.0 }; + double inplace_surfvol[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol); + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); +#endif + std::vector fractional_flows; + std::vector well_resflows_phase; + std::fstream tstep_os; + if (output_) { + std::string filename = output_dir_ + "/step_timing.param"; + tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); + } + while (!timer.done()) { + // Report timestep and (optionally) write state to disk. + step_timer.start(); + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + // outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); + + } + + SimulatorReport sreport; + + bool well_control_passed = !check_well_controls_; + int well_control_iteration = 0; + do { + // Run solver. + const double current_time = timer.currentTime(); + double stepsize = timer.currentStepLength(); + polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c); + solver_timer.start(); + std::vector initial_pressure = state.pressure(); + solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src); + + // Stop timer and report. + solver_timer.stop(); + const double st = solver_timer.secsSinceStart(); + std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; + + stime += st; + sreport.pressure_time = st; + + // Optionally, check if well controls are satisfied. + if (check_well_controls_) { + Opm::computePhaseFlowRatesPerWell(*wells_, + well_state.perfRates(), + fractional_flows, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { + OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); + } + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; + } + } + } while (!well_control_passed); + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double polyinj = 0; + double polyprod = 0; + + Opm::computeInjectedProduced(props_, polymer_props_, + state, + transport_src, polymer_inflow_c, timer.currentStepLength(), + injected, produced, + polyinj, polyprod); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + std::cout.precision(5); + const int width = 18; + std::cout << "\nMass balance report.\n"; + std::cout << " Injected reservoir volumes: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] << std::endl; + std::cout << " Produced reservoir volumes: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] << std::endl; + std::cout << " Total inj reservoir volumes: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] << std::endl; + std::cout << " Total prod reservoir volumes: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] << std::endl; + + + // Update pore volumes if rock is compressible. + + // The reports below are geared towards two phases only. +#if 0 + // Report mass balances. + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + Opm::computeInjectedProduced(props_, state, transport_src, stepsize, + injected, produced); + Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + std::cout.precision(5); + const int width = 18; + std::cout << "\nMass balance report.\n"; + std::cout << " Injected surface volumes: " + << std::setw(width) << injected[0] + << std::setw(width) << injected[1] << std::endl; + std::cout << " Produced surface volumes: " + << std::setw(width) << produced[0] + << std::setw(width) << produced[1] << std::endl; + std::cout << " Total inj surface volumes: " + << std::setw(width) << tot_injected[0] + << std::setw(width) << tot_injected[1] << std::endl; + std::cout << " Total prod surface volumes: " + << std::setw(width) << tot_produced[0] + << std::setw(width) << tot_produced[1] << std::endl; + const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0], + init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] }; + std::cout << " Initial - inplace + inj - prod: " + << std::setw(width) << balance[0] + << std::setw(width) << balance[1] + << std::endl; + std::cout << " Relative mass error: " + << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0]) + << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1]) + << std::endl; + std::cout.precision(8); + + // Make well reports. + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, + state.pressure(), state.surfacevol(), state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } +#endif + sreport.total_time = step_timer.secsSinceStart(); + if (output_) { + sreport.reportParam(tstep_os); + + if (output_vtk_) { + outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); + } + outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); + outputWaterCut(watercut, output_dir_); +#if 0 + outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } +#endif + tstep_os.close(); + } + + // advance to next timestep before reporting at this location + ++timer; + + // write an output file for later inspection + } + + total_timer.stop(); + + SimulatorReport report; + report.pressure_time = stime; + report.transport_time = 0.0; + report.total_time = total_timer.secsSinceStart(); + return report; + } + + + + + +} // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp new file mode 100644 index 000000000..d6cd62e6b --- /dev/null +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp @@ -0,0 +1,97 @@ +/* + Copyright 2014 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 . +*/ + +#ifndef OPM_SIMULATORFULLYIMPLICITTWOPHASEPOLYMER_HEADER_INCLUDED +#define OPM_SIMULATORFULLYIMPLICITTWOPHASEPOLYMER_HEADER_INCLUDED + +#include +#include + +struct UnstructuredGrid; +struct Wells; +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class IncompPropsAdInterface; + class LinearSolverInterface; + class SimulatorTimer; + class PolymerState; + class PolymerPropsAd; + class PolymerInflowInterface; + class WellsManager; + class WellState; + struct SimulatorReport; + + /// Class collecting all necessary components for a two-phase simulation. + class SimulatorFullyImplicitTwophasePolymer + { + public: + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) + /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] polymer_props polymer properties + /// \param[in] linsolver linear solver + /// \param[in] well_manager well manager, may manage no (null) wells + /// \param[in] polymer_inflow polymer influx. + /// \param[in] gravity if non-null, gravity vector + SimulatorFullyImplicitTwophasePolymer(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + LinearSolverInterface& linsolver, + WellsManager& wells_manager, + PolymerInflowInterface& polymer_inflow, + const double* gravity); + + /// Run the simulation. + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + /// \return simulation report, with timing data + SimulatorReport run(SimulatorTimer& timer, + PolymerState& state, + WellState& well_state); + + private: + class Impl; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + boost::shared_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORFULLYIMPLICITTWOPHASEPOLYMER_HEADER_INCLUDED