From abf5443f334ed4eae8e22f939a129dfec9fb33a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 09:56:09 +0100 Subject: [PATCH 01/19] Replace flowMain() function with FlowMain class. --- examples/flow.cpp | 3 ++- examples/flow_mpi.cpp | 3 ++- examples/flow_multisegment.cpp | 3 ++- examples/flow_multisegment_mpi.cpp | 3 ++- opm/autodiff/flowMain.hpp | 7 ++++++- 5 files changed, 14 insertions(+), 5 deletions(-) diff --git a/examples/flow.cpp b/examples/flow.cpp index 9e601ff66..f47bee9f0 100644 --- a/examples/flow.cpp +++ b/examples/flow.cpp @@ -37,5 +37,6 @@ main(int argc, char** argv) typedef UnstructuredGrid Grid; typedef Opm::SimulatorFullyImplicitBlackoil Simulator; - return Opm::flowMain(argc, argv); + Opm::FlowMain mainfunc; + return mainfunc.execute(argc, argv); } diff --git a/examples/flow_mpi.cpp b/examples/flow_mpi.cpp index c68e850d4..3ae0b57ab 100644 --- a/examples/flow_mpi.cpp +++ b/examples/flow_mpi.cpp @@ -37,5 +37,6 @@ main(int argc, char** argv) typedef Dune::CpGrid Grid; typedef Opm::SimulatorFullyImplicitBlackoil Simulator; - return Opm::flowMain(argc, argv); + Opm::FlowMain mainfunc; + return mainfunc.execute(argc, argv); } diff --git a/examples/flow_multisegment.cpp b/examples/flow_multisegment.cpp index feec4346d..447a963c7 100644 --- a/examples/flow_multisegment.cpp +++ b/examples/flow_multisegment.cpp @@ -37,5 +37,6 @@ main(int argc, char** argv) typedef UnstructuredGrid Grid; typedef Opm::SimulatorFullyImplicitBlackoilMultiSegment Simulator; - return Opm::flowMain(argc, argv); + Opm::FlowMain mainfunc; + return mainfunc.execute(argc, argv); } diff --git a/examples/flow_multisegment_mpi.cpp b/examples/flow_multisegment_mpi.cpp index 5525e0658..64b6fdc85 100644 --- a/examples/flow_multisegment_mpi.cpp +++ b/examples/flow_multisegment_mpi.cpp @@ -37,5 +37,6 @@ main(int argc, char** argv) typedef Dune::CpGrid Grid; typedef Opm::SimulatorFullyImplicitBlackoilMultiSegment Simulator; - return Opm::flowMain(argc, argv); + Opm::FlowMain mainfunc; + return mainfunc.execute(argc, argv); } diff --git a/opm/autodiff/flowMain.hpp b/opm/autodiff/flowMain.hpp index 76647abd4..8b3b08cfa 100644 --- a/opm/autodiff/flowMain.hpp +++ b/opm/autodiff/flowMain.hpp @@ -119,7 +119,10 @@ namespace Opm /// content of this function used to be in the main() function of /// flow.cpp. template - inline int flowMain(int argc, char** argv) + class FlowMain + { + public: + int execute(int argc, char** argv) try { using namespace Opm; @@ -443,6 +446,8 @@ namespace Opm return EXIT_FAILURE; } + }; // class FlowMain + } // namespace Opm From e3ceac44a6b4a9c21936954893e12f79d3400518 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 09:58:00 +0100 Subject: [PATCH 02/19] Renamed flowMain.hpp -> FlowMain.hpp. --- CMakeLists_files.cmake | 2 +- examples/flow.cpp | 2 +- examples/flow_mpi.cpp | 2 +- examples/flow_multisegment.cpp | 2 +- examples/flow_multisegment_mpi.cpp | 2 +- opm/autodiff/{flowMain.hpp => FlowMain.hpp} | 0 6 files changed, 5 insertions(+), 5 deletions(-) rename opm/autodiff/{flowMain.hpp => FlowMain.hpp} (100%) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 61a6a2aba..543bf187e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -141,7 +141,7 @@ 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/GeoProps.hpp opm/autodiff/GridHelpers.hpp opm/autodiff/GridInit.hpp diff --git a/examples/flow.cpp b/examples/flow.cpp index f47bee9f0..53a5b2d6d 100644 --- a/examples/flow.cpp +++ b/examples/flow.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include // ----------------- Main program ----------------- diff --git a/examples/flow_mpi.cpp b/examples/flow_mpi.cpp index 3ae0b57ab..e3f65b970 100644 --- a/examples/flow_mpi.cpp +++ b/examples/flow_mpi.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include // ----------------- Main program ----------------- diff --git a/examples/flow_multisegment.cpp b/examples/flow_multisegment.cpp index 447a963c7..cfb113bf3 100644 --- a/examples/flow_multisegment.cpp +++ b/examples/flow_multisegment.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include // ----------------- Main program ----------------- diff --git a/examples/flow_multisegment_mpi.cpp b/examples/flow_multisegment_mpi.cpp index 64b6fdc85..e9a78fcdf 100644 --- a/examples/flow_multisegment_mpi.cpp +++ b/examples/flow_multisegment_mpi.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include // ----------------- Main program ----------------- diff --git a/opm/autodiff/flowMain.hpp b/opm/autodiff/FlowMain.hpp similarity index 100% rename from opm/autodiff/flowMain.hpp rename to opm/autodiff/FlowMain.hpp From f7c96b4fc5de092547caacb1e6804cc3b77997c9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 10:51:03 +0100 Subject: [PATCH 03/19] Split out function setupParallelism(). --- opm/autodiff/FlowMain.hpp | 86 +++++++++++++++++++++++---------------- 1 file changed, 51 insertions(+), 35 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 8b3b08cfa..32e46734e 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -124,19 +124,9 @@ namespace Opm public: int execute(int argc, char** argv) try { - using namespace Opm; + setupParallelism(argc, argv); - // 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) + if (output_cout_) { std::string version = moduleVersionName(); std::cout << "**********************************************************************\n"; @@ -151,29 +141,14 @@ namespace Opm 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 ) + if ( output_cout_ ) { std::cout << "--------------- Reading parameters ---------------" << std::endl; } - parameter::ParameterGroup param(argc, argv, false, output_cout); - if( !output_cout ) + parameter::ParameterGroup param(argc, argv, false, output_cout_); + if( !output_cout_ ) { param.disableOutput(); } @@ -203,7 +178,7 @@ namespace Opm std::string deck_filename = param.get("deck_filename"); // Write parameters used for later reference. (only if rank is zero) - bool output = ( mpi_rank == 0 ) && param.getDefault("output", true); + bool output = output_cout_ && param.getDefault("output", true); std::string output_dir; if (output) { // Create output directory if needed. @@ -346,7 +321,7 @@ namespace Opm // 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 ) + if( must_distribute_ ) { Opm::distributeGridAndData( grid, deck, eclipseState, state, new_props, geoprops, materialLawManager, parallel_information, use_local_perm ); } @@ -412,7 +387,7 @@ namespace Opm threshold_pressures); if (!schedule->initOnly()){ - if( output_cout ) + if( output_cout_ ) { std::cout << "\n\n================ Starting main simulation loop ===============\n" << std::flush; @@ -420,7 +395,7 @@ namespace Opm SimulatorReport fullReport = simulator.run(simtimer, state); - if( output_cout ) + if( output_cout_ ) { std::cout << "\n\n================ End of simulation ===============\n\n"; fullReport.reportFullyImplicit(std::cout); @@ -434,7 +409,7 @@ namespace Opm } } else { outputWriter.writeInit( simtimer ); - if ( output_cout ) + if ( output_cout_ ) { std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush; } @@ -446,6 +421,47 @@ namespace Opm return EXIT_FAILURE; } + + + private: + + // ------------ Data members ------------ + + bool output_cout_ = false; + bool must_distribute_ = false; + + // ------------ Methods ------------ + + 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) { + // opm_get_num_threads() only works as expected within a parallel region. + std::cout << "OpenMP using " << omp_get_num_threads() << " threads." << std::endl; + } +#endif + } + + }; // class FlowMain From c1bd7670c676692583743c07ee714da81f5645f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 11:12:43 +0100 Subject: [PATCH 04/19] Refine OpenMP number-of-threads output. Now each MPI rank will report, making it easier to see how MPI and OpenMP interacts w.r.t. processes and threads. --- opm/autodiff/FlowMain.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 32e46734e..fdb082d0a 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -455,8 +455,13 @@ namespace Opm } #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; + // 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 } From 966ace2c536929c9f3b9cf5581a0908e065e8f11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 11:20:53 +0100 Subject: [PATCH 05/19] Split out printStartupMessage() method. --- opm/autodiff/FlowMain.hpp | 57 ++++++++++++++++++++++++++++----------- 1 file changed, 41 insertions(+), 16 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index fdb082d0a..63f5321a2 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -125,34 +125,20 @@ namespace Opm int execute(int argc, char** argv) try { setupParallelism(argc, argv); - - 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"; - } + printStartupMessage(); // 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"; @@ -427,11 +413,27 @@ namespace Opm // ------------ Data members ------------ + + + + bool output_cout_ = false; bool must_distribute_ = false; + + + + // ------------ Methods ------------ + + + + + // Set up MPI and OpenMP. + // Writes to: + // output_cout_ + // must_distribute_ void setupParallelism(int argc, char** argv) { // MPI setup. @@ -467,6 +469,29 @@ namespace Opm } + + + + // 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"; + } + } + + + }; // class FlowMain From 64eef0609a3252101854862930a4d177bfb287f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 13:20:32 +0100 Subject: [PATCH 06/19] Split out setupParameters() method. --- opm/autodiff/FlowMain.hpp | 110 +++++++++++++++++++++----------------- 1 file changed, 60 insertions(+), 50 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 63f5321a2..4e606cf23 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -126,50 +126,20 @@ namespace Opm try { setupParallelism(argc, argv); printStartupMessage(); - - // 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=, or\n" - " c) as a parameter in a parameter file (.param or .xml) passed to the program.\n"; + const bool ok = setupParameters(argc, argv); + if (!ok) { 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("deck_filename"); // Write parameters used for later reference. (only if rank is zero) - bool output = output_cout_ && param.getDefault("output", true); + bool output = output_cout_ && param_.getDefault("output", true); std::string output_dir; if (output) { // Create output directory if needed. output_dir = - param.getDefault("output_dir", std::string("output")); + param_.getDefault("output_dir", std::string("output")); boost::filesystem::path fpath(output_dir); try { create_directories(fpath); @@ -179,7 +149,7 @@ namespace Opm return EXIT_FAILURE; } // Write simulation parameters. - param.writeParam(output_dir + "/simulation.param"); + param_.writeParam(output_dir + "/simulation.param"); } std::string logFile = output_dir + "/LOGFILE.txt"; @@ -195,6 +165,7 @@ namespace Opm Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); Opm::DeckConstPtr deck; std::shared_ptr eclipseState; + std::string deck_filename = param_.get("deck_filename"); try { deck = parser->parseFile(deck_filename, parseMode); Opm::checkDeck(deck); @@ -211,8 +182,8 @@ namespace Opm 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("output_interval"); + if (param_.has("output_interval")) { + int output_interval = param_.get("output_interval"); IOConfigPtr ioConfig = eclipseState->getIOConfig(); ioConfig->overrideRestartWriteInterval((size_t)output_interval); } @@ -231,7 +202,7 @@ namespace Opm Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid), Opm::UgGridHelpers::cartDims(grid), - param); + param_); BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid ); // check_well_controls = param.getDefault("check_well_controls", false); @@ -245,7 +216,7 @@ namespace Opm typename Simulator::ReservoirState state; // Init state variables (saturation and pressure). - if (param.has("init_saturation")) { + if (param_.has("init_saturation")) { initStateBasic(Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid), Opm::UgGridHelpers::cartDims(grid), @@ -254,7 +225,7 @@ namespace Opm Opm::UgGridHelpers::beginFaceCentroids(grid), Opm::UgGridHelpers::beginCellCentroids(grid), Opm::UgGridHelpers::dimensions(grid), - props, param, gravity[2], state); + props, param_, gravity[2], state); initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state); @@ -271,7 +242,7 @@ namespace Opm state.init(Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::numFaces(grid), props.numPhases()); - const double grav = param.getDefault("gravity", unit::gravity); + 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 { @@ -299,7 +270,7 @@ namespace Opm 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); + 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; @@ -315,7 +286,7 @@ namespace Opm // 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() ); + Opm::BlackoilOutputWriter outputWriter(grid, param_, eclipseState, pu, new_props.permeability() ); // Solver for Newton iterations. std::unique_ptr fis_solver; @@ -328,8 +299,8 @@ namespace Opm std::shared_ptr simCfg = eclipseState->getSimulationConfig(); std::string solver_approach = flowDefaultSolver; - if (param.has("solver_approach")) { - solver_approach = param.get("solver_approach"); + if (param_.has("solver_approach")) { + solver_approach = param_.get("solver_approach"); } else { if (simCfg->useCPR()) { solver_approach = cprSolver; @@ -337,11 +308,11 @@ namespace Opm } if (solver_approach == cprSolver) { - fis_solver.reset(new NewtonIterationBlackoilCPR(param, parallel_information)); + fis_solver.reset(new NewtonIterationBlackoilCPR(param_, parallel_information)); } else if (solver_approach == interleavedSolver) { - fis_solver.reset(new NewtonIterationBlackoilInterleaved(param, parallel_information)); + fis_solver.reset(new NewtonIterationBlackoilInterleaved(param_, parallel_information)); } else if (solver_approach == directSolver) { - fis_solver.reset(new NewtonIterationBlackoilSimple(param, parallel_information)); + fis_solver.reset(new NewtonIterationBlackoilSimple(param_, parallel_information)); } else { OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized."); } @@ -359,7 +330,7 @@ namespace Opm computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]); std::vector threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp); - Simulator simulator(param, + Simulator simulator(param_, grid, geoprops, new_props, @@ -391,7 +362,7 @@ namespace Opm 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); + warnIfUnusedParams(param_); } } else { outputWriter.writeInit( simtimer ); @@ -419,6 +390,7 @@ namespace Opm bool output_cout_ = false; bool must_distribute_ = false; + parameter::ParameterGroup param_; @@ -492,6 +464,44 @@ namespace Opm + + + // 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=, or\n" + " c) as a parameter in a parameter file (.param or .xml) passed to the program.\n"; + return false; + } + return true; + } + + }; // class FlowMain From 4fae28f3d3ccc01d0e052c8fb3ce3ad062431191 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 13:37:20 +0100 Subject: [PATCH 07/19] Split out setupOutput() method. --- opm/autodiff/FlowMain.hpp | 61 +++++++++++++++++++++++---------------- 1 file changed, 36 insertions(+), 25 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 4e606cf23..01fd07ecb 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -130,29 +130,9 @@ namespace Opm if (!ok) { return EXIT_FAILURE; } + setupOutput(); - double gravity[3] = { 0.0 }; - - // Write parameters used for later reference. (only if rank is zero) - bool output = output_cout_ && 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"; + std::string logFile = output_dir_ + "/LOGFILE.txt"; Opm::ParserPtr parser(new Opm::Parser()); { std::shared_ptr streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); @@ -212,6 +192,7 @@ namespace Opm RockCompressibility rock_comp(deck, eclipseState); // Gravity. + double gravity[3] = { 0.0 }; gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; typename Simulator::ReservoirState state; @@ -358,8 +339,8 @@ namespace Opm fullReport.reportFullyImplicit(std::cout); } - if (output) { - std::string filename = output_dir + "/walltime.txt"; + 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); warnIfUnusedParams(param_); @@ -391,7 +372,8 @@ namespace Opm bool output_cout_ = false; bool must_distribute_ = false; parameter::ParameterGroup param_; - + bool output_to_files_ = false; + std::string output_dir_ = "output"; @@ -502,6 +484,35 @@ namespace Opm } + + + + // 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"); + } + } + + }; // class FlowMain From 7a536570dc637925cfd9881e12a5990797545bf1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 9 Dec 2015 15:06:44 +0100 Subject: [PATCH 08/19] Split out readDeckInput() method. --- opm/autodiff/FlowMain.hpp | 128 +++++++++++++++++++++----------------- 1 file changed, 72 insertions(+), 56 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 01fd07ecb..de65698c1 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -131,69 +131,38 @@ namespace Opm return EXIT_FAILURE; } setupOutput(); + readDeckInput(); - std::string logFile = output_dir_ + "/LOGFILE.txt"; - Opm::ParserPtr parser(new Opm::Parser()); - { - std::shared_ptr streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); - std::shared_ptr counterLog = std::make_shared(Opm::Log::DefaultMessageTypes); - - Opm::OpmLog::addBackend( "STREAM" , streamLog ); - Opm::OpmLog::addBackend( "COUNTER" , counterLog ); - } - - Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); - Opm::DeckConstPtr deck; - std::shared_ptr eclipseState; - std::string deck_filename = param_.get("deck_filename"); - 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 porv = eclipseState->getDoubleGridProperty("PORV")->getData(); - GridInit grid_init(deck, eclipseState, porv); + std::vector porv = eclipse_state_->getDoubleGridProperty("PORV")->getData(); + GridInit grid_init(deck_, eclipse_state_, 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("output_interval"); - IOConfigPtr ioConfig = eclipseState->getIOConfig(); - ioConfig->overrideRestartWriteInterval((size_t)output_interval); - } - - const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); + const PhaseUsage pu = Opm::phaseUsageFromDeck(deck_); std::vector compressedToCartesianIdx; Opm::createGlobalCellArray(grid, compressedToCartesianIdx); typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager; auto materialLawManager = std::make_shared(); - materialLawManager->initFromDeck(deck, eclipseState, compressedToCartesianIdx); + materialLawManager->initFromDeck(deck_, eclipse_state_, compressedToCartesianIdx); // Rock and fluid init - BlackoilPropertiesFromDeck props( deck, eclipseState, materialLawManager, + BlackoilPropertiesFromDeck props( deck_, eclipse_state_, materialLawManager, Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid), Opm::UgGridHelpers::cartDims(grid), param_); - BlackoilPropsAdFromDeck new_props( deck, eclipseState, materialLawManager, grid ); + BlackoilPropsAdFromDeck new_props( deck_, eclipse_state_, 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); + RockCompressibility rock_comp(deck_, eclipse_state_); // Gravity. double gravity[3] = { 0.0 }; - gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; + gravity[2] = deck_->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; typename Simulator::ReservoirState state; // Init state variables (saturation and pressure). @@ -219,12 +188,12 @@ namespace Opm / state.surfacevol()[c*numPhases + pu.phase_pos[Oil]]; } } - } else if (deck->hasKeyword("EQUIL") && props.numPhases() == 3) { + } 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); + initStateEquil(grid, props, deck_, eclipse_state_, grav, state); state.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); } else { initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid), @@ -234,12 +203,12 @@ namespace Opm Opm::UgGridHelpers::beginFaceCentroids(grid), Opm::UgGridHelpers::beginCellCentroids(grid), Opm::UgGridHelpers::dimensions(grid), - props, deck, gravity[2], state); + 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")) { + if (deck_->hasKeyword("SWATINIT")) { const int numCells = Opm::UgGridHelpers::numCells(grid); std::vector cells(numCells); for (int c = 0; c < numCells; ++c) { cells[c] = c; } @@ -253,7 +222,7 @@ namespace Opm const bool use_local_perm = param_.getDefault("use_local_perm", true); - DerivedGeology geoprops(grid, new_props, eclipseState, use_local_perm, grav); + DerivedGeology geoprops(grid, new_props, eclipse_state_, use_local_perm, grav); boost::any parallel_information; // At this point all properties and state variables are correctly initialized @@ -261,13 +230,13 @@ namespace Opm // and initilialize new properties and states for it. if( must_distribute_ ) { - Opm::distributeGridAndData( grid, deck, eclipseState, state, new_props, geoprops, materialLawManager, parallel_information, use_local_perm ); + Opm::distributeGridAndData( grid, deck_, eclipse_state_, 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() ); + Opm::BlackoilOutputWriter outputWriter(grid, param_, eclipse_state_, pu, new_props.permeability() ); // Solver for Newton iterations. std::unique_ptr fis_solver; @@ -277,7 +246,7 @@ namespace Opm const std::string directSolver = "direct"; const std::string flowDefaultSolver = interleavedSolver; - std::shared_ptr simCfg = eclipseState->getSimulationConfig(); + std::shared_ptr simCfg = eclipse_state_->getSimulationConfig(); std::string solver_approach = flowDefaultSolver; if (param_.has("solver_approach")) { @@ -300,7 +269,7 @@ namespace Opm } - Opm::ScheduleConstPtr schedule = eclipseState->getSchedule(); + Opm::ScheduleConstPtr schedule = eclipse_state_->getSchedule(); Opm::TimeMapConstPtr timeMap(schedule->getTimeMap()); SimulatorTimer simtimer; @@ -308,8 +277,8 @@ namespace Opm simtimer.init(timeMap); std::map, double> maxDp; - computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]); - std::vector threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp); + computeMaxDp(maxDp, deck_, eclipse_state_, grid, state, props, gravity[2]); + std::vector threshold_pressures = thresholdPressures(deck_, eclipse_state_, grid, maxDp); Simulator simulator(param_, grid, @@ -318,9 +287,9 @@ namespace Opm rock_comp.isActive() ? &rock_comp : 0, *fis_solver, grav, - deck->hasKeyword("DISGAS"), - deck->hasKeyword("VAPOIL"), - eclipseState, + deck_->hasKeyword("DISGAS"), + deck_->hasKeyword("VAPOIL"), + eclipse_state_, outputWriter, threshold_pressures); @@ -368,13 +337,17 @@ namespace Opm - + // setupParallelism() bool output_cout_ = false; bool must_distribute_ = false; + // setupParameters() parameter::ParameterGroup param_; + // setupOutput() bool output_to_files_ = false; std::string output_dir_ = "output"; - + // readDeckInput() + std::shared_ptr deck_; + std::shared_ptr eclipse_state_; @@ -513,6 +486,49 @@ namespace Opm } + + + + // 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 = std::make_shared(logFile , Log::DefaultMessageTypes); + std::shared_ptr counterLog = std::make_shared(Log::DefaultMessageTypes); + + OpmLog::addBackend( "STREAM" , streamLog ); + OpmLog::addBackend( "COUNTER" , counterLog ); + } + + // Create Deck and EclipseState. + try { + std::string deck_filename = param_.get("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("output_interval"); + IOConfigPtr ioConfig = eclipse_state_->getIOConfig(); + ioConfig->overrideRestartWriteInterval(static_cast(output_interval)); + } + } + }; // class FlowMain From 4711fd3dfe11844881f31ec4b9cac7e1bc2e6e4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 11 Dec 2015 17:01:08 +0100 Subject: [PATCH 09/19] Complete refactoring of FlowMain. --- opm/autodiff/FlowMain.hpp | 493 +++++++++++++++++++++++--------------- 1 file changed, 303 insertions(+), 190 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index de65698c1..d96dafa93 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -122,8 +122,11 @@ namespace Opm class FlowMain { public: + + int execute(int argc, char** argv) try { + // Setup. setupParallelism(argc, argv); printStartupMessage(); const bool ok = setupParameters(argc, argv); @@ -132,196 +135,15 @@ namespace Opm } setupOutput(); readDeckInput(); + setupGridAndProps(); + setupState(); + distributeData(); + setupOutputWriter(); + setupLinearSolver(); + createSimulator(); - std::vector porv = eclipse_state_->getDoubleGridProperty("PORV")->getData(); - GridInit grid_init(deck_, eclipse_state_, porv); - auto&& grid = grid_init.grid(); - - const PhaseUsage pu = Opm::phaseUsageFromDeck(deck_); - - std::vector compressedToCartesianIdx; - Opm::createGlobalCellArray(grid, compressedToCartesianIdx); - - typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager; - auto materialLawManager = std::make_shared(); - materialLawManager->initFromDeck(deck_, eclipse_state_, compressedToCartesianIdx); - - // Rock and fluid init - BlackoilPropertiesFromDeck props( deck_, eclipse_state_, materialLawManager, - Opm::UgGridHelpers::numCells(grid), - Opm::UgGridHelpers::globalCell(grid), - Opm::UgGridHelpers::cartDims(grid), - param_); - - BlackoilPropsAdFromDeck new_props( deck_, eclipse_state_, 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_, eclipse_state_); - - // Gravity. - double gravity[3] = { 0.0 }; - 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_, eclipse_state_, 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 cells(numCells); - for (int c = 0; c < numCells; ++c) { cells[c] = c; } - std::vector 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, eclipse_state_, 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( must_distribute_ ) - { - Opm::distributeGridAndData( grid, deck_, eclipse_state_, 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_, eclipse_state_, pu, new_props.permeability() ); - - // Solver for Newton iterations. - std::unique_ptr 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 simCfg = eclipse_state_->getSimulationConfig(); - std::string solver_approach = flowDefaultSolver; - - if (param_.has("solver_approach")) { - solver_approach = param_.get("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 = eclipse_state_->getSchedule(); - Opm::TimeMapConstPtr timeMap(schedule->getTimeMap()); - SimulatorTimer simtimer; - - // initialize variables - simtimer.init(timeMap); - - std::map, double> maxDp; - computeMaxDp(maxDp, deck_, eclipse_state_, grid, state, props, gravity[2]); - std::vector threshold_pressures = thresholdPressures(deck_, eclipse_state_, grid, maxDp); - - Simulator simulator(param_, - grid, - geoprops, - new_props, - rock_comp.isActive() ? &rock_comp : 0, - *fis_solver, - grav, - deck_->hasKeyword("DISGAS"), - deck_->hasKeyword("VAPOIL"), - eclipse_state_, - 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_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); - warnIfUnusedParams(param_); - } - } else { - outputWriter.writeInit( simtimer ); - if ( output_cout_ ) - { - std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush; - } - } - return EXIT_SUCCESS; + // Run. + return runSimulator(); } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; @@ -332,6 +154,10 @@ namespace Opm private: + typedef BlackoilPropsAdFromDeck FluidProps; + typedef FluidProps::MaterialLawManager MaterialLawManager; + typedef typename Simulator::ReservoirState ReservoirState; + // ------------ Data members ------------ @@ -348,7 +174,25 @@ namespace Opm // readDeckInput() std::shared_ptr deck_; std::shared_ptr eclipse_state_; - + // setupGridAndProps() + std::unique_ptr> grid_init_; + std::shared_ptr material_law_manager_; + std::unique_ptr fluidprops_; + std::unique_ptr rock_comp_; + std::array gravity_; + bool use_local_perm_ = true; + std::unique_ptr geoprops_; + // setupState() + ReservoirState state_; + std::vector threshold_pressures_; + // distributeData() + boost::any parallel_information_; + // setupOutputWriter() + std::unique_ptr output_writer_; + // setupLinearSolver + std::unique_ptr fis_solver_; + // createSimulator() + std::unique_ptr simulator_; // ------------ Methods ------------ @@ -529,6 +373,275 @@ namespace Opm } } + + + + + // 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& porv = eclipse_state_->getDoubleGridProperty("PORV")->getData(); + grid_init_.reset(new GridInit(deck_, eclipse_state_, porv)); + const Grid& grid = grid_init_->grid(); + + // Create material law manager. + std::vector 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, double> maxDp; + computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), state_, props, gravity_[2]); + threshold_pressures_ = thresholdPressures(deck_, eclipse_state_, grid, maxDp); + + // 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 cells(numCells); + for (int c = 0; c < numCells; ++c) { cells[c] = c; } + std::vector 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 simCfg = eclipse_state_->getSimulationConfig(); + std::string solver_approach = flowDefaultSolver; + + if (param_.has("solver_approach")) { + solver_approach = param_.get("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."); + } + } + + + + + + // Create simulator instance. + // Writes to: + // simulator_ + void createSimulator() + { + // Create the simulator instance. + simulator_.reset(new Simulator(param_, + grid_init_->grid(), + *geoprops_, + *fluidprops_, + rock_comp_->isActive() ? rock_comp_.get() : nullptr, + *fis_solver_, + gravity_.data(), + deck_->hasKeyword("DISGAS"), + deck_->hasKeyword("VAPOIL"), + eclipse_state_, + *output_writer_, + threshold_pressures_)); + } + + + + + + // 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 (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); + warnIfUnusedParams(param_); + } + } else { + output_writer_->writeInit( simtimer ); + if (output_cout_) { + std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush; + } + } + return EXIT_SUCCESS; + } + + }; // class FlowMain From 4d31a99f5b9d086e797b60b5128abe2727cea6a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Dec 2015 10:21:21 +0100 Subject: [PATCH 10/19] Moved content of warnIfUnusedParams() into runSimulator(). --- opm/autodiff/FlowMain.hpp | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index d96dafa93..414d4791e 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -100,19 +100,6 @@ 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 @@ -624,13 +611,19 @@ namespace Opm 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); - warnIfUnusedParams(param_); } } else { output_writer_->writeInit( simtimer ); From bb49ddbfc45acad0d5f4ff49b9ed584bf20dca15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Dec 2015 10:46:00 +0100 Subject: [PATCH 11/19] Improve comments. --- opm/autodiff/FlowMain.hpp | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 414d4791e..e5597c798 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -100,17 +100,19 @@ namespace Opm { - /// 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. + /// This class encapsulates the setup and running of + /// a simulator based on an input deck. template class FlowMain { 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. @@ -141,14 +143,23 @@ namespace Opm private: + + + + + // ------------ 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; @@ -185,9 +196,6 @@ namespace Opm // ------------ Methods ------------ - - - // Set up MPI and OpenMP. // Writes to: // output_cout_ From 39900761a8a7839893b1ffcea4d4bc12e418adf0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Dec 2015 13:51:59 +0100 Subject: [PATCH 12/19] Transform FlowMain to use CRTP. Main content is now in FlowMainBase, which takes Implementation as a template parameter. FlowMain inherits from FlowMainBase with itself as Implementation template parameter. Only the createSimulator() method is implemented in FlowMain (as opposed to in FlowMainBase). All subclasses must implement createSimulator(). --- opm/autodiff/FlowMain.hpp | 102 +++++++++++++++++++++++--------------- 1 file changed, 62 insertions(+), 40 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index e5597c798..d8e775875 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -102,8 +102,8 @@ namespace Opm /// This class encapsulates the setup and running of /// a simulator based on an input deck. - template - class FlowMain + template + class FlowMainBase { public: @@ -116,23 +116,23 @@ namespace Opm int execute(int argc, char** argv) try { // Setup. - setupParallelism(argc, argv); - printStartupMessage(); - const bool ok = setupParameters(argc, argv); + asImpl().setupParallelism(argc, argv); + asImpl().printStartupMessage(); + const bool ok = asImpl().setupParameters(argc, argv); if (!ok) { return EXIT_FAILURE; } - setupOutput(); - readDeckInput(); - setupGridAndProps(); - setupState(); - distributeData(); - setupOutputWriter(); - setupLinearSolver(); - createSimulator(); + asImpl().setupOutput(); + asImpl().readDeckInput(); + asImpl().setupGridAndProps(); + asImpl().setupState(); + asImpl().distributeData(); + asImpl().setupOutputWriter(); + asImpl().setupLinearSolver(); + asImpl().createSimulator(); // Run. - return runSimulator(); + return asImpl().runSimulator(); } catch (const std::exception &e) { std::cerr << "Program threw an exception: " << e.what() << "\n"; @@ -141,7 +141,7 @@ namespace Opm - private: + protected: @@ -571,30 +571,6 @@ namespace Opm - // Create simulator instance. - // Writes to: - // simulator_ - void createSimulator() - { - // Create the simulator instance. - simulator_.reset(new Simulator(param_, - grid_init_->grid(), - *geoprops_, - *fluidprops_, - rock_comp_->isActive() ? rock_comp_.get() : nullptr, - *fis_solver_, - gravity_.data(), - deck_->hasKeyword("DISGAS"), - deck_->hasKeyword("VAPOIL"), - eclipse_state_, - *output_writer_, - threshold_pressures_)); - } - - - - - // Run the simulator. // Returns EXIT_SUCCESS if it does not throw. int runSimulator() @@ -643,7 +619,53 @@ namespace Opm } - }; // class FlowMain + + + + // Access the most-derived class used for + // static polymorphism (CRTP). + Implementation& asImpl() + { + return static_cast(*this); + } + + + }; // class FlowMainBase + + + + + + + // The FlowMain class is the basic black-oil simulator case. + template + class FlowMain : public FlowMainBase, Grid, Simulator> + { + protected: + using Base = FlowMainBase, 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 From b156ce0b556026d90040a9c6092f4367ab1194a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Dec 2015 14:04:20 +0100 Subject: [PATCH 13/19] Create and use FlowMainSolvent class. --- CMakeLists_files.cmake | 1 + examples/flow_solvent.cpp | 406 +------------------------------ opm/autodiff/FlowMainSolvent.hpp | 117 +++++++++ 3 files changed, 123 insertions(+), 401 deletions(-) create mode 100644 opm/autodiff/FlowMainSolvent.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 543bf187e..8eb52c67b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -142,6 +142,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/DuneMatrix.hpp opm/autodiff/ExtractParallelGridInformationToISTL.hpp opm/autodiff/FlowMain.hpp + opm/autodiff/FlowMainSolvent.hpp opm/autodiff/GeoProps.hpp opm/autodiff/GridHelpers.hpp opm/autodiff/GridInit.hpp diff --git a/examples/flow_solvent.cpp b/examples/flow_solvent.cpp index 1827bc90b..c7459320c 100644 --- a/examples/flow_solvent.cpp +++ b/examples/flow_solvent.cpp @@ -25,414 +25,18 @@ #endif // HAVE_CONFIG_H -#include - -#include - -#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 3) -#include -#else -#include -#endif - -#if HAVE_DUNE_CORNERPOINT && WANT_DUNE_CORNERPOINTGRID -#define USE_DUNE_CORNERPOINTGRID 1 -#include -#include -#else -#undef USE_DUNE_CORNERPOINTGRID -#endif - -#include - -#include - #include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix. - -#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 - +#include // ----------------- 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 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=, 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("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 streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); - std::shared_ptr counterLog = std::make_shared(Opm::Log::DefaultMessageTypes); - - Opm::OpmLog::addBackend( "STREAM" , streamLog ); - Opm::OpmLog::addBackend( "COUNTER" , counterLog ); - } - - Opm::DeckConstPtr deck; - std::shared_ptr eclipseState; - 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 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("output_interval"); - IOConfigPtr ioConfig = eclipseState->getIOConfig(); - ioConfig->overrideRestartWriteInterval((size_t)output_interval); - } - - const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); - - std::vector compressedToCartesianIdx; - Opm::createGlobalCellArray(grid, compressedToCartesianIdx); - - typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager; - auto materialLawManager = std::make_shared(); - 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 cells(numCells); - for (int c = 0; c < numCells; ++c) { cells[c] = c; } - std::vector 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 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, double> maxDp; - computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]); - std::vector threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp); - std::vector 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 mainfunc; + return mainfunc.execute(argc, argv); } -catch (const std::exception &e) { - std::cerr << "Program threw an exception: " << e.what() << "\n"; - return EXIT_FAILURE; -} - diff --git a/opm/autodiff/FlowMainSolvent.hpp b/opm/autodiff/FlowMainSolvent.hpp new file mode 100644 index 000000000..d3307f31e --- /dev/null +++ b/opm/autodiff/FlowMainSolvent.hpp @@ -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 . +*/ + +#ifndef OPM_FLOWMAINSOLVENT_HEADER_INCLUDED +#define OPM_FLOWMAINSOLVENT_HEADER_INCLUDED + + + +#include +#include + + + +namespace Opm +{ + + // The FlowMainSolvent class is for a black-oil simulator with solvent. + template + class FlowMainSolvent : public FlowMainBase, Grid, Simulator> + { + protected: + using Base = FlowMainBase, Grid, Simulator>; + friend Base; + + // Set in setupGridAndProps() + std::unique_ptr 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 From b420cccfc33bdc45986204258e13169d54fd388e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Dec 2015 14:37:53 +0100 Subject: [PATCH 14/19] Add warning about simulating without WPOLYMER. Warning was in the main() of flow_polymer.cpp, but will disappear from there after refactoring. --- .../fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index b6d78b74f..050d79e38 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -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)); From bfcbd094884e26bbf1b095823083b094049cb39a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 14 Dec 2015 14:42:26 +0100 Subject: [PATCH 15/19] Create and use FlowMainPolymer class. --- CMakeLists_files.cmake | 1 + examples/flow_polymer.cpp | 320 ++----------------------------- opm/autodiff/FlowMainPolymer.hpp | 124 ++++++++++++ 3 files changed, 136 insertions(+), 309 deletions(-) create mode 100644 opm/autodiff/FlowMainPolymer.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 8eb52c67b..3818ab8c3 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -142,6 +142,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/DuneMatrix.hpp opm/autodiff/ExtractParallelGridInformationToISTL.hpp opm/autodiff/FlowMain.hpp + opm/autodiff/FlowMainPolymer.hpp opm/autodiff/FlowMainSolvent.hpp opm/autodiff/GeoProps.hpp opm/autodiff/GridHelpers.hpp diff --git a/examples/flow_polymer.cpp b/examples/flow_polymer.cpp index eac231e1c..9506c0e76 100644 --- a/examples/flow_polymer.cpp +++ b/examples/flow_polymer.cpp @@ -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 . */ -#include "config.h" -#include + +#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 -#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 - +#include // ----------------- Main program ----------------- int main(int argc, char** argv) -try { - using namespace Opm; + typedef UnstructuredGrid Grid; + typedef Opm::SimulatorFullyImplicitBlackoilPolymer 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=, 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 grid; - std::shared_ptr props; - std::shared_ptr new_props; - std::shared_ptr rock_comp; - PolymerBlackoilState state; - // bool check_well_controls = false; - // int max_well_control_iterations = 0; - double gravity[3] = { 0.0 }; - std::string deck_filename = param.get("deck_filename"); - - // Write parameters used for later reference. - bool output = param.getDefault("output", true); - std::string output_dir; - if (output) { - // Create output directory if needed. - output_dir = - param.getDefault("output_dir", std::string("output")); - boost::filesystem::path fpath(output_dir); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - // Write simulation parameters. - param.writeParam(output_dir + "/simulation.param"); - } - - std::string logFile = output_dir + "/LOGFILE.txt"; - Opm::ParserPtr parser(new Opm::Parser()); - { - std::shared_ptr streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); - std::shared_ptr counterLog = std::make_shared(Opm::Log::DefaultMessageTypes); - - Opm::OpmLog::addBackend( "STREAM" , streamLog ); - Opm::OpmLog::addBackend( "COUNTER" , counterLog ); - } - - - Opm::DeckConstPtr deck; - std::shared_ptr eclipseState; - 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 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 compressedToCartesianIdx; - Opm::createGlobalCellArray(*grid->c_grid(), compressedToCartesianIdx); - - typedef BlackoilPropsAdFromDeck::MaterialLawManager MaterialLawManager; - auto materialLawManager = std::make_shared(); - 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 cells(nc); - for (int c = 0; c < nc; ++c) { cells[c] = c; } - std::vector pc = state.saturation(); - props->capPress(nc, state.saturation().data(), cells.data(), pc.data(),NULL); - new_props->setSwatInitScaling(state.saturation(),pc); - } - - bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); - const double *grav = use_gravity ? &gravity[0] : 0; - - // Solver for Newton iterations. - std::unique_ptr fis_solver; - if (param.getDefault("use_cpr", true)) { - fis_solver.reset(new NewtonIterationBlackoilCPR(param)); - } else { - fis_solver.reset(new NewtonIterationBlackoilSimple(param)); - } - - Opm::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, double> maxDp; - computeMaxDp(maxDp, deck, eclipseState, *grid->c_grid(), state, *props, gravity[2]); - std::vector threshold_pressures = thresholdPressures(deck, eclipseState, *grid->c_grid(), maxDp); - std::vector 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 - 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 mainfunc; + return mainfunc.execute(argc, argv); } diff --git a/opm/autodiff/FlowMainPolymer.hpp b/opm/autodiff/FlowMainPolymer.hpp new file mode 100644 index 000000000..a917eb8e6 --- /dev/null +++ b/opm/autodiff/FlowMainPolymer.hpp @@ -0,0 +1,124 @@ +/* + 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 . +*/ + +#ifndef OPM_FLOWMAINPOLYMER_HEADER_INCLUDED +#define OPM_FLOWMAINPOLYMER_HEADER_INCLUDED + + + +#include +#include +#include + + + +namespace Opm +{ + + // The FlowMainPolymer class is for a black-oil simulator with polymer. + template + class FlowMainPolymer : public FlowMainBase, Grid, Simulator> + { + protected: + using Base = FlowMainBase, Grid, Simulator>; + friend Base; + + // Set in setupGridAndProps() + std::unique_ptr polymer_props_legacy_; // Held by reference in polymer_props_ + std::unique_ptr 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_)); + } + + + + + + // 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 From 78dbb79ea4d8ec0f8be42135a8eeac34786aa357 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Dec 2015 11:12:01 +0100 Subject: [PATCH 16/19] Fix compile error for in-class member init. Reported by Kai Bao. --- opm/autodiff/FlowMain.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index d8e775875..6f84e2c05 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -168,7 +168,7 @@ namespace Opm parameter::ParameterGroup param_; // setupOutput() bool output_to_files_ = false; - std::string output_dir_ = "output"; + std::string output_dir_ = std::string("output"); // readDeckInput() std::shared_ptr deck_; std::shared_ptr eclipse_state_; From dc4274f4a2eba25eca2bf4ba1316cfed0ae27ea3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Dec 2015 11:12:55 +0100 Subject: [PATCH 17/19] Include nncs in threshold pressures. This was left out when rebasing (it had been introduced in the flow etc. mains). This also fixes a bug: the nnc thresholds should be inserted at the end of the vector to be consistent with the operators, not at the beginning. --- opm/autodiff/FlowMain.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 6f84e2c05..2844cecef 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -477,6 +477,8 @@ namespace Opm std::map, double> maxDp; computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), state_, props, gravity_[2]); threshold_pressures_ = thresholdPressures(deck_, eclipse_state_, grid, maxDp); + std::vector 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")) { From 8507427d83daa46777acb78fe757dfa704eb3a02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Dec 2015 14:55:33 +0100 Subject: [PATCH 18/19] Use correct syntax for logical not. --- opm/autodiff/GeoProps.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/autodiff/GeoProps.hpp b/opm/autodiff/GeoProps.hpp index a154ba708..e375c6bcf 100644 --- a/opm/autodiff/GeoProps.hpp +++ b/opm/autodiff/GeoProps.hpp @@ -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(); From 6839b811f2a278b11905484970c6ef588fa52d0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Dec 2015 15:05:32 +0100 Subject: [PATCH 19/19] Make polymer solver use direct solver only. The interleaved solver should work for this case, but does not currently. --- opm/autodiff/FlowMainPolymer.hpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/opm/autodiff/FlowMainPolymer.hpp b/opm/autodiff/FlowMainPolymer.hpp index a917eb8e6..b594b3c0b 100644 --- a/opm/autodiff/FlowMainPolymer.hpp +++ b/opm/autodiff/FlowMainPolymer.hpp @@ -87,6 +87,19 @@ namespace Opm + // 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_