diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 9bf20a207..6ee8469b2 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -84,6 +84,10 @@ namespace Properties { NEW_TYPE_TAG(EclFlowProblem, INHERITS_FROM(BlackOilModel, EclBaseProblem)); SET_BOOL_PROP(EclFlowProblem, DisableWells, true); SET_BOOL_PROP(EclFlowProblem, EnableDebuggingChecks, false); + +// SWATINIT is done by the flow part of flow_ebos. this can be removed once the legacy +// code for fluid and satfunc handling gets fully retired. +SET_BOOL_PROP(EclFlowProblem, EnableSwatinit, false); }} namespace Opm { @@ -144,7 +148,8 @@ namespace Opm { FIP_PV = 5, //< Pore volume FIP_WEIGHTED_PRESSURE = 6 }; - std::array, 7> fip; + static const int fipValues = FIP_WEIGHTED_PRESSURE + 1 ; + std::array, fipValues> fip; }; // --------- Public methods --------- @@ -156,7 +161,6 @@ namespace Opm { /// \param[in] grid grid data structure /// \param[in] fluid fluid properties /// \param[in] geo rock properties - /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] wells well structure /// \param[in] vfp_properties Vertical flow performance tables /// \param[in] linsolver linear solver @@ -166,7 +170,6 @@ namespace Opm { const ModelParameters& param, const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , - const RockCompressibility* rock_comp_props, const StandardWellsDense& well_model, const NewtonIterationBlackoilInterface& linsolver, const bool terminal_output) @@ -189,7 +192,6 @@ namespace Opm { , isBeginReportStep_(false) , invalidateIntensiveQuantitiesCache_(true) { - DUNE_UNUSED_PARAMETER(rock_comp_props); const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); const std::vector pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); const std::vector depth(geo_.z().data(), geo_.z().data() + geo_.z().size()); @@ -226,7 +228,7 @@ namespace Opm { const EclipseState& eclState() const - { return *ebosSimulator_.gridManager().eclState(); } + { return ebosSimulator_.gridManager().eclState(); } /// Called once before each time step. /// \param[in] timer simulation timer @@ -300,7 +302,7 @@ namespace Opm { // Compute the nonlinear update. const int nc = AutoDiffGrid::numCells(grid_); - const int nw = wellModel().wells().number_of_wells; + const int nw = numWells(); BVector x(nc); BVector xw(nw); @@ -447,7 +449,7 @@ namespace Opm { int sizeNonLinear() const { const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int nw = wellModel().wells().number_of_wells; + const int nw = numWells(); return numPhases() * (nc + nw); } @@ -476,8 +478,11 @@ namespace Opm { const auto& ebosJac = ebosSimulator_.model().linearizer().matrix(); auto& ebosResid = ebosSimulator_.model().linearizer().residual(); - // apply well residual to the residual. - wellModel().apply(ebosResid); + if( xw.size() > 0 ) + { + // apply well residual to the residual. + wellModel().apply(ebosResid); + } // set initial guess x = 0.0; @@ -497,9 +502,12 @@ namespace Opm { istlSolver().solve( opA, x, ebosResid ); } - // recover wells. - xw = 0.0; - wellModel().recoverVariable(x, xw); + if( xw.size() > 0 ) + { + // recover wells. + xw = 0.0; + wellModel().recoverVariable(x, xw); + } } //===================================================================== @@ -757,7 +765,6 @@ namespace Opm { } - /// Return true if output to cout is wanted. bool terminalOutputEnabled() const { @@ -1015,7 +1022,7 @@ namespace Opm { const auto& pv = geo_.poreVolume(); const int maxnp = Opm::BlackoilPhases::MaxNumPhases; - for (int i = 0; i<7; i++) { + for (int i = 0; i> values(dims, std::vector(7,0.0)); + std::vector> values(dims, std::vector(FIPData::fipValues,0.0)); std::vector hcpv(dims, 0.0); std::vector pres(dims, 0.0); @@ -1107,11 +1114,12 @@ namespace Opm { // mask[c] is 1 if we need to compute something in parallel const auto & pinfo = boost::any_cast(istlSolver().parallelInformation()); - const auto& mask = pinfo.getOwnerMask(); + const auto& mask = pinfo.updateOwnerMask( fipnum ); + auto comm = pinfo.communicator(); // Compute the global dims value and resize values accordingly. dims = comm.max(dims); - values.resize(dims, std::vector(7,0.0)); + values.resize(dims, std::vector(FIPData::fipValues,0.0)); //Accumulate phases for each region for (int phase = 0; phase < maxnp; ++phase) { @@ -1249,6 +1257,8 @@ namespace Opm { /// return true if wells are available in the reservoir bool wellsActive() const { return well_model_.wellsActive(); } + int numWells() const { return wellsActive() ? wells().number_of_wells : 0; } + /// return true if wells are available on this process bool localWellsActive() const { return well_model_.localWellsActive(); } diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 61bd1f710..2e3d07d7b 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -152,7 +152,7 @@ namespace Opm asImpl().setupOutputWriter(); asImpl().setupLinearSolver(); asImpl().createSimulator(); - + // Run. auto ret = asImpl().runSimulator(); @@ -388,12 +388,12 @@ namespace Opm - // Setup OpmLog backend with output_dir. + // Setup OpmLog backend with output_dir. void setupLogging() { std::string deck_filename = param_.get("deck_filename"); // create logFile - using boost::filesystem::path; + using boost::filesystem::path; path fpath(deck_filename); std::string baseName; std::ostringstream debugFileStream; diff --git a/opm/autodiff/FlowMainEbos.hpp b/opm/autodiff/FlowMainEbos.hpp index f41fd098a..13eb8eaf2 100644 --- a/opm/autodiff/FlowMainEbos.hpp +++ b/opm/autodiff/FlowMainEbos.hpp @@ -23,23 +23,137 @@ #ifndef OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED #define OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED -#include +#include + #include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include #include namespace Opm { // The FlowMain class is the ebos based black-oil simulator. - class FlowMainEbos : public FlowMainBase + class FlowMainEbos { - protected: - typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator; - typedef FlowMainBase Base; - friend Base; - - typedef typename TTAG(EclFlowProblem) TypeTag; + public: + typedef TTAG(EclFlowProblem) TypeTag; + typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager MaterialLawManager; typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator; + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator; + typedef typename Simulator::ReservoirState ReservoirState; + typedef typename Simulator::OutputWriter OutputWriter; + + /// 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 { + setupParallelism(argc, argv); + printStartupMessage(); + const bool ok = setupParameters(argc, argv); + if (!ok) { + return EXIT_FAILURE; + } + + setupOutput(); + setupEbosSimulator(); + setupLogging(); + extractMessages(); + setupGridAndProps(); + runDiagnostics(); + setupState(); + writeInit(); + setupOutputWriter(); + setupLinearSolver(); + createSimulator(); + + // Run. + auto ret = runSimulator(); + + mergeParallelLogFiles(); + + return ret; + } + catch (const std::exception &e) { + std::ostringstream message; + message << "Program threw an exception: " << e.what(); + + if( output_cout_ ) + { + // in some cases exceptions are thrown before the logging system is set + // up. + if (OpmLog::hasBackend("STREAMLOG")) { + OpmLog::error(message.str()); + } + else { + std::cout << message.str() << "\n"; + } + } + + return EXIT_FAILURE; + } + } + + protected: + 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); + mpi_rank_ = mpi_helper.rank(); + const int mpi_size = mpi_helper.size(); + output_cout_ = ( mpi_rank_ == 0 ); + must_distribute_ = ( mpi_size > 1 ); + +#ifdef _OPENMP + // OpenMP setup. + if (!getenv("OMP_NUM_THREADS")) { + // Default to at most 4 threads, regardless of + // number of cores (unless ENV(OMP_NUM_THREADS) is defined) + int num_cores = omp_get_num_procs(); + int num_threads = std::min(4, num_cores); + omp_set_num_threads(num_threads); + } +#pragma omp parallel + if (omp_get_thread_num() == 0) { + // omp_get_num_threads() only works as expected within a parallel region. + const int num_omp_threads = omp_get_num_threads(); + if (mpi_size == 1) { + std::cout << "OpenMP using " << num_omp_threads << " threads." << std::endl; + } else { + std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank_ << "." << std::endl; + } + } +#endif + } // Print startup message if on output rank. void printStartupMessage() @@ -66,12 +180,157 @@ namespace Opm } } - // Parser the input and creates the Deck and EclipseState objects. + // Read parameters, see if a deck was specified on the command line, and if + // it was, insert it into parameters. // Writes to: - // deck_ - // eclipse_state_ - // May throw if errors are encountered, here configured to be somewhat tolerant. - void readDeckInput() + // param_ + // Returns true if ok, false if not. + bool setupParameters(int argc, char** argv) + { + 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 { + const auto casename = this->simulationCaseName( param_.unhandledArguments()[ 0 ] ); + param_.insertParameter("deck_filename", casename.string() ); + } + } + + // 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; + } + + // 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); + // Always read output_dir as it will be set unconditionally later. + // Not doing this might cause files to be created in the current + // directory. + output_dir_ = + param_.getDefault("output_dir", std::string(".")); + + if (output_to_files_) { + // Create output directory if needed. + boost::filesystem::path fpath(output_dir_); + if (!is_directory(fpath)) { + try { + create_directories(fpath); + } + catch (...) { + OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); + } + } + // Write simulation parameters. + param_.writeParam(output_dir_ + "/simulation.param"); + } + } + + // Setup OpmLog backend with output_dir. + void setupLogging() + { + std::string deck_filename = param_.get("deck_filename"); + // create logFile + using boost::filesystem::path; + path fpath(deck_filename); + std::string baseName; + std::ostringstream debugFileStream; + std::ostringstream logFileStream; + + if (boost::to_upper_copy(path(fpath.extension()).string()) == ".DATA") { + baseName = path(fpath.stem()).string(); + } else { + baseName = path(fpath.filename()).string(); + } + if (param_.has("output_dir")) { + logFileStream << output_dir_ << "/"; + debugFileStream << output_dir_ + "/"; + } + + logFileStream << baseName; + debugFileStream << "." << baseName; + + if ( must_distribute_ && mpi_rank_ != 0 ) + { + // Added rank to log file for non-zero ranks. + // This prevents message loss. + debugFileStream << "."<< mpi_rank_; + // If the following file appears then there is a bug. + logFileStream << "." << mpi_rank_; + } + logFileStream << ".PRT"; + debugFileStream << ".DEBUG"; + + std::string debugFile = debugFileStream.str(); + logFile_ = logFileStream.str(); + + std::shared_ptr prtLog = std::make_shared(logFile_ , Log::NoDebugMessageTypes, false, output_cout_); + std::shared_ptr streamLog = std::make_shared(std::cout, Log::StdoutMessageTypes); + OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog ); + OpmLog::addBackend( "STREAMLOG", streamLog); + std::shared_ptr debugLog = std::make_shared(debugFile, Log::DefaultMessageTypes, false, output_cout_); + OpmLog::addBackend( "DEBUGLOG" , debugLog); + const auto& msgLimits = eclState().getSchedule().getMessageLimits(); + const std::map limits = {{Log::MessageType::Note, msgLimits.getCommentPrintLimit(0)}, + {Log::MessageType::Info, msgLimits.getMessagePrintLimit(0)}, + {Log::MessageType::Warning, msgLimits.getWarningPrintLimit(0)}, + {Log::MessageType::Error, msgLimits.getErrorPrintLimit(0)}, + {Log::MessageType::Problem, msgLimits.getProblemPrintLimit(0)}, + {Log::MessageType::Bug, msgLimits.getBugPrintLimit(0)}}; + prtLog->setMessageLimiter(std::make_shared()); + prtLog->setMessageFormatter(std::make_shared(false)); + streamLog->setMessageLimiter(std::make_shared(10, limits)); + streamLog->setMessageFormatter(std::make_shared(true)); + + // Read parameters. + if ( output_cout_ ) + { + OpmLog::debug("\n--------------- Reading parameters ---------------\n"); + } + } + + void mergeParallelLogFiles() + { + // force closing of all log files. + OpmLog::removeAllBackends(); + + if( mpi_rank_ != 0 || !must_distribute_ ) + { + return; + } + + namespace fs = boost::filesystem; + fs::path output_path("."); + if ( param_.has("output_dir") ) + { + output_path = fs::path(output_dir_); + } + + fs::path deck_filename(param_.get("deck_filename")); + + std::for_each(fs::directory_iterator(output_path), + fs::directory_iterator(), + detail::ParallelFileMerger(output_path, deck_filename.stem().string())); + } + + void setupEbosSimulator() { std::string progName("flow_ebos"); std::string deckFile("--ecl-deck-file-name="); @@ -84,23 +343,286 @@ namespace Opm ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false)); ebosSimulator_->model().applyInitialSolution(); - Base::deck_ = ebosSimulator_->gridManager().deck(); - Base::eclipse_state_ = ebosSimulator_->gridManager().eclState(); - IOConfig& ioConfig = Base::eclipse_state_->getIOConfig(); - ioConfig.setOutputDir(Base::output_dir_); + try { + if (output_cout_) { + MissingFeatures::checkKeywords(deck()); + } + + IOConfig& ioConfig = eclState().getIOConfig(); + ioConfig.setOutputDir(output_dir_); + + // Possible to force initialization only behavior (NOSIM). + if (param_.has("nosim")) { + const bool nosim = param_.get("nosim"); + ioConfig.overrideNOSIM( nosim ); + } + } + 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 (Base::param_.has("output_interval")) { - const int output_interval = Base::param_.get("output_interval"); - eclipse_state_->getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) ); + if (param_.has("output_interval")) { + const int output_interval = param_.get("output_interval"); + eclState().getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) ); + } + } + // Create grid and property objects. + // Writes to: + // material_law_manager_ + // fluidprops_ + // gravity_ + void setupGridAndProps() + { + Dune::CpGrid& grid = ebosSimulator_->gridManager().grid(); + material_law_manager_ = ebosSimulator_->problem().materialLawManager(); + + // create the legacy properties objects + fluidprops_.reset(new BlackoilPropsAdFromDeck(deck(), + eclState(), + material_law_manager_, + grid)); + + // Gravity. + static_assert(Grid::dimension == 3, "Only 3D grids are supported"); + gravity_.fill(0.0); + if (!deck().hasKeyword("NOGRAV")) + gravity_[2] = + param_.getDefault("gravity", unit::gravity); + + // Geological properties + bool use_local_perm = param_.getDefault("use_local_perm", true); + geoprops_.reset(new DerivedGeology(grid, *fluidprops_, eclState(), use_local_perm, gravity_.data())); + } + + const Deck& deck() const + { return ebosSimulator_->gridManager().deck(); } + + Deck& deck() + { return ebosSimulator_->gridManager().deck(); } + + const EclipseState& eclState() const + { return ebosSimulator_->gridManager().eclState(); } + + EclipseState& eclState() + { return ebosSimulator_->gridManager().eclState(); } + + // 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 = this->grid(); + + // Need old-style fluid object for init purposes (only). + BlackoilPropertiesFromDeck props(deck(), + eclState(), + material_law_manager_, + grid.size(/*codim=*/0), + grid.globalCell().data(), + grid.logicalCartesianSize().data(), + param_); + + + // Init state variables (saturation and pressure). + if (param_.has("init_saturation")) { + state_.reset(new ReservoirState(grid.size(/*codim=*/0), + grid.numFaces(), + props.numPhases())); + + initStateBasic(grid.size(/*codim=*/0), + grid.globalCell().data(), + grid.logicalCartesianSize().data(), + grid.numFaces(), + Opm::UgGridHelpers::faceCells(grid), + Opm::UgGridHelpers::beginFaceCentroids(grid), + Opm::UgGridHelpers::beginCellCentroids(grid), + Grid::dimension, + 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); + + // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState. + auto& gor = state_->getCellData( BlackoilState::GASOILRATIO ); + const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL ); + for (int c = 0; c < numCells; ++c) { + // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field. + gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]]; + } + } + } else if (deck().hasKeyword("EQUIL")) { + // Which state class are we really using - what a f... mess? + state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases())); + + initStateEquil(grid, props, deck(), eclState(), gravity_[2], *state_); + //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); + } else { + state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases())); + 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_); } - // Possible to force initialization only behavior (NOSIM). - if (Base::param_.has("nosim")) { - const bool nosim = Base::param_.get("nosim"); - ioConfig.overrideNOSIM( nosim ); + // 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); } + initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck().hasKeyword("DISGAS"), deck().hasKeyword("VAPOIL")); + } + + // Extract messages from parser. + // Writes to: + // OpmLog singleton. + void extractMessages() + { + if ( !output_cout_ ) + { + return; + } + + auto extractMessage = [this](const Message& msg) { + auto log_type = this->convertMessageType(msg.mtype); + const auto& location = msg.location; + if (location) { + OpmLog::addMessage(log_type, Log::fileMessage(location.filename, location.lineno, msg.message)); + } else { + OpmLog::addMessage(log_type, msg.message); + } + }; + + // Extract messages from Deck. + for(const auto& msg : deck().getMessageContainer()) { + extractMessage(msg); + } + + // Extract messages from EclipseState. + for (const auto& msg : eclState().getMessageContainer()) { + extractMessage(msg); + } + } + + // Run diagnostics. + // Writes to: + // OpmLog singleton. + void runDiagnostics() + { + if( ! output_cout_ ) + { + return; + } + + // Run relperm diagnostics + RelpermDiagnostics diagnostic; + diagnostic.diagnosis(eclState(), deck(), this->grid()); + } + + void writeInit() + { + bool output = param_.getDefault("output", true); + bool output_ecl = param_.getDefault("output_ecl", true); + const Grid& grid = this->grid(); + if( output && output_ecl && output_cout_) + { + const EclipseGrid& inputGrid = eclState().getInputGrid(); + eclipse_writer_.reset(new EclipseWriter(eclState(), UgGridHelpers::createEclipseGrid( grid , inputGrid ))); + +#warning TODO +#if 0 + eclipse_writer_->writeInitial(geoprops_->simProps(grid), + geoprops_->nonCartesianConnections()); +#endif + } + } + + // 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 OutputWriter(grid(), + param_, + eclState(), + std::move(eclipse_writer_), + Opm::phaseUsageFromDeck(deck()), + fluidprops_->permeability())); + } + + // Run the simulator. + // Returns EXIT_SUCCESS if it does not throw. + int runSimulator() + { + const auto& schedule = eclState().getSchedule(); + const auto& timeMap = schedule.getTimeMap(); + auto& ioConfig = eclState().getIOConfig(); + SimulatorTimer simtimer; + + // initialize variables + const auto& initConfig = eclState().getInitConfig(); + simtimer.init(timeMap, (size_t)initConfig.getRestartStep()); + + if (!ioConfig.initOnly()) { + if (output_cout_) { + std::string msg; + msg = "\n\n================ Starting main simulation loop ===============\n"; + OpmLog::info(msg); + } + + SimulatorReport fullReport = simulator_->run(simtimer, *state_); + + if (output_cout_) { + std::ostringstream ss; + ss << "\n\n================ End of simulation ===============\n\n"; + fullReport.reportFullyImplicit(ss); + OpmLog::info(ss.str()); + if (param_.anyUnused()) { + // This allows a user to catch typos and misunderstandings in the + // use of simulator parameters. + std::cout << "-------------------- Unused parameters: --------------------\n"; + param_.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + } + + if (output_to_files_) { + std::string filename = output_dir_ + "/walltime.txt"; + std::fstream tot_os(filename.c_str(), std::fstream::trunc | std::fstream::out); + fullReport.reportParam(tot_os); + } + } else { + if (output_cout_) { + std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush; + } + + } + return EXIT_SUCCESS; } // Setup linear solver. @@ -109,7 +631,9 @@ namespace Opm void setupLinearSolver() { typedef typename BlackoilModelEbos :: ISTLSolverType ISTLSolverType; - Base::fis_solver_.reset( new ISTLSolverType( Base::param_, Base::parallel_information_ ) ); + + extractParallelGridInformationToISTL(grid(), parallel_information_); + fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) ); } /// This is the main function of Flow. @@ -119,22 +643,93 @@ namespace Opm void createSimulator() { // Create the simulator instance. - Base::simulator_.reset(new Simulator(*ebosSimulator_, - Base::param_, - *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_)); + simulator_.reset(new Simulator(*ebosSimulator_, + param_, + *geoprops_, + *fluidprops_, + *fis_solver_, + gravity_.data(), + FluidSystem::enableDissolvedGas(), + FluidSystem::enableVaporizedOil(), + eclState(), + *output_writer_, + defunctWellNames())); } private: + boost::filesystem::path simulationCaseName( const std::string& casename ) { + namespace fs = boost::filesystem; + + const auto exists = []( const fs::path& f ) -> bool { + if( !fs::exists( f ) ) return false; + + if( fs::is_regular_file( f ) ) return true; + + return fs::is_symlink( f ) + && fs::is_regular_file( fs::read_symlink( f ) ); + }; + + auto simcase = fs::path( casename ); + + if( exists( simcase ) ) { + return simcase; + } + + for( const auto& ext : { std::string("data"), std::string("DATA") } ) { + if( exists( simcase.replace_extension( ext ) ) ) { + return simcase; + } + } + + throw std::invalid_argument( "Cannot find input case " + casename ); + } + + + int64_t convertMessageType(const Message::type& mtype) + { + switch (mtype) { + case Message::type::Debug: + return Log::MessageType::Debug; + case Message::type::Info: + return Log::MessageType::Info; + case Message::type::Warning: + return Log::MessageType::Warning; + case Message::type::Error: + return Log::MessageType::Error; + case Message::type::Problem: + return Log::MessageType::Problem; + case Message::type::Bug: + return Log::MessageType::Bug; + case Message::type::Note: + return Log::MessageType::Note; + } + throw std::logic_error("Invalid messages type!\n"); + } + + Grid& grid() + { return ebosSimulator_->gridManager().grid(); } + + std::unordered_set defunctWellNames() const + { return ebosSimulator_->gridManager().defunctWellNames(); } + std::unique_ptr ebosSimulator_; + int mpi_rank_ = 0; + bool output_cout_ = false; + bool must_distribute_ = false; + parameter::ParameterGroup param_; + bool output_to_files_ = false; + std::string output_dir_ = std::string("."); + std::shared_ptr material_law_manager_; + std::unique_ptr fluidprops_; + std::array gravity_; + std::unique_ptr geoprops_; + std::unique_ptr state_; + std::unique_ptr eclipse_writer_; + std::unique_ptr output_writer_; + boost::any parallel_information_; + std::unique_ptr fis_solver_; + std::unique_ptr simulator_; + std::string logFile_; }; } // namespace Opm diff --git a/opm/autodiff/GridInit.hpp b/opm/autodiff/GridInit.hpp index 5b900ba22..1cbed0a49 100644 --- a/opm/autodiff/GridInit.hpp +++ b/opm/autodiff/GridInit.hpp @@ -73,18 +73,42 @@ namespace Opm class GridInit { public: + GridInit() + { + gridSelfManaged_ = false; + } + /// Initialize from a deck and/or an eclipse state and (logical cartesian) specified pore volumes. GridInit(const EclipseState& eclipse_state, const std::vector& porv) { - grid_.processEclipseFormat(eclipse_state.getInputGrid(), false, false, false, porv); + gridSelfManaged_ = true; + + grid_ = new Dune::CpGrid; + grid_->processEclipseFormat(eclipse_state.getInputGrid(), false, false, false, porv); } + + ~GridInit() + { + if (gridSelfManaged_) + delete grid_; + } + /// Access the created grid. Note that mutable access may be required for load balancing. Dune::CpGrid& grid() { - return grid_; + return *grid_; } + + /// set the grid from the outside + void setGrid(Dune::CpGrid& newGrid) + { + gridSelfManaged_ = false; + grid_ = &newGrid; + } + private: - Dune::CpGrid grid_; + Dune::CpGrid* grid_; + bool gridSelfManaged_; }; #endif // HAVE_OPM_GRID diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp index fcb91801c..32fd09bed 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp @@ -84,7 +84,6 @@ public: /// /// \param[in] geo derived geological properties /// \param[in] props fluid and rock properties - /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] linsolver linear solver /// \param[in] gravity if non-null, gravity vector /// \param[in] has_disgas true for dissolved gas option @@ -96,20 +95,18 @@ public: const parameter::ParameterGroup& param, DerivedGeology& geo, BlackoilPropsAdInterface& props, - const RockCompressibility* rock_comp_props, NewtonIterationBlackoilInterface& linsolver, const double* gravity, const bool has_disgas, const bool has_vapoil, - std::shared_ptr eclipse_state, + const EclipseState& eclState, BlackoilOutputWriterEbos& output_writer, - const std::vector& threshold_pressures_by_face) + const std::unordered_set& defunct_well_names) : ebosSimulator_(ebosSimulator), param_(param), model_param_(param), solver_param_(param), props_(props), - rock_comp_props_(rock_comp_props), gravity_(gravity), geo_(geo), solver_(linsolver), @@ -117,10 +114,10 @@ public: has_vapoil_(has_vapoil), terminal_output_(param.getDefault("output_terminal", true)), output_writer_(output_writer), - threshold_pressures_by_face_(threshold_pressures_by_face), + defunct_well_names_( defunct_well_names ), is_parallel_run_( false ) { - DUNE_UNUSED_PARAMETER(eclipse_state); + DUNE_UNUSED_PARAMETER(eclState); // Misc init. const int num_cells = AutoDiffGrid::numCells(grid()); allcells_.resize(num_cells); @@ -229,7 +226,8 @@ public: props_.permeability(), dynamic_list_econ_limited, is_parallel_run_, - well_potentials ); + well_potentials, + defunct_well_names_ ); const Wells* wells = wells_manager.c_wells(); WellState well_state; @@ -395,7 +393,6 @@ protected: model_param_, props_, geo_, - rock_comp_props_, well_model, solver_, terminal_output_)); @@ -699,10 +696,10 @@ protected: const EclipseState& eclState() const - { return *ebosSimulator_.gridManager().eclState(); } + { return ebosSimulator_.gridManager().eclState(); } EclipseState& eclState() - { return *ebosSimulator_.gridManager().eclState(); } + { return ebosSimulator_.gridManager().eclState(); } // Data. Simulator& ebosSimulator_; @@ -718,7 +715,6 @@ protected: // Observed objects. BlackoilPropsAdInterface& props_; - const RockCompressibility* rock_comp_props_; const double* gravity_; // Solvers DerivedGeology& geo_; @@ -731,8 +727,11 @@ protected: // output_writer OutputWriter& output_writer_; std::unique_ptr rateConverter_; - // Threshold pressures. - std::vector threshold_pressures_by_face_; + // The names of wells that should be defunct + // (e.g. in a parallel run when they are handeled by + // a different process) + std::unordered_set defunct_well_names_; + // Whether this a parallel simulation or not bool is_parallel_run_; diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index 0a7951085..c865713b4 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -92,14 +92,17 @@ enum WellVariablePositions { , fluid_(nullptr) , active_(nullptr) , vfp_properties_(nullptr) - , well_perforation_densities_(wells_arg->well_connpos[wells_arg->number_of_wells]) - , well_perforation_pressure_diffs_(wells_arg->well_connpos[wells_arg->number_of_wells]) - , wellVariables_(wells_arg->number_of_wells * wells_arg->number_of_phases) - , F0_(wells_arg->number_of_wells * wells_arg->number_of_phases) + , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) + , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) + , wellVariables_( wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0) + , F0_(wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0 ) { - invDuneD_.setBuildMode( Mat::row_wise ); - duneC_.setBuildMode( Mat::row_wise ); - duneB_.setBuildMode( Mat::row_wise ); + if( wells_ ) + { + invDuneD_.setBuildMode( Mat::row_wise ); + duneC_.setBuildMode( Mat::row_wise ); + duneB_.setBuildMode( Mat::row_wise ); + } } void init(const BlackoilPropsAdInterface* fluid_arg, @@ -709,6 +712,10 @@ enum WellVariablePositions { std::vector residual() { + if( ! wellsActive() ) + { + return std::vector(); + } const int np = numPhases(); const int nw = wells().number_of_wells;