Merge pull request from andlaus/no_opm-parser_pointers

flow_ebos: do not use (smart) pointers for opm-parser objects anymore
This commit is contained in:
Andreas Lauser 2016-12-08 10:03:24 +01:00 committed by GitHub
commit d9c01406ea
6 changed files with 716 additions and 81 deletions

View File

@ -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<std::vector<double>, 7> fip;
static const int fipValues = FIP_WEIGHTED_PRESSURE + 1 ;
std::array<std::vector<double>, 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<FluidSystem, BlackoilIndices>& 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<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
const std::vector<double> 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<FIPData::fipValues; i++) {
fip_.fip[i].resize(nc,0.0);
}
@ -1039,7 +1046,7 @@ namespace Opm {
// For a parallel run this is just a local maximum and needs to be updated later
int dims = *std::max_element(fipnum.begin(), fipnum.end());
std::vector<std::vector<double>> values(dims, std::vector<double>(7,0.0));
std::vector<std::vector<double>> values(dims, std::vector<double>(FIPData::fipValues,0.0));
std::vector<double> hcpv(dims, 0.0);
std::vector<double> 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<const ParallelISTLInformation&>(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<double>(7,0.0));
values.resize(dims, std::vector<double>(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(); }

View File

@ -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<std::string>("deck_filename");
// create logFile
using boost::filesystem::path;
using boost::filesystem::path;
path fpath(deck_filename);
std::string baseName;
std::ostringstream debugFileStream;

View File

@ -23,23 +23,137 @@
#ifndef OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED
#define OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED
#include <opm/autodiff/FlowMain.hpp>
#include <opm/simulators/ParallelFileMerger.hpp>
#include <opm/autodiff/BlackoilModelEbos.hpp>
#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
#include <opm/autodiff/MissingFeatures.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/moduleVersion.hpp>
#include <opm/autodiff/ExtractParallelGridInformationToISTL.hpp>
#include <opm/core/props/satfunc/RelpermDiagnostics.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/OpmLog/EclipsePRTLog.hpp>
#include <opm/common/OpmLog/LogUtil.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/IOConfig/IOConfig.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
#include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
#include <ewoms/version.hh>
namespace Opm
{
// The FlowMain class is the ebos based black-oil simulator.
class FlowMainEbos : public FlowMainBase<FlowMainEbos, Dune::CpGrid, Opm::SimulatorFullyImplicitBlackoilEbos>
class FlowMainEbos
{
protected:
typedef Opm::SimulatorFullyImplicitBlackoilEbos Simulator;
typedef FlowMainBase<FlowMainEbos, Dune::CpGrid, Simulator> 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=<path to your deck>, or\n"
" c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
return false;
}
return true;
}
// Set output_to_files_ and set/create output dir. Write parameter file.
// Writes to:
// output_to_files_
// output_dir_
// Throws std::runtime_error if failed to create (if requested) output dir.
void setupOutput()
{
// Write parameters used for later reference. (only if rank is zero)
output_to_files_ = output_cout_ && param_.getDefault("output", true);
// 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<std::string>("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<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile_ , Log::NoDebugMessageTypes, false, output_cout_);
std::shared_ptr<StreamLog> streamLog = std::make_shared<StreamLog>(std::cout, Log::StdoutMessageTypes);
OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
OpmLog::addBackend( "STREAMLOG", streamLog);
std::shared_ptr<StreamLog> debugLog = std::make_shared<EclipsePRTLog>(debugFile, Log::DefaultMessageTypes, false, output_cout_);
OpmLog::addBackend( "DEBUGLOG" , debugLog);
const auto& msgLimits = eclState().getSchedule().getMessageLimits();
const std::map<int64_t, int> 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<MessageLimiter>());
prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(false));
streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10, limits));
streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(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<std::string>("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<bool>("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<int>("output_interval");
eclipse_state_->getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) );
if (param_.has("output_interval")) {
const int output_interval = param_.get<int>("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<bool>("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<int> cells(numCells);
for (int c = 0; c < numCells; ++c) { cells[c] = c; }
std::vector<double> pc = state_->saturation();
props.capPress(numCells, state_->saturation().data(), cells.data(), pc.data(), nullptr);
fluidprops_->setSwatInitScaling(state_->saturation(), pc);
}
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<std::string> defunctWellNames() const
{ return ebosSimulator_->gridManager().defunctWellNames(); }
std::unique_ptr<EbosSimulator> 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<MaterialLawManager> material_law_manager_;
std::unique_ptr<BlackoilPropsAdFromDeck> fluidprops_;
std::array<double, 3> gravity_;
std::unique_ptr<DerivedGeology> geoprops_;
std::unique_ptr<ReservoirState> state_;
std::unique_ptr<EclipseWriter> eclipse_writer_;
std::unique_ptr<OutputWriter> output_writer_;
boost::any parallel_information_;
std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
std::unique_ptr<Simulator> simulator_;
std::string logFile_;
};
} // namespace Opm

View File

@ -73,18 +73,42 @@ namespace Opm
class GridInit<Dune::CpGrid>
{
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<double>& 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

View File

@ -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<EclipseState> eclipse_state,
const EclipseState& eclState,
BlackoilOutputWriterEbos& output_writer,
const std::vector<double>& threshold_pressures_by_face)
const std::unordered_set<std::string>& 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<RateConverterType> rateConverter_;
// Threshold pressures.
std::vector<double> 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<std::string> defunct_well_names_;
// Whether this a parallel simulation or not
bool is_parallel_run_;

View File

@ -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<double> residual() {
if( ! wellsActive() )
{
return std::vector<double>();
}
const int np = numPhases();
const int nw = wells().number_of_wells;