OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #561 from atgeirr/refactor-extended-mains

Browse Source 
Refactor main functions of all Flow variants.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-12-22 15:20:59 +01:00
commit 5315a0dfba
13 changed files with 962 additions and 1169 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@ -141,7 +141,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/fastSparseProduct.hpp
opm/autodiff/DuneMatrix.hpp
opm/autodiff/ExtractParallelGridInformationToISTL.hpp
opm/autodiff/flowMain.hpp
opm/autodiff/FlowMain.hpp
opm/autodiff/FlowMainPolymer.hpp
opm/autodiff/FlowMainSolvent.hpp
opm/autodiff/GeoProps.hpp
opm/autodiff/GridHelpers.hpp
opm/autodiff/GridInit.hpp

5
examples/flow.cpp
View File

@ -27,7 +27,7 @@
#include <opm/core/grid.h>
#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
#include <opm/autodiff/flowMain.hpp>
#include <opm/autodiff/FlowMain.hpp>
// ----------------- Main program -----------------
@ -37,5 +37,6 @@ main(int argc, char** argv)
typedef UnstructuredGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoil<Grid> Simulator;
return Opm::flowMain<Grid, Simulator>(argc, argv);
Opm::FlowMain<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}

5
examples/flow_mpi.cpp
View File

@ -27,7 +27,7 @@
#include <dune/grid/CpGrid.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
#include <opm/autodiff/flowMain.hpp>
#include <opm/autodiff/FlowMain.hpp>
// ----------------- Main program -----------------
@ -37,5 +37,6 @@ main(int argc, char** argv)
typedef Dune::CpGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoil<Grid> Simulator;
return Opm::flowMain<Grid, Simulator>(argc, argv);
Opm::FlowMain<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}

5
examples/flow_multisegment.cpp
View File

@ -27,7 +27,7 @@
#include <opm/core/grid.h>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp>
#include <opm/autodiff/flowMain.hpp>
#include <opm/autodiff/FlowMain.hpp>
// ----------------- Main program -----------------
@ -37,5 +37,6 @@ main(int argc, char** argv)
typedef UnstructuredGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoilMultiSegment<Grid> Simulator;
return Opm::flowMain<Grid, Simulator>(argc, argv);
Opm::FlowMain<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}

5
examples/flow_multisegment_mpi.cpp
View File

@ -27,7 +27,7 @@
#include <dune/grid/CpGrid.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment.hpp>
#include <opm/autodiff/flowMain.hpp>
#include <opm/autodiff/FlowMain.hpp>
// ----------------- Main program -----------------
@ -37,5 +37,6 @@ main(int argc, char** argv)
typedef Dune::CpGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoilMultiSegment<Grid> Simulator;
return Opm::flowMain<Grid, Simulator>(argc, argv);
Opm::FlowMain<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}

320
examples/flow_polymer.cpp
View File

@ -1,5 +1,5 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014 STATOIL ASA.
This file is part of the Open Porous Media project (OPM).
@ -17,324 +17,26 @@
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/pressure/FlowBCManager.hpp>
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/thresholdPressures.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/moduleVersion.hpp>
#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/PolymerInflow.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#include <memory>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
#include <opm/autodiff/FlowMainPolymer.hpp>
// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
using namespace Opm;
typedef UnstructuredGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoilPolymer<Grid> Simulator;
{
std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow-Polymer (version " << version << ")"
<< std::string(18 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow-Polymer is a simulator for fully implicit three-phase, *\n";
std::cout << "* four-component (black-oil + polymer) flow, and is part of OPM. *\n";
std::cout << "* For more information see http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
// Read parameters, see if a deck was specified on the command line.
std::cout << "--------------- Reading parameters ---------------" << std::endl;
parameter::ParameterGroup param(argc, argv, false);
if (!param.unhandledArguments().empty()) {
if (param.unhandledArguments().size() != 1) {
OPM_THROW(std::runtime_error, "You can only specify a single input deck on the command line.");
} else {
param.insertParameter("deck_filename", param.unhandledArguments()[0]);
}
}
// We must have an input deck. Grid and props will be read from that.
if (!param.has("deck_filename")) {
std::cerr << "This program must be run with an input deck.\n"
"Specify the deck filename either\n"
" a) as a command line argument by itself\n"
" b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
" c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
OPM_THROW(std::runtime_error, "Input deck required.");
}
std::shared_ptr<GridManager> grid;
std::shared_ptr<BlackoilPropertiesFromDeck> props;
std::shared_ptr<BlackoilPropsAdFromDeck> new_props;
std::shared_ptr<RockCompressibility> 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<std::string>("deck_filename");
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::string output_dir;
if (output) {
// Create output directory if needed.
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
// Write simulation parameters.
param.writeParam(output_dir + "/simulation.param");
}
std::string logFile = output_dir + "/LOGFILE.txt";
Opm::ParserPtr parser(new Opm::Parser());
{
std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
std::shared_ptr<Opm::CounterLog> counterLog = std::make_shared<Opm::CounterLog>(Opm::Log::DefaultMessageTypes);
Opm::OpmLog::addBackend( "STREAM" , streamLog );
Opm::OpmLog::addBackend( "COUNTER" , counterLog );
}
Opm::DeckConstPtr deck;
std::shared_ptr<EclipseState> eclipseState;
Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
try {
deck = parser->parseFile(deck_filename , parseMode);
Opm::checkDeck(deck);
eclipseState.reset(new Opm::EclipseState(deck , parseMode));
}
catch (const std::invalid_argument& e) {
std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl;
std::cerr << "Exception caught: " << e.what() << std::endl;
return EXIT_FAILURE;
}
// Grid init
std::vector<double> porv = eclipseState->getDoubleGridProperty("PORV")->getData();
grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv));
auto &cGrid = *grid->c_grid();
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
// Rock and fluid init
std::vector<int> compressedToCartesianIdx;
Opm::createGlobalCellArray(*grid->c_grid(), compressedToCartesianIdx);
typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
props.reset(new BlackoilPropertiesFromDeck( deck, eclipseState, materialLawManager,
Opm::UgGridHelpers::numCells(cGrid),
Opm::UgGridHelpers::globalCell(cGrid),
Opm::UgGridHelpers::cartDims(cGrid),
param));
new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, materialLawManager, cGrid));
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);
state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0);
} else {
initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
// The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
if (deck->hasKeyword("SWATINIT")) {
const int nc = grid->c_grid()->number_of_cells;
std::vector<int> cells(nc);
for (int c = 0; c < nc; ++c) { cells[c] = c; }
std::vector<double> pc = state.saturation();
props->capPress(nc, state.saturation().data(), cells.data(), pc.data(),NULL);
new_props->setSwatInitScaling(state.saturation(),pc);
}
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
// Solver for Newton iterations.
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
if (param.getDefault("use_cpr", true)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param));
} else {
fis_solver.reset(new NewtonIterationBlackoilSimple(param));
}
Opm::ScheduleConstPtr schedule = eclipseState->getSchedule();
Opm::TimeMapConstPtr timeMap(schedule->getTimeMap());
SimulatorTimer simtimer;
// initialize variables
simtimer.init(timeMap);
if (polymer){
if (!use_wpolymer) {
OPM_MESSAGE("Warning: simulate polymer injection without WPOLYMER.");
} else {
if (param.has("polymer_start_days")) {
OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck."
"You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
}
}
} else {
if (use_wpolymer) {
OPM_MESSAGE("Warning: use WPOLYMER in a non-polymer scenario.");
}
}
bool use_local_perm = param.getDefault("use_local_perm", true);
Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, use_local_perm, grav);
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, *grid->c_grid(), state, *props, gravity[2]);
std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, *grid->c_grid(), maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geology.nnc(), maxDp);
threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
Opm::BlackoilOutputWriter
outputWriter(cGrid, param, eclipseState, pu,
new_props->permeability());
SimulatorFullyImplicitBlackoilPolymer<UnstructuredGrid>
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,
deck->hasKeyword("PLYSHLOG"),
deck->hasKeyword("SHRATE"),
eclipseState,
outputWriter,
deck,
threshold_pressures);
if (!schedule->initOnly()){
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
SimulatorReport fullReport = simulator.run(simtimer, state);
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.report(std::cout);
if (output) {
std::string filename = output_dir + "/walltime.txt";
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
fullReport.reportParam(tot_os);
warnIfUnusedParams(param);
}
} else {
outputWriter.writeInit( simtimer );
std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
}
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
Opm::FlowMainPolymer<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}

406
examples/flow_solvent.cpp
View File

@ -25,414 +25,18 @@
#endif // HAVE_CONFIG_H
#include <dune/common/version.hh>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
#include <dune/common/parallel/mpihelper.hh>
#else
#include <dune/common/mpihelper.hh>
#endif
#if HAVE_DUNE_CORNERPOINT && WANT_DUNE_CORNERPOINTGRID
#define USE_DUNE_CORNERPOINTGRID 1
#include <dune/grid/CpGrid.hpp>
#include <dune/grid/common/GridAdapter.hpp>
#else
#undef USE_DUNE_CORNERPOINTGRID
#endif
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
#include <opm/autodiff/RedistributeDataHandles.hpp>
#include <opm/autodiff/moduleVersion.hpp>
#include <opm/core/utility/share_obj.hpp>
#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#include <memory>
#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <numeric>
#include <cstdlib>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
#include <opm/autodiff/FlowMainSolvent.hpp>
// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
using namespace Opm;
#if USE_DUNE_CORNERPOINTGRID
// Must ensure an instance of the helper is created to initialise MPI.
const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv);
const int mpi_rank = mpi_helper.rank();
const int mpi_size = mpi_helper.size();
#else
// default values for serial run
const int mpi_rank = 0;
const int mpi_size = 1;
#endif
typedef UnstructuredGrid Grid;
typedef Opm::SimulatorFullyImplicitBlackoilSolvent<Grid> Simulator;
// Write parameters used for later reference. (only if rank is zero)
const bool output_cout = ( mpi_rank == 0 );
if (output_cout)
{
std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow-Solvent (version " << version << ")"
<< std::string(18 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow-Solvent is a simulator for fully implicit three-phase, *\n";
std::cout << "* four-component (black-oil + solvent) flow, and is part of OPM. *\n";
std::cout << "* For more information see http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
// Read parameters, see if a deck was specified on the command line.
if ( output_cout )
{
std::cout << "--------------- Reading parameters ---------------" << std::endl;
}
parameter::ParameterGroup param(argc, argv, false, output_cout);
if( !output_cout )
{
param.disableOutput();
}
if (!param.unhandledArguments().empty()) {
if (param.unhandledArguments().size() != 1) {
std::cerr << "You can only specify a single input deck on the command line.\n";
return EXIT_FAILURE;
} else {
param.insertParameter("deck_filename", param.unhandledArguments()[0]);
}
}
// We must have an input deck. Grid and props will be read from that.
if (!param.has("deck_filename")) {
std::cerr << "This program must be run with an input deck.\n"
"Specify the deck filename either\n"
" a) as a command line argument by itself\n"
" b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
" c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
return EXIT_FAILURE;
}
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
std::string deck_filename = param.get<std::string>("deck_filename");
// Write parameters used for later reference. (only if rank is zero)
bool output = ( mpi_rank == 0 ) && 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 (...) {
std::cerr << "Creating directories failed: " << fpath << std::endl;
return EXIT_FAILURE;
}
// Write simulation parameters.
param.writeParam(output_dir + "/simulation.param");
}
std::string logFile = output_dir + "/LOGFILE.txt";
Opm::ParserPtr parser(new Opm::Parser());
{
std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
std::shared_ptr<Opm::CounterLog> counterLog = std::make_shared<Opm::CounterLog>(Opm::Log::DefaultMessageTypes);
Opm::OpmLog::addBackend( "STREAM" , streamLog );
Opm::OpmLog::addBackend( "COUNTER" , counterLog );
}
Opm::DeckConstPtr deck;
std::shared_ptr<EclipseState> eclipseState;
Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
try {
deck = parser->parseFile(deck_filename, parseMode);
Opm::checkDeck(deck);
eclipseState.reset(new Opm::EclipseState(deck , parseMode));
}
catch (const std::invalid_argument& e) {
std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl;
std::cerr << "Exception caught: " << e.what() << std::endl;
return EXIT_FAILURE;
}
std::vector<double> porv = eclipseState->getDoubleGridProperty("PORV")->getData();
// UnstructuredGrid as grid manager
typedef UnstructuredGrid Grid;
GridManager gridManager( eclipseState->getEclipseGrid(), porv );
const Grid& grid = *(gridManager.c_grid());
// Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
if (param.has("output_interval")) {
int output_interval = param.get<int>("output_interval");
IOConfigPtr ioConfig = eclipseState->getIOConfig();
ioConfig->overrideRestartWriteInterval((size_t)output_interval);
}
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
std::vector<int> compressedToCartesianIdx;
Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
// Rock and fluid init
BlackoilPropertiesFromDeck props( deck, eclipseState, materialLawManager,
Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
param);
BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid );
SolventPropsAdFromDeck solvent_props( deck, eclipseState, Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid));
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
RockCompressibility rock_comp(deck, eclipseState);
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
BlackoilSolventState state;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
Opm::UgGridHelpers::numFaces(grid),
Opm::UgGridHelpers::faceCells(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, param, gravity[2], state);
initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state);
enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
const int numPhases = props.numPhases();
const int numCells = Opm::UgGridHelpers::numCells(grid);
for (int c = 0; c < numCells; ++c) {
state.gasoilratio()[c] = state.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
/ state.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
}
}
} else if (deck->hasKeyword("EQUIL") && props.numPhases() == 3) {
state.init(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::numFaces(grid),
props.numPhases());
const double grav = param.getDefault("gravity", unit::gravity);
initStateEquil(grid, props, deck, eclipseState, grav, state);
state.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
} else {
initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::numFaces(grid),
Opm::UgGridHelpers::faceCells(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, deck, gravity[2], state);
}
// The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
if (deck->hasKeyword("SWATINIT")) {
const int numCells = Opm::UgGridHelpers::numCells(grid);
std::vector<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> pc = state.saturation();
props.capPress(numCells, state.saturation().data(), cells.data(), pc.data(),NULL);
new_props.setSwatInitScaling(state.saturation(),pc);
}
//state.solvent_saturation()[49] = 0.1;
//state.saturation()[49*3 + 2] -= 0.1;
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
#if USE_DUNE_CORNERPOINTGRID
if(output_cout)
{
std::cout << std::endl << "Warning: use of local perm is not yet implemented for CpGrid!" << std::endl << std::endl;
}
const bool use_local_perm = false;
#else
const bool use_local_perm = param.getDefault("use_local_perm", true);
#endif
DerivedGeology geoprops(grid, new_props, eclipseState, use_local_perm, grav);
boost::any parallel_information;
// At this point all properties and state variables are correctly initialized
// If there are more than one processors involved, we now repartition the grid
// and initilialize new properties and states for it.
if( mpi_size > 1 )
{
Opm::distributeGridAndData( grid, deck, eclipseState, state, new_props, geoprops, materialLawManager, parallel_information, use_local_perm );
}
// create output writer after grid is distributed, otherwise the parallel output
// won't work correctly since we need to create a mapping from the distributed to
// the global view
Opm::BlackoilOutputWriter outputWriter(grid, param, eclipseState, pu, new_props.permeability() );
// Solver for Newton iterations.
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
if (param.getDefault("use_interleaved", true)) {
fis_solver.reset(new NewtonIterationBlackoilInterleaved(param));
} else if (param.getDefault("use_cpr", true)) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param));
} else {
fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information));
}
Opm::ScheduleConstPtr schedule = eclipseState->getSchedule();
Opm::TimeMapConstPtr timeMap(schedule->getTimeMap());
SimulatorTimer simtimer;
// initialize variables
simtimer.init(timeMap);
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geoprops.nnc(), maxDp);
threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
SimulatorFullyImplicitBlackoilSolvent< Grid > simulator(param,
grid,
geoprops,
new_props,
solvent_props,
rock_comp.isActive() ? &rock_comp : 0,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL"),
eclipseState,
outputWriter,
deck,
threshold_pressures,
deck->hasKeyword("SOLVENT") );
if (!schedule->initOnly()){
if( output_cout )
{
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
}
SimulatorReport fullReport = simulator.run(simtimer, state);
if( output_cout )
{
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.reportFullyImplicit(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);
}
} else {
outputWriter.writeInit( simtimer );
if ( output_cout )
{
std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
}
}
Opm::FlowMainSolvent<Grid, Simulator> mainfunc;
return mainfunc.execute(argc, argv);
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
return EXIT_FAILURE;
}

675
opm/autodiff/FlowMain.hpp Normal file
View File

@ -0,0 +1,675 @@
/*
Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 IRIS AS
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLOWMAIN_HEADER_INCLUDED
#define OPM_FLOWMAIN_HEADER_INCLUDED
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/version.hh>
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
#include <dune/common/parallel/mpihelper.hh>
#else
#include <dune/common/mpihelper.hh>
#endif
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/RedistributeDataHandles.hpp>
#include <opm/autodiff/moduleVersion.hpp>
#include <opm/core/utility/share_obj.hpp>
#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <memory>
#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <numeric>
#include <cstdlib>
namespace Opm
{
/// This class encapsulates the setup and running of
/// a simulator based on an input deck.
template <class Implementation, class Grid, class Simulator>
class FlowMainBase
{
public:
/// This is the main function of Flow.
/// It runs a complete simulation, with the given grid and
/// simulator classes, based on user command-line input. The
/// content of this function used to be in the main() function of
/// flow.cpp.
int execute(int argc, char** argv)
try {
// Setup.
asImpl().setupParallelism(argc, argv);
asImpl().printStartupMessage();
const bool ok = asImpl().setupParameters(argc, argv);
if (!ok) {
return EXIT_FAILURE;
}
asImpl().setupOutput();
asImpl().readDeckInput();
asImpl().setupGridAndProps();
asImpl().setupState();
asImpl().distributeData();
asImpl().setupOutputWriter();
asImpl().setupLinearSolver();
asImpl().createSimulator();
// Run.
return asImpl().runSimulator();
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
return EXIT_FAILURE;
}
protected:
// ------------ Types ------------
typedef BlackoilPropsAdFromDeck FluidProps;
typedef FluidProps::MaterialLawManager MaterialLawManager;
typedef typename Simulator::ReservoirState ReservoirState;
// ------------ Data members ------------
// The comments indicate in which method the
// members first occur.
// setupParallelism()
bool output_cout_ = false;
bool must_distribute_ = false;
// setupParameters()
parameter::ParameterGroup param_;
// setupOutput()
bool output_to_files_ = false;
std::string output_dir_ = std::string("output");
// readDeckInput()
std::shared_ptr<const Deck> deck_;
std::shared_ptr<EclipseState> eclipse_state_;
// setupGridAndProps()
std::unique_ptr<GridInit<Grid>> grid_init_;
std::shared_ptr<MaterialLawManager> material_law_manager_;
std::unique_ptr<FluidProps> fluidprops_;
std::unique_ptr<RockCompressibility> rock_comp_;
std::array<double, 3> gravity_;
bool use_local_perm_ = true;
std::unique_ptr<DerivedGeology> geoprops_;
// setupState()
ReservoirState state_;
std::vector<double> threshold_pressures_;
// distributeData()
boost::any parallel_information_;
// setupOutputWriter()
std::unique_ptr<BlackoilOutputWriter> output_writer_;
// setupLinearSolver
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
// createSimulator()
std::unique_ptr<Simulator> simulator_;
// ------------ Methods ------------
// Set up MPI and OpenMP.
// Writes to:
// output_cout_
// must_distribute_
void setupParallelism(int argc, char** argv)
{
// MPI setup.
// Must ensure an instance of the helper is created to initialise MPI.
// For a build without MPI the Dune::FakeMPIHelper is used, so rank will
// be 0 and size 1.
const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv);
const int mpi_rank = mpi_helper.rank();
const int mpi_size = mpi_helper.size();
output_cout_ = ( mpi_rank == 0 );
must_distribute_ = ( mpi_size > 1 );
#ifdef _OPENMP
// OpenMP setup.
if (!getenv("OMP_NUM_THREADS")) {
// Default to at most 4 threads, regardless of
// number of cores (unless ENV(OMP_NUM_THREADS) is defined)
int num_cores = omp_get_num_procs();
int num_threads = std::min(4, num_cores);
omp_set_num_threads(num_threads);
}
#pragma omp parallel
if (omp_get_thread_num() == 0) {
// omp_get_num_threads() only works as expected within a parallel region.
const int num_omp_threads = omp_get_num_threads();
if (mpi_size == 1) {
std::cout << "OpenMP using " << num_omp_threads << " threads." << std::endl;
} else {
std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank << "." << std::endl;
}
}
#endif
}
// Print startup message if on output rank.
void printStartupMessage()
{
if (output_cout_) {
const std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow (version " << version << ")"
<< std::string(26 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n";
std::cout << "* and is part of OPM. For more information see: *\n";
std::cout << "* http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
}
// Read parameters, see if a deck was specified on the command line, and if
// it was, insert it into parameters.
// Writes to:
// param_
// Returns true if ok, false if not.
bool setupParameters(int argc, char** argv)
{
// Read parameters.
if ( output_cout_ ) {
std::cout << "--------------- Reading parameters ---------------" << std::endl;
}
param_ = parameter::ParameterGroup(argc, argv, false, output_cout_);
// See if a deck was specified on the command line.
if (!param_.unhandledArguments().empty()) {
if (param_.unhandledArguments().size() != 1) {
std::cerr << "You can only specify a single input deck on the command line.\n";
return false;
} else {
param_.insertParameter("deck_filename", param_.unhandledArguments()[0]);
}
}
// We must have an input deck. Grid and props will be read from that.
if (!param_.has("deck_filename")) {
std::cerr << "This program must be run with an input deck.\n"
"Specify the deck filename either\n"
" a) as a command line argument by itself\n"
" b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
" c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
return false;
}
return true;
}
// Set output_to_files_ and set/create output dir. Write parameter file.
// Writes to:
// output_to_files_
// output_dir_
// Throws std::runtime_error if failed to create (if requested) output dir.
void setupOutput()
{
// Write parameters used for later reference. (only if rank is zero)
output_to_files_ = output_cout_ && param_.getDefault("output", true);
if (output_to_files_) {
// 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");
}
}
// Parser the input and creates the Deck and EclipseState objects.
// Writes to:
// deck_
// eclipse_state_
// May throw if errors are encountered, here configured to be somewhat tolerant.
void readDeckInput()
{
// Create Parser
std::string logFile = output_dir_ + "/LOGFILE.txt";
ParserPtr parser(new Parser());
{
std::shared_ptr<StreamLog> streamLog = std::make_shared<StreamLog>(logFile , Log::DefaultMessageTypes);
std::shared_ptr<CounterLog> counterLog = std::make_shared<CounterLog>(Log::DefaultMessageTypes);
OpmLog::addBackend( "STREAM" , streamLog );
OpmLog::addBackend( "COUNTER" , counterLog );
}
// Create Deck and EclipseState.
try {
std::string deck_filename = param_.get<std::string>("deck_filename");
ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
deck_ = parser->parseFile(deck_filename, parseMode);
checkDeck(deck_);
eclipse_state_.reset(new EclipseState(deck_, parseMode));
}
catch (const std::invalid_argument& e) {
std::cerr << "Failed to create valid EclipseState object. See logfile: " << logFile << std::endl;
std::cerr << "Exception caught: " << e.what() << std::endl;
throw;
}
// Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
if (param_.has("output_interval")) {
const int output_interval = param_.get<int>("output_interval");
IOConfigPtr ioConfig = eclipse_state_->getIOConfig();
ioConfig->overrideRestartWriteInterval(static_cast<size_t>(output_interval));
}
}
// Create grid and property objects.
// Writes to:
// grid_init_
// material_law_manager_
// fluidprops_
// rock_comp_
// gravity_
// use_local_perm_
// geoprops_
void setupGridAndProps()
{
// Create grid.
const std::vector<double>& porv = eclipse_state_->getDoubleGridProperty("PORV")->getData();
grid_init_.reset(new GridInit<Grid>(deck_, eclipse_state_, porv));
const Grid& grid = grid_init_->grid();
// Create material law manager.
std::vector<int> compressedToCartesianIdx;
Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
material_law_manager_.reset(new MaterialLawManager());
material_law_manager_->initFromDeck(deck_, eclipse_state_, compressedToCartesianIdx);
// Rock and fluid properties.
fluidprops_.reset(new BlackoilPropsAdFromDeck(deck_, eclipse_state_, material_law_manager_, grid));
// Rock compressibility.
rock_comp_.reset(new RockCompressibility(deck_, eclipse_state_));
// Gravity.
assert(UgGridHelpers::dimensions(grid) == 3);
gravity_.fill(0.0);
gravity_[2] = deck_->hasKeyword("NOGRAV")
? param_.getDefault("gravity", 0.0)
: param_.getDefault("gravity", unit::gravity);
// Geological properties
use_local_perm_ = param_.getDefault("use_local_perm", use_local_perm_);
geoprops_.reset(new DerivedGeology(grid, *fluidprops_, eclipse_state_, use_local_perm_, gravity_.data()));
}
// Initialise the reservoir state. Updated fluid props for SWATINIT.
// Writes to:
// state_
// threshold_pressures_
// fluidprops_ (if SWATINIT is used)
void setupState()
{
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck_);
const Grid& grid = grid_init_->grid();
// Need old-style fluid object for init purposes (only).
BlackoilPropertiesFromDeck props( deck_, eclipse_state_, material_law_manager_,
Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
param_);
// Init state variables (saturation and pressure).
if (param_.has("init_saturation")) {
initStateBasic(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
Opm::UgGridHelpers::numFaces(grid),
Opm::UgGridHelpers::faceCells(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, param_, gravity_[2], state_);
initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state_);
enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
const int numPhases = props.numPhases();
const int numCells = Opm::UgGridHelpers::numCells(grid);
for (int c = 0; c < numCells; ++c) {
state_.gasoilratio()[c] = state_.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
/ state_.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
}
}
} else if (deck_->hasKeyword("EQUIL") && props.numPhases() == 3) {
state_.init(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::numFaces(grid),
props.numPhases());
initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], state_);
state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
} else {
initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::numFaces(grid),
Opm::UgGridHelpers::faceCells(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, deck_, gravity_[2], state_);
}
// Threshold pressures.
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), state_, props, gravity_[2]);
threshold_pressures_ = thresholdPressures(deck_, eclipse_state_, grid, maxDp);
std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipse_state_, geoprops_->nnc(), maxDp);
threshold_pressures_.insert(threshold_pressures_.end(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
// The capillary pressure is scaled in fluidprops_ to match the scaled capillary pressure in props.
if (deck_->hasKeyword("SWATINIT")) {
const int numCells = Opm::UgGridHelpers::numCells(grid);
std::vector<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> pc = state_.saturation();
props.capPress(numCells, state_.saturation().data(), cells.data(), pc.data(), nullptr);
fluidprops_->setSwatInitScaling(state_.saturation(), pc);
}
}
// Distribute the grid, properties and state.
// Writes to:
// grid_init_->grid()
// state_
// fluidprops_
// geoprops_
// material_law_manager_
// parallel_information_
void distributeData()
{
// At this point all properties and state variables are correctly initialized
// If there are more than one processors involved, we now repartition the grid
// and initilialize new properties and states for it.
if (must_distribute_) {
distributeGridAndData(grid_init_->grid(), deck_, eclipse_state_, state_, *fluidprops_, *geoprops_,
material_law_manager_, parallel_information_, use_local_perm_);
}
}
// Setup output writer.
// Writes to:
// output_writer_
void setupOutputWriter()
{
// create output writer after grid is distributed, otherwise the parallel output
// won't work correctly since we need to create a mapping from the distributed to
// the global view
output_writer_.reset(new BlackoilOutputWriter(grid_init_->grid(),
param_,
eclipse_state_,
Opm::phaseUsageFromDeck(deck_),
fluidprops_->permeability()));
}
// Setup linear solver.
// Writes to:
// fis_solver_
void setupLinearSolver()
{
const std::string cprSolver = "cpr";
const std::string interleavedSolver = "interleaved";
const std::string directSolver = "direct";
const std::string flowDefaultSolver = interleavedSolver;
std::shared_ptr<const Opm::SimulationConfig> simCfg = eclipse_state_->getSimulationConfig();
std::string solver_approach = flowDefaultSolver;
if (param_.has("solver_approach")) {
solver_approach = param_.get<std::string>("solver_approach");
} else {
if (simCfg->useCPR()) {
solver_approach = cprSolver;
}
}
if (solver_approach == cprSolver) {
fis_solver_.reset(new NewtonIterationBlackoilCPR(param_, parallel_information_));
} else if (solver_approach == interleavedSolver) {
fis_solver_.reset(new NewtonIterationBlackoilInterleaved(param_, parallel_information_));
} else if (solver_approach == directSolver) {
fis_solver_.reset(new NewtonIterationBlackoilSimple(param_, parallel_information_));
} else {
OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized.");
}
}
// Run the simulator.
// Returns EXIT_SUCCESS if it does not throw.
int runSimulator()
{
Opm::ScheduleConstPtr schedule = eclipse_state_->getSchedule();
Opm::TimeMapConstPtr timeMap(schedule->getTimeMap());
SimulatorTimer simtimer;
// initialize variables
simtimer.init(timeMap);
if (!schedule->initOnly()) {
if (output_cout_) {
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
}
SimulatorReport fullReport = simulator_->run(simtimer, state_);
if (output_cout_) {
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.reportFullyImplicit(std::cout);
if (param_.anyUnused()) {
// This allows a user to catch typos and misunderstandings in the
// use of simulator parameters.
std::cout << "-------------------- Unused parameters: --------------------\n";
param_.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
if (output_to_files_) {
std::string filename = output_dir_ + "/walltime.txt";
std::fstream tot_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
fullReport.reportParam(tot_os);
}
} else {
output_writer_->writeInit( simtimer );
if (output_cout_) {
std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
}
}
return EXIT_SUCCESS;
}
// Access the most-derived class used for
// static polymorphism (CRTP).
Implementation& asImpl()
{
return static_cast<Implementation&>(*this);
}
}; // class FlowMainBase
// The FlowMain class is the basic black-oil simulator case.
template <class Grid, class Simulator>
class FlowMain : public FlowMainBase<FlowMain<Grid, Simulator>, Grid, Simulator>
{
protected:
using Base = FlowMainBase<FlowMain<Grid, Simulator>, Grid, Simulator>;
friend Base;
// Create simulator instance.
// Writes to:
// simulator_
void createSimulator()
{
// Create the simulator instance.
Base::simulator_.reset(new Simulator(Base::param_,
Base::grid_init_->grid(),
*Base::geoprops_,
*Base::fluidprops_,
Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr,
*Base::fis_solver_,
Base::gravity_.data(),
Base::deck_->hasKeyword("DISGAS"),
Base::deck_->hasKeyword("VAPOIL"),
Base::eclipse_state_,
*Base::output_writer_,
Base::threshold_pressures_));
}
};
} // namespace Opm
#endif // OPM_FLOWMAIN_HEADER_INCLUDED

137
opm/autodiff/FlowMainPolymer.hpp Normal file
View File

@ -0,0 +1,137 @@
/*
Copyright 2014, 2015 STATOIL ASA.
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLOWMAINPOLYMER_HEADER_INCLUDED
#define OPM_FLOWMAINPOLYMER_HEADER_INCLUDED
#include <opm/autodiff/FlowMain.hpp>
#include <opm/polymer/PolymerProperties.hpp>
#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
namespace Opm
{
// The FlowMainPolymer class is for a black-oil simulator with polymer.
template <class Grid, class Simulator>
class FlowMainPolymer : public FlowMainBase<FlowMainPolymer<Grid, Simulator>, Grid, Simulator>
{
protected:
using Base = FlowMainBase<FlowMainPolymer<Grid, Simulator>, Grid, Simulator>;
friend Base;
// Set in setupGridAndProps()
std::unique_ptr<PolymerProperties> polymer_props_legacy_; // Held by reference in polymer_props_
std::unique_ptr<PolymerPropsAd> polymer_props_;
// ------------ Methods ------------
// Print startup message if on output rank.
void printStartupMessage()
{
if (Base::output_cout_) {
const std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow-Polymer (version " << version << ")"
<< std::string(18 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow-Polymer is a simulator for fully implicit three-phase, *\n";
std::cout << "* four-component (black-oil + polymer) flow, and is part of OPM. *\n";
std::cout << "* For more information see http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
}
// Set up grid and property objects, by calling base class
// version and then creating polymer property objects.
// Writes to:
// polymer_props_legacy_
// polymer_props_
void setupGridAndProps()
{
Base::setupGridAndProps();
polymer_props_legacy_.reset(new PolymerProperties(Base::deck_, Base::eclipse_state_));
polymer_props_.reset(new PolymerPropsAd(*polymer_props_legacy_));
}
// Setup linear solver.
// Writes to:
// fis_solver_
void setupLinearSolver()
{
OPM_MESSAGE("Caution: polymer solver currently only works with direct linear solver.");
Base::fis_solver_.reset(new NewtonIterationBlackoilSimple(Base::param_, Base::parallel_information_));
}
// Create simulator instance.
// Writes to:
// simulator_
void createSimulator()
{
// Create the simulator instance.
Base::simulator_.reset(new Simulator(Base::param_,
Base::grid_init_->grid(),
*Base::geoprops_,
*Base::fluidprops_,
*polymer_props_,
Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr,
*Base::fis_solver_,
Base::gravity_.data(),
Base::deck_->hasKeyword("DISGAS"),
Base::deck_->hasKeyword("VAPOIL"),
Base::deck_->hasKeyword("POLYMER"),
Base::deck_->hasKeyword("PLYSHLOG"),
Base::deck_->hasKeyword("SHRATE"),
Base::eclipse_state_,
*Base::output_writer_,
Base::deck_,
Base::threshold_pressures_));
}
};
} // namespace Opm
#endif // OPM_FLOWMAINPOLYMER_HEADER_INCLUDED

117
opm/autodiff/FlowMainSolvent.hpp Normal file
View File

@ -0,0 +1,117 @@
/*
Copyright 2015 IRIS AS
Copyright 2015 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLOWMAINSOLVENT_HEADER_INCLUDED
#define OPM_FLOWMAINSOLVENT_HEADER_INCLUDED
#include <opm/autodiff/FlowMain.hpp>
#include <opm/autodiff/SolventPropsAdFromDeck.hpp>
namespace Opm
{
// The FlowMainSolvent class is for a black-oil simulator with solvent.
template <class Grid, class Simulator>
class FlowMainSolvent : public FlowMainBase<FlowMainSolvent<Grid, Simulator>, Grid, Simulator>
{
protected:
using Base = FlowMainBase<FlowMainSolvent<Grid, Simulator>, Grid, Simulator>;
friend Base;
// Set in setupGridAndProps()
std::unique_ptr<SolventPropsAdFromDeck> solvent_props_;
// ------------ Methods ------------
// Print startup message if on output rank.
void printStartupMessage()
{
if (Base::output_cout_) {
const std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow-Solvent (version " << version << ")"
<< std::string(18 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow-Solvent is a simulator for fully implicit three-phase, *\n";
std::cout << "* four-component (black-oil + solvent) flow, and is part of OPM. *\n";
std::cout << "* For more information see http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
}
// Set up grid and property objects, by calling base class
// version and then creating solvent property object.
void setupGridAndProps()
{
Base::setupGridAndProps();
const Grid& grid = Base::grid_init_->grid();
solvent_props_.reset(new SolventPropsAdFromDeck(Base::deck_,
Base::eclipse_state_,
UgGridHelpers::numCells(grid),
UgGridHelpers::globalCell(grid)));
}
// Create simulator instance.
// Writes to:
// simulator_
void createSimulator()
{
// Create the simulator instance.
Base::simulator_.reset(new Simulator(Base::param_,
Base::grid_init_->grid(),
*Base::geoprops_,
*Base::fluidprops_,
*solvent_props_,
Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr,
*Base::fis_solver_,
Base::gravity_.data(),
Base::deck_->hasKeyword("DISGAS"),
Base::deck_->hasKeyword("VAPOIL"),
Base::eclipse_state_,
*Base::output_writer_,
Base::deck_,
Base::threshold_pressures_,
Base::deck_->hasKeyword("SOLVENT")));
}
};
} // namespace Opm
#endif // OPM_FLOWMAINSOLVENT_HEADER_INCLUDED

2
opm/autodiff/GeoProps.hpp
View File

@ -153,7 +153,7 @@ namespace Opm
}
// opmfil is hardcoded to be true. i.e the pinch processor is never used
if (~opmfil && eclgrid->isPinchActive()) {
if (!opmfil && eclgrid->isPinchActive()) {
const double minpv = eclgrid->getMinpvValue();
const double thickness = eclgrid->getPinchThresholdThickness();
auto transMode = eclgrid->getPinchOption();

449
opm/autodiff/flowMain.hpp
View File

@ -1,449 +0,0 @@
/*
Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2015 IRIS AS
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_FLOWMAIN_HEADER_INCLUDED
#define OPM_FLOWMAIN_HEADER_INCLUDED
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/version.hh>
#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3)
#include <dune/common/parallel/mpihelper.hh>
#else
#include <dune/common/mpihelper.hh>
#endif
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/RedistributeDataHandles.hpp>
#include <opm/autodiff/moduleVersion.hpp>
#include <opm/core/utility/share_obj.hpp>
#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <memory>
#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <numeric>
#include <cstdlib>
namespace Opm
{
/// Calling this will print the unused parameters, if any.
/// This allows a user to catch typos and misunderstandings in the
/// use of simulator parameters.
inline void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
/// This is the main function of Flow.
/// It runs a complete simulation, with the given grid and
/// simulator classes, based on user command-line input. The
/// content of this function used to be in the main() function of
/// flow.cpp.
template <class Grid, class Simulator>
inline int flowMain(int argc, char** argv)
try {
using namespace Opm;
// Must ensure an instance of the helper is created to initialise MPI.
// For a build without MPI the Dune::FakeMPIHelper is used, so rank will
// be 0 and size 1.
const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv);
const int mpi_rank = mpi_helper.rank();
const int mpi_size = mpi_helper.size();
// Write parameters used for later reference. (only if rank is zero)
const bool output_cout = ( mpi_rank == 0 );
if (output_cout)
{
std::string version = moduleVersionName();
std::cout << "**********************************************************************\n";
std::cout << "* *\n";
std::cout << "* This is Flow (version " << version << ")"
<< std::string(26 - version.size(), ' ') << "*\n";
std::cout << "* *\n";
std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n";
std::cout << "* and is part of OPM. For more information see: *\n";
std::cout << "* http://opm-project.org *\n";
std::cout << "* *\n";
std::cout << "**********************************************************************\n\n";
}
#ifdef _OPENMP
if (!getenv("OMP_NUM_THREADS")) {
//Default to at most 4 threads, regardless of
//number of cores (unless ENV(OMP_NUM_THREADS) is defined)
int num_cores = omp_get_num_procs();
int num_threads = std::min(4, num_cores);
omp_set_num_threads(num_threads);
}
#pragma omp parallel
if (omp_get_thread_num() == 0){
//opm_get_num_threads() only works as expected within a parallel region.
std::cout << "OpenMP using " << omp_get_num_threads() << " threads." << std::endl;
}
#endif
// Read parameters, see if a deck was specified on the command line.
if ( output_cout )
{
std::cout << "--------------- Reading parameters ---------------" << std::endl;
}
parameter::ParameterGroup param(argc, argv, false, output_cout);
if( !output_cout )
{
param.disableOutput();
}
if (!param.unhandledArguments().empty()) {
if (param.unhandledArguments().size() != 1) {
std::cerr << "You can only specify a single input deck on the command line.\n";
return EXIT_FAILURE;
} else {
param.insertParameter("deck_filename", param.unhandledArguments()[0]);
}
}
// We must have an input deck. Grid and props will be read from that.
if (!param.has("deck_filename")) {
std::cerr << "This program must be run with an input deck.\n"
"Specify the deck filename either\n"
" a) as a command line argument by itself\n"
" b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
" c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
return EXIT_FAILURE;
}
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
std::string deck_filename = param.get<std::string>("deck_filename");
// Write parameters used for later reference. (only if rank is zero)
bool output = ( mpi_rank == 0 ) && 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 (...) {
std::cerr << "Creating directories failed: " << fpath << std::endl;
return EXIT_FAILURE;
}
// Write simulation parameters.
param.writeParam(output_dir + "/simulation.param");
}
std::string logFile = output_dir + "/LOGFILE.txt";
Opm::ParserPtr parser(new Opm::Parser());
{
std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
std::shared_ptr<Opm::CounterLog> counterLog = std::make_shared<Opm::CounterLog>(Opm::Log::DefaultMessageTypes);
Opm::OpmLog::addBackend( "STREAM" , streamLog );
Opm::OpmLog::addBackend( "COUNTER" , counterLog );
}
Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
Opm::DeckConstPtr deck;
std::shared_ptr<EclipseState> eclipseState;
try {
deck = parser->parseFile(deck_filename, parseMode);
Opm::checkDeck(deck);
eclipseState.reset(new Opm::EclipseState(deck , parseMode));
}
catch (const std::invalid_argument& e) {
std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl;
std::cerr << "Exception caught: " << e.what() << std::endl;
return EXIT_FAILURE;
}
std::vector<double> porv = eclipseState->getDoubleGridProperty("PORV")->getData();
GridInit<Grid> grid_init(deck, eclipseState, porv);
auto&& grid = grid_init.grid();
// Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
if (param.has("output_interval")) {
int output_interval = param.get<int>("output_interval");
IOConfigPtr ioConfig = eclipseState->getIOConfig();
ioConfig->overrideRestartWriteInterval((size_t)output_interval);
}
const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
std::vector<int> compressedToCartesianIdx;
Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager;
auto materialLawManager = std::make_shared<MaterialLawManager>();
materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx);
// Rock and fluid init
BlackoilPropertiesFromDeck props( deck, eclipseState, materialLawManager,
Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
param);
BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid );
// check_well_controls = param.getDefault("check_well_controls", false);
// max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Rock compressibility.
RockCompressibility rock_comp(deck, eclipseState);
// Gravity.
gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
typename Simulator::ReservoirState state;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::cartDims(grid),
Opm::UgGridHelpers::numFaces(grid),
Opm::UgGridHelpers::faceCells(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, param, gravity[2], state);
initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state);
enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
const int numPhases = props.numPhases();
const int numCells = Opm::UgGridHelpers::numCells(grid);
for (int c = 0; c < numCells; ++c) {
state.gasoilratio()[c] = state.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
/ state.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
}
}
} else if (deck->hasKeyword("EQUIL") && props.numPhases() == 3) {
state.init(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::numFaces(grid),
props.numPhases());
const double grav = param.getDefault("gravity", unit::gravity);
initStateEquil(grid, props, deck, eclipseState, grav, state);
state.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
} else {
initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
Opm::UgGridHelpers::globalCell(grid),
Opm::UgGridHelpers::numFaces(grid),
Opm::UgGridHelpers::faceCells(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
Opm::UgGridHelpers::beginCellCentroids(grid),
Opm::UgGridHelpers::dimensions(grid),
props, deck, gravity[2], state);
}
// The capillary pressure is scaled in new_props to match the scaled capillary pressure in props.
if (deck->hasKeyword("SWATINIT")) {
const int numCells = Opm::UgGridHelpers::numCells(grid);
std::vector<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> pc = state.saturation();
props.capPress(numCells, state.saturation().data(), cells.data(), pc.data(),NULL);
new_props.setSwatInitScaling(state.saturation(),pc);
}
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
const double *grav = use_gravity ? &gravity[0] : 0;
const bool use_local_perm = param.getDefault("use_local_perm", true);
DerivedGeology geoprops(grid, new_props, eclipseState, use_local_perm, grav);
boost::any parallel_information;
// At this point all properties and state variables are correctly initialized
// If there are more than one processors involved, we now repartition the grid
// and initilialize new properties and states for it.
if( mpi_size > 1 )
{
Opm::distributeGridAndData( grid, deck, eclipseState, state, new_props, geoprops, materialLawManager, parallel_information, use_local_perm );
}
// create output writer after grid is distributed, otherwise the parallel output
// won't work correctly since we need to create a mapping from the distributed to
// the global view
Opm::BlackoilOutputWriter outputWriter(grid, param, eclipseState, pu, new_props.permeability() );
// Solver for Newton iterations.
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
{
const std::string cprSolver = "cpr";
const std::string interleavedSolver = "interleaved";
const std::string directSolver = "direct";
const std::string flowDefaultSolver = interleavedSolver;
std::shared_ptr<const Opm::SimulationConfig> simCfg = eclipseState->getSimulationConfig();
std::string solver_approach = flowDefaultSolver;
if (param.has("solver_approach")) {
solver_approach = param.get<std::string>("solver_approach");
} else {
if (simCfg->useCPR()) {
solver_approach = cprSolver;
}
}
if (solver_approach == cprSolver) {
fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information));
} else if (solver_approach == interleavedSolver) {
fis_solver.reset(new NewtonIterationBlackoilInterleaved(param, parallel_information));
} else if (solver_approach == directSolver) {
fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information));
} else {
OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized.");
}
}
Opm::ScheduleConstPtr schedule = eclipseState->getSchedule();
Opm::TimeMapConstPtr timeMap(schedule->getTimeMap());
SimulatorTimer simtimer;
// initialize variables
simtimer.init(timeMap);
std::map<std::pair<int, int>, double> maxDp;
computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
Simulator simulator(param,
grid,
geoprops,
new_props,
rock_comp.isActive() ? &rock_comp : 0,
*fis_solver,
grav,
deck->hasKeyword("DISGAS"),
deck->hasKeyword("VAPOIL"),
eclipseState,
outputWriter,
threshold_pressures);
if (!schedule->initOnly()){
if( output_cout )
{
std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush;
}
SimulatorReport fullReport = simulator.run(simtimer, state);
if( output_cout )
{
std::cout << "\n\n================ End of simulation ===============\n\n";
fullReport.reportFullyImplicit(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);
}
} else {
outputWriter.writeInit( simtimer );
if ( output_cout )
{
std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
}
}
return EXIT_SUCCESS;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
return EXIT_FAILURE;
}
} // namespace Opm
#endif // OPM_FLOWMAIN_HEADER_INCLUDED

1
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
View File

@ -112,6 +112,7 @@ namespace Opm
}
polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum()));
} else {
OPM_MESSAGE("Warning: simulating with no WPOLYMER in deck (no polymer will be injected).");
polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day,
1.0*Opm::unit::day,
0.0));