Merge branch 'master' into frankenstein

* master: (42 commits)
  Let only one rank write to step_timing.txt
  Do not refer users to issue tracker if multiple procs log.
  Remove unused variable.
  Use vector instead of VLA, also add missing includes.
  changed: bundle eigen3 in the original tarball for debian
  update redhat6 packaging
  Bugfix parallel computation of weighted pressure etc.
  Fixed uninitialized bug, and added logging/comment
  Removed superfluous std::move
  Refactoring
  Initial version of summary data
  Do not store collective communication in the wells object.
  Make sure that updateWellControls is called on each process.
  Make WellSwitchingLogger work with DUNE 2.3
  Schedule::getGroup returns reference, not pointer
  Removed warning in WellSwitchLogger::calculateMessageSize
  Correctly initialize MPI for multisegment wells test
  Changed some names in WellSwitchingLogger
  Use speaking name for bool in getCellData
  Whitespace and other formatting changes
  ...
This commit is contained in:
Andreas Lauser 2016-10-14 19:31:56 +02:00
commit 44d3d5b536
26 changed files with 1151 additions and 184 deletions

View File

@ -67,6 +67,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/polymer/fullyimplicit/PolymerPropsAd.cpp
opm/simulators/SimulatorCompressibleTwophase.cpp
opm/simulators/SimulatorIncompTwophase.cpp
opm/simulators/WellSwitchingLogger.cpp
)
@ -88,6 +89,7 @@ list (APPEND TEST_SOURCE_FILES
tests/test_solventprops_ad.cpp
tests/test_multisegmentwells.cpp
# tests/test_thresholdpressure.cpp
tests/test_wellswitchlogger.cpp
)
list (APPEND TEST_DATA_FILES
@ -241,7 +243,9 @@ list (APPEND PUBLIC_HEADER_FILES
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp
opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp
opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp
opm/simulators/ParallelFileMerger.hpp
opm/simulators/SimulatorCompressibleTwophase.hpp
opm/simulators/SimulatorIncompTwophase.hpp
opm/simulators/WellSwitchingLogger.hpp
)

2
debian/control vendored
View File

@ -7,7 +7,7 @@ Build-Depends: build-essential, debhelper (>= 9), libboost-filesystem-dev,
libdune-common-dev, libdune-istl-dev, cmake, libtinyxml-dev, bc,
libert.ecl-dev, git, zlib1g-dev, libtool, doxygen,
texlive-latex-extra, texlive-latex-recommended, ghostscript,
libopm-core-dev, libeigen3-dev (>= 3.2.0), libopm-material-dev,
libopm-core-dev, libopm-material-dev,
libopm-parser-dev, libboost-iostreams-dev, libopm-common-dev,
libopm-grid-dev, libdune-grid-dev, libug-dev, libopm-output-dev,
libtrilinos-zoltan-dev, libopenmpi-dev, mpi-default-bin

17
debian/mk_orig_source vendored Executable file
View File

@ -0,0 +1,17 @@
#!/bin/bash
ORIGDIR=`pwd`
TOPDIR=`dirname $0`
tmp=`mktemp -d`
pushd .
cd $TOPDIR/..
git archive --prefix opm-simulators-$1/ -o $tmp/opm-simulators.tar.gz $2
cd $tmp
wget https://github.com/OPM/eigen3/archive/master.tar.gz -O eigen3-master.tar.gz
tar zxvf opm-simulators.tar.gz
cd opm-simulators-$1
tar zxvf ../eigen3-master.tar.gz
cd ..
tar zcvf $ORIGDIR/opm-simulators-$1.tar.gz opm-simulators-$1/
popd
rm $tmp -Rf

5
debian/rules vendored
View File

@ -21,7 +21,10 @@ override_dh_auto_build:
# consider using -DUSE_VERSIONED_DIR=ON if backporting
override_dh_auto_configure:
dh_auto_configure --buildsystem=cmake -- -DCMAKE_BUILD_TYPE=Release -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_INSTALL_DOCDIR=share/doc/libopm-simulators1 -DWHOLE_PROG_OPTIM=OFF -DUSE_RUNPATH=OFF
mkdir build-eigen
cd build-eigen && cmake $(CURDIR)/eigen3-master -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=/tmp/install-eigen3
cd build-eigen && make install
dh_auto_configure --buildsystem=cmake -- -DCMAKE_BUILD_TYPE=Release -DUSE_MPI=1 -DBUILD_SHARED_LIBS=1 -DCMAKE_INSTALL_DOCDIR=share/doc/libopm-simulators1 -DWHOLE_PROG_OPTIM=OFF -DUSE_RUNPATH=OFF -DCMAKE_PREFIX_PATH=/tmp/install-eigen3
override_dh_auto_install:
dh_auto_install -- install-html

View File

@ -92,9 +92,20 @@ namespace Opm {
struct SimulatorData {
SimulatorData(int num_phases);
enum FipId {
FIP_AQUA = Opm::Phases::Water,
FIP_LIQUID = Opm::Phases::Oil,
FIP_VAPOUR = Opm::Phases::Gas,
FIP_DISSOLVED_GAS = 3,
FIP_VAPORIZED_OIL = 4,
FIP_PV = 5, //< Pore volume
FIP_WEIGHTED_PRESSURE = 6
};
std::vector<ReservoirResidualQuant> rq;
ADB rsSat;
ADB rvSat;
ADB rsSat; // Saturated gas-oil ratio
ADB rvSat; // Saturated oil-gas ratio
std::array<V, 7> fip;
};
typedef typename ModelTraits<Implementation>::ReservoirState ReservoirState;

View File

@ -424,6 +424,7 @@ typedef Eigen::Array<double,
: rq(num_phases)
, rsSat(ADB::null())
, rvSat(ADB::null())
, fip()
{
}
@ -1011,12 +1012,19 @@ typedef Eigen::Array<double,
const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
assert(dx.size() == total_residual_v.size());
asImpl().wellModel().updateWellState(dx.array(), dpMaxRel(), well_state);
asImpl().wellModel().updateWellControls(well_state);
}
// We have to update the well controls regardless whether there are local
// wells active or not as parallel logging will take place that needs to
// communicate with all processes.
asImpl().wellModel().updateWellControls(well_state);
} while (it < 15);
if (converged) {
OpmLog::note("well converged iter: " + std::to_string(it));
if (terminalOutputEnabled())
{
OpmLog::note("well converged iter: " + std::to_string(it));
}
const int nw = wells().number_of_wells;
{
// We will set the bhp primary variable to the new ones,
@ -2164,21 +2172,18 @@ typedef Eigen::Array<double,
const ADB pv_mult = poroMult(pressure);
const V& pv = geo_.poreVolume();
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
std::vector<V> fip(5, V::Zero(nc));
for (int phase = 0; phase < maxnp; ++phase) {
if (active_[ phase ]) {
const int pos = pu.phase_pos[ phase ];
const auto& b = asImpl().fluidReciprocFVF(phase, canonical_phase_pressures[phase], temperature, rs, rv, cond);
fip[phase] = ((pv_mult * b * saturation[pos] * pv).value());
sd_.fip[phase] = ((pv_mult * b * saturation[pos] * pv).value());
}
}
if (active_[ Oil ] && active_[ Gas ]) {
// Account for gas dissolved in oil and vaporized oil
const int po = pu.phase_pos[Oil];
const int pg = pu.phase_pos[Gas];
fip[3] = rs.value() * fip[po];
fip[4] = rv.value() * fip[pg];
sd_.fip[SimulatorData::FIP_DISSOLVED_GAS] = rs.value() * sd_.fip[SimulatorData::FIP_LIQUID];
sd_.fip[SimulatorData::FIP_VAPORIZED_OIL] = rv.value() * sd_.fip[SimulatorData::FIP_VAPOUR];
}
// For a parallel run this is just a local maximum and needs to be updated later
@ -2191,23 +2196,57 @@ typedef Eigen::Array<double,
if ( !isParallel() )
{
for (int i = 0; i < 5; ++i) {
//Accumulate phases for each region
for (int phase = 0; phase < maxnp; ++phase) {
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) {
values[fipnum[c]-1][i] += fip[i][c];
const int region = fipnum[c] - 1;
if (region != -1) {
values[region][phase] += sd_.fip[phase][c];
}
}
}
//Accumulate RS and RV-volumes for each region
if (active_[ Oil ] && active_[ Gas ]) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
}
}
}
hcpv = V::Zero(dims);
pres = V::Zero(dims);
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0) {
hcpv[fipnum[c]-1] += pv[c] * hydrocarbon[c];
pres[fipnum[c]-1] += pv[c] * pressure.value()[c];
values[fipnum[c]-1][5] += pv[c];
values[fipnum[c]-1][6] += pv[c] * pressure.value()[c] * hydrocarbon[c];
const int region = fipnum[c] - 1;
if (region != -1) {
hcpv[region] += pv[c] * hydrocarbon[c];
pres[region] += pv[c] * pressure.value()[c];
}
}
sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
//Compute hydrocarbon pore volume weighted average pressure.
//If we have no hydrocarbon in region, use pore volume weighted average pressure instead
if (hcpv[region] != 0) {
sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
} else {
sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
}
values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
}
}
}
@ -2223,10 +2262,23 @@ typedef Eigen::Array<double,
dims = comm.max(dims);
values.resize(dims, V::Zero(7));
for (int i = 0; i < 5; ++i) {
//Accumulate phases for each region
for (int phase = 0; phase < maxnp; ++phase) {
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0 && mask[c]) {
values[fipnum[c]-1][i] += fip[i][c];
const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) {
values[region][phase] += sd_.fip[phase][c];
}
}
}
//Accumulate RS and RV-volumes for each region
if (active_[ Oil ] && active_[ Gas ]) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) {
values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
}
}
}
@ -2235,11 +2287,32 @@ typedef Eigen::Array<double,
pres = V::Zero(dims);
for (int c = 0; c < nc; ++c) {
if (fipnum[c] != 0 && mask[c]) {
hcpv[fipnum[c]-1] += pv[c] * hydrocarbon[c];
pres[fipnum[c]-1] += pv[c] * pressure.value()[c];
values[fipnum[c]-1][5] += pv[c];
values[fipnum[c]-1][6] += pv[c] * pressure.value()[c] * hydrocarbon[c];
const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) {
hcpv[region] += pv[c] * hydrocarbon[c];
pres[region] += pv[c] * pressure.value()[c];
}
}
comm.sum(hcpv.data(), hcpv.size());
comm.sum(pres.data(), pres.size());
sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1 && mask[c]) {
sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
if (hcpv[region] != 0) {
sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
} else {
sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
}
values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
}
}
@ -2250,24 +2323,12 @@ typedef Eigen::Array<double,
{
comm.sum(values[reg].data(), values[reg].size());
}
comm.sum(hcpv.data(), hcpv.size());
comm.sum(pres.data(), pres.size());
#else
// This should never happen!
OPM_THROW(std::logic_error, "HAVE_MPI should be defined if we are running in parallel");
#endif
}
// compute PAV and PORV for every regions.
for (int reg = 0; reg < dims; ++reg) {
if (hcpv[reg] != 0) {
values[reg][6] /= hcpv[reg];
} else {
values[reg][6] = pres[reg] / values[reg][5];
}
}
return values;
}

View File

@ -253,7 +253,6 @@ namespace Opm {
/// Compute fluid in place.
/// \param[in] ReservoirState
/// \param[in] WellState
/// \param[in] FIPNUM for active cells not global cells.
/// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas.
std::vector<V>

View File

@ -625,6 +625,21 @@ namespace Opm {
wellModel().updateWellState(dwells, dpMaxRel(), well_state);
for( auto w = 0; w < wells().number_of_wells; ++w ) {
if (wells().type[w] == INJECTOR) {
continue;
}
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf ) {
int wc = wells().well_cells[perf];
if ( (ss[wc] + sg[wc]) > 0) {
well_state.solventFraction()[perf] = ss[wc] / (ss[wc] + sg[wc]);
}
}
}
// Update phase conditions used for property calculations.
updatePhaseCondFromPrimalVariable(reservoir_state);
}

View File

@ -40,6 +40,7 @@
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
#include <opm/simulators/ParallelFileMerger.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
@ -85,6 +86,7 @@
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#ifdef _OPENMP
#include <omp.h>
#endif
@ -152,10 +154,21 @@ namespace Opm
asImpl().createSimulator();
// Run.
return asImpl().runSimulator();
auto ret = asImpl().runSimulator();
asImpl().mergeParallelLogFiles();
return ret;
}
catch (const std::exception &e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
std::ostringstream message;
message << "Program threw an exception: " << e.what();
if( output_cout_ )
{
OpmLog::error(message.str());
}
return EXIT_FAILURE;
}
@ -183,6 +196,7 @@ namespace Opm
// members first occur.
// setupParallelism()
int mpi_rank_ = 0;
bool output_cout_ = false;
bool must_distribute_ = false;
// setupParameters()
@ -232,9 +246,9 @@ namespace Opm
// 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();
mpi_rank_ = mpi_helper.rank();
const int mpi_size = mpi_helper.size();
output_cout_ = ( mpi_rank == 0 );
output_cout_ = ( mpi_rank_ == 0 );
must_distribute_ = ( mpi_size > 1 );
#ifdef _OPENMP
@ -253,7 +267,7 @@ namespace Opm
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;
std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank_ << "." << std::endl;
}
}
#endif
@ -340,10 +354,14 @@ namespace Opm
{
// 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.
output_dir_ =
param_.getDefault("output_dir", std::string("."));
boost::filesystem::path fpath(output_dir_);
if (!is_directory(fpath)) {
try {
@ -370,36 +388,80 @@ namespace Opm
using boost::filesystem::path;
path fpath(deck_filename);
std::string baseName;
std::string debugFile;
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")) {
logFile_ = output_dir_ + "/" + baseName + ".PRT";
debugFile = output_dir_ + "/." + baseName + ".DEBUG";
} else {
logFile_ = baseName + ".PRT";
debugFile = "." + baseName + ".DEBUG";
logFileStream << output_dir_ << "/";
debugFileStream << output_dir_ + "/";
}
std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile_ , Log::NoDebugMessageTypes);
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);
std::shared_ptr<StreamLog> debugLog = std::make_shared<EclipsePRTLog>(debugFile, Log::DefaultMessageTypes, false, output_cout_);
OpmLog::addBackend( "DEBUGLOG" , debugLog);
prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(false));
streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10));
streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true));
// Read parameters.
OpmLog::debug("\n--------------- Reading parameters ---------------\n");
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()));
}
// Parser the input and creates the Deck and EclipseState objects.
// Writes to:
@ -418,7 +480,12 @@ namespace Opm
ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }});
deck_ = parser->parseFile(deck_filename, parseContext);
checkDeck(deck_, parser);
MissingFeatures::checkKeywords(*deck_);
if ( output_cout_)
{
MissingFeatures::checkKeywords(*deck_);
}
eclipse_state_.reset(new EclipseState(*deck_, parseContext));
auto ioConfig = eclipse_state_->getIOConfig();
ioConfig->setOutputDir(output_dir_);
@ -619,6 +686,11 @@ namespace Opm
// OpmLog singleton.
void extractMessages()
{
if ( !output_cout_ )
{
return;
}
auto extractMessage = [](const Message& msg) {
auto log_type = detail::convertMessageType(msg.mtype);
const auto& location = msg.location;
@ -649,6 +721,11 @@ namespace Opm
// OpmLog singleton.
void runDiagnostics()
{
if( ! output_cout_ )
{
return;
}
// Run relperm diagnostics
RelpermDiagnostics diagnostic;
diagnostic.diagnosis(eclipse_state_, deck_, grid_init_->grid());

View File

@ -22,6 +22,8 @@
#ifndef OPM_MULTISEGMENTWELLS_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELLS_HEADER_INCLUDED
#include <dune/common/parallel/mpihelper.hh>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
@ -41,6 +43,7 @@
#include <opm/autodiff/WellMultiSegment.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/simulators/WellSwitchingLogger.hpp>
@ -78,6 +81,9 @@ namespace Opm {
Eigen::Dynamic,
Eigen::Dynamic,
Eigen::RowMajor>;
using Communication =
Dune::CollectiveCommunication<typename Dune::MPIHelper
::MPICommunicator>;
// --------- Public methods ---------
// TODO: using a vector of WellMultiSegmentConstPtr for now

View File

@ -824,9 +824,10 @@ namespace Opm
MultisegmentWells::
updateWellControls(WellState& xw) const
{
wellhelpers::WellSwitchingLogger logger;
if( msWells().empty() ) return ;
std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = numPhases();
@ -860,9 +861,9 @@ namespace Opm
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
// Each well is only active on one process. Therefore we always print the sitch info.
std::cout << "Switching control mode for well " << msWells()[w]->name()
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
logger.wellSwitched(msWells()[w]->name(),
well_controls_iget_type(wc, current),
well_controls_iget_type(wc, ctrl_index));
xw.currentControls()[w] = ctrl_index;
current = xw.currentControls()[w];
}

View File

@ -112,7 +112,12 @@ namespace Opm
std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
if( param_.getDefault("timestep.adaptive", true ) )
{
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
if (param_.getDefault("use_TUNING", false)) {
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule->getTuning(), timer.currentStepNum(), param_, terminal_output_ ) );
} else {
adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
}
}
std::string restorefilename = param_.getDefault("restorefile", std::string("") );
@ -276,11 +281,14 @@ namespace Opm
FIPUnitConvert(eclipse_state_->getUnits(), COIP);
V OOIP_totals = FIPTotals(OOIP, state);
V COIP_totals = FIPTotals(COIP, state);
outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
for (size_t reg = 0; reg < OOIP.size(); ++reg) {
outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
}
if ( terminal_output_ )
{
outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
for (size_t reg = 0; reg < OOIP.size(); ++reg) {
outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
}
}
// accumulate total time
stime += st;
@ -299,7 +307,11 @@ namespace Opm
step_report.total_newton_iterations = solver->nonlinearIterations();
step_report.total_linear_iterations = solver->linearIterations();
step_report.total_linearizations = solver->linearizations();
step_report.reportParam(tstep_os);
if ( output_writer_.isIORank() )
{
step_report.reportParam(tstep_os);
}
}
// Increment timer, remember well state.

View File

@ -273,6 +273,12 @@ namespace Opm
/** \brief return true if output is enabled */
bool output () const { return output_; }
/** \brief Whether this process does write to disk */
bool isIORank () const
{
parallelOutput_->isIORank();
}
void restore(SimulatorTimerInterface& timer,
BlackoilState& state,
WellStateFullyImplicitBlackoil& wellState,
@ -423,33 +429,44 @@ namespace Opm
namespace detail {
/**
* Converts an ADB into a standard vector by copy
* Converts an ADB::V into a standard vector by copy
*/
inline std::vector<double> adbToDoubleVector(const Opm::AutoDiffBlock<double>& adb) {
const auto& adb_v = adb.value();
inline std::vector<double> adbVToDoubleVector(const Opm::AutoDiffBlock<double>::V& adb_v) {
std::vector<double> vec(adb_v.data(), adb_v.data() + adb_v.size());
return vec;
}
/**
* Converts an ADB into a standard vector by copy
*/
inline std::vector<double> adbToDoubleVector(const Opm::AutoDiffBlock<double>& adb) {
return adbVToDoubleVector(adb.value());
}
/**
* Returns the data requested in the restartConfig
*/
template<class Model>
std::vector<data::CellData> getCellData(
void getRestartData(
std::vector<data::CellData>& output,
const Opm::PhaseUsage& phaseUsage,
const Model& model,
const Model& physicalModel,
const RestartConfig& restartConfig,
const int reportStepNum) {
const int reportStepNum,
const bool log) {
typedef Opm::AutoDiffBlock<double> ADB;
std::vector<data::CellData> simProps;
const typename Model::SimulatorData& sd = physicalModel.getSimulatorData();
//Get the value of each of the keys
std::map<std::string, int> outKeywords = restartConfig.getRestartKeywords(reportStepNum);
for (auto& keyValue : outKeywords) {
//Get the value of each of the keys for the restart keywords
std::map<std::string, int> rstKeywords = restartConfig.getRestartKeywords(reportStepNum);
for (auto& keyValue : rstKeywords) {
keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
}
const typename Model::SimulatorData& sd = model.getSimulatorData();
//Get shorthands for water, oil, gas
const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
@ -463,161 +480,305 @@ namespace Opm
/**
* Formation volume factors for water, oil, gas
*/
if (aqua_active && outKeywords["BW"] > 0) {
outKeywords["BW"] = 0;
simProps.emplace_back(data::CellData{
if (aqua_active && rstKeywords["BW"] > 0) {
rstKeywords["BW"] = 0;
output.emplace_back(data::CellData{
"1OVERBW",
Opm::UnitSystem::measure::water_inverse_formation_volume_factor,
std::move(adbToDoubleVector(sd.rq[aqua_idx].b))});
adbToDoubleVector(sd.rq[aqua_idx].b),
true});
}
if (liquid_active && outKeywords["BO"] > 0) {
outKeywords["BO"] = 0;
simProps.emplace_back(data::CellData{
if (liquid_active && rstKeywords["BO"] > 0) {
rstKeywords["BO"] = 0;
output.emplace_back(data::CellData{
"1OVERBO",
Opm::UnitSystem::measure::oil_inverse_formation_volume_factor,
std::move(adbToDoubleVector(sd.rq[liquid_idx].b))});
adbToDoubleVector(sd.rq[liquid_idx].b),
true});
}
if (vapour_active && outKeywords["BG"] > 0) {
outKeywords["BG"] = 0;
simProps.emplace_back(data::CellData{
if (vapour_active && rstKeywords["BG"] > 0) {
rstKeywords["BG"] = 0;
output.emplace_back(data::CellData{
"1OVERBG",
Opm::UnitSystem::measure::gas_inverse_formation_volume_factor,
std::move(adbToDoubleVector(sd.rq[vapour_idx].b))});
adbToDoubleVector(sd.rq[vapour_idx].b),
true});
}
/**
* Densities for water, oil gas
*/
if (outKeywords["DEN"] > 0) {
outKeywords["DEN"] = 0;
if (rstKeywords["DEN"] > 0) {
rstKeywords["DEN"] = 0;
if (aqua_active) {
simProps.emplace_back(data::CellData{
output.emplace_back(data::CellData{
"WAT_DEN",
Opm::UnitSystem::measure::density,
std::move(adbToDoubleVector(sd.rq[aqua_idx].rho))});
adbToDoubleVector(sd.rq[aqua_idx].rho),
true});
}
if (liquid_active) {
simProps.emplace_back(data::CellData{
output.emplace_back(data::CellData{
"OIL_DEN",
Opm::UnitSystem::measure::density,
std::move(adbToDoubleVector(sd.rq[liquid_idx].rho))});
adbToDoubleVector(sd.rq[liquid_idx].rho),
true});
}
if (vapour_active) {
simProps.emplace_back(data::CellData{
output.emplace_back(data::CellData{
"GAS_DEN",
Opm::UnitSystem::measure::density,
std::move(adbToDoubleVector(sd.rq[vapour_idx].rho))});
adbToDoubleVector(sd.rq[vapour_idx].rho),
true});
}
}
/**
* Viscosities for water, oil gas
*/
if (outKeywords["VISC"] > 0) {
outKeywords["VISC"] = 0;
if (rstKeywords["VISC"] > 0) {
rstKeywords["VISC"] = 0;
if (aqua_active) {
simProps.emplace_back(data::CellData{
output.emplace_back(data::CellData{
"WAT_VISC",
Opm::UnitSystem::measure::viscosity,
std::move(adbToDoubleVector(sd.rq[aqua_idx].mu))});
adbToDoubleVector(sd.rq[aqua_idx].mu),
true});
}
if (liquid_active) {
simProps.emplace_back(data::CellData{
output.emplace_back(data::CellData{
"OIL_VISC",
Opm::UnitSystem::measure::viscosity,
std::move(adbToDoubleVector(sd.rq[liquid_idx].mu))});
adbToDoubleVector(sd.rq[liquid_idx].mu),
true});
}
if (vapour_active) {
simProps.emplace_back(data::CellData{
output.emplace_back(data::CellData{
"GAS_VISC",
Opm::UnitSystem::measure::viscosity,
std::move(adbToDoubleVector(sd.rq[vapour_idx].mu))});
adbToDoubleVector(sd.rq[vapour_idx].mu),
true});
}
}
/**
* Relative permeabilities for water, oil, gas
*/
if (aqua_active && outKeywords["KRW"] > 0) {
if (aqua_active && rstKeywords["KRW"] > 0) {
if (sd.rq[aqua_idx].kr.size() > 0) {
outKeywords["KRW"] = 0;
simProps.emplace_back(data::CellData{
rstKeywords["KRW"] = 0;
output.emplace_back(data::CellData{
"WATKR",
Opm::UnitSystem::measure::permeability,
std::move(adbToDoubleVector(sd.rq[aqua_idx].kr))});
adbToDoubleVector(sd.rq[aqua_idx].kr),
true});
}
else {
Opm::OpmLog::warning("Empty:WATKR",
"Not emitting empty Water Rel-Perm");
if ( log )
{
Opm::OpmLog::warning("Empty:WATKR",
"Not emitting empty Water Rel-Perm");
}
}
}
if (liquid_active && outKeywords["KRO"] > 0) {
if (liquid_active && rstKeywords["KRO"] > 0) {
if (sd.rq[liquid_idx].kr.size() > 0) {
outKeywords["KRO"] = 0;
simProps.emplace_back(data::CellData{
rstKeywords["KRO"] = 0;
output.emplace_back(data::CellData{
"OILKR",
Opm::UnitSystem::measure::permeability,
std::move(adbToDoubleVector(sd.rq[liquid_idx].kr))});
adbToDoubleVector(sd.rq[liquid_idx].kr),
true});
}
else {
Opm::OpmLog::warning("Empty:OILKR",
"Not emitting empty Oil Rel-Perm");
if ( log )
{
Opm::OpmLog::warning("Empty:OILKR",
"Not emitting empty Oil Rel-Perm");
}
}
}
if (vapour_active && outKeywords["KRG"] > 0) {
if (vapour_active && rstKeywords["KRG"] > 0) {
if (sd.rq[vapour_idx].kr.size() > 0) {
outKeywords["KRG"] = 0;
simProps.emplace_back(data::CellData{
rstKeywords["KRG"] = 0;
output.emplace_back(data::CellData{
"GASKR",
Opm::UnitSystem::measure::permeability,
std::move(adbToDoubleVector(sd.rq[vapour_idx].kr))});
adbToDoubleVector(sd.rq[vapour_idx].kr),
true});
}
else {
Opm::OpmLog::warning("Empty:GASKR",
"Not emitting empty Gas Rel-Perm");
if ( log )
{
Opm::OpmLog::warning("Empty:GASKR",
"Not emitting empty Gas Rel-Perm");
}
}
}
/**
* Vaporized and dissolved gas/oil ratio
*/
if (vapour_active && liquid_active && outKeywords["RSSAT"] > 0) {
outKeywords["RSSAT"] = 0;
simProps.emplace_back(data::CellData{
if (vapour_active && liquid_active && rstKeywords["RSSAT"] > 0) {
rstKeywords["RSSAT"] = 0;
output.emplace_back(data::CellData{
"RSSAT",
Opm::UnitSystem::measure::gas_oil_ratio,
std::move(adbToDoubleVector(sd.rsSat))});
adbToDoubleVector(sd.rsSat),
true});
}
if (vapour_active && liquid_active && outKeywords["RVSAT"] > 0) {
outKeywords["RVSAT"] = 0;
simProps.emplace_back(data::CellData{
if (vapour_active && liquid_active && rstKeywords["RVSAT"] > 0) {
rstKeywords["RVSAT"] = 0;
output.emplace_back(data::CellData{
"RVSAT",
Opm::UnitSystem::measure::oil_gas_ratio,
std::move(adbToDoubleVector(sd.rvSat))});
adbToDoubleVector(sd.rvSat),
true});
}
/**
* Bubble point and dew point pressures
*/
if (vapour_active && liquid_active && outKeywords["PBPD"] > 0) {
outKeywords["PBPD"] = 0;
if (log && vapour_active &&
liquid_active && rstKeywords["PBPD"] > 0) {
rstKeywords["PBPD"] = 0;
Opm::OpmLog::warning("Bubble/dew point pressure output unsupported",
"Writing bubble points and dew points (PBPD) to file is unsupported, "
"as the simulator does not use these internally.");
}
//Warn for any unhandled keyword
for (auto& keyValue : outKeywords) {
if (keyValue.second > 0) {
std::string logstring = "Keyword '";
logstring.append(keyValue.first);
logstring.append("' is unhandled for output to file.");
Opm::OpmLog::warning("Unhandled output keyword", logstring);
if (log) {
for (auto& keyValue : rstKeywords) {
if (keyValue.second > 0) {
std::string logstring = "Keyword '";
logstring.append(keyValue.first);
logstring.append("' is unhandled for output to file.");
Opm::OpmLog::warning("Unhandled output keyword", logstring);
}
}
}
}
return simProps;
/**
* Checks if the summaryConfig has a keyword with the standardized field, region, or block prefixes.
*/
inline bool hasFRBKeyword(const SummaryConfig& summaryConfig, const std::string keyword) {
std::string field_kw = "F" + keyword;
std::string region_kw = "R" + keyword;
std::string block_kw = "B" + keyword;
return summaryConfig.hasKeyword(field_kw)
|| summaryConfig.hasKeyword(region_kw)
|| summaryConfig.hasKeyword(block_kw);
}
/**
* Returns the data as asked for in the summaryConfig
*/
template<class Model>
void getSummaryData(
std::vector<data::CellData>& output,
const Opm::PhaseUsage& phaseUsage,
const Model& physicalModel,
const SummaryConfig& summaryConfig) {
typedef Opm::AutoDiffBlock<double> ADB;
const typename Model::SimulatorData& sd = physicalModel.getSimulatorData();
//Get shorthands for water, oil, gas
const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
/**
* Now process all of the summary config files
*/
// Water in place
if (aqua_active && hasFRBKeyword(summaryConfig, "WIP")) {
output.emplace_back(data::CellData{
"WIP",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_AQUA]),
false});
}
if (liquid_active) {
//Oil in place (liquid phase only)
if (hasFRBKeyword(summaryConfig, "OIPL")) {
output.emplace_back(data::CellData{
"OIPL",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_LIQUID]),
false});
}
//Oil in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "OIPG")) {
output.emplace_back(data::CellData{
"OIPG",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL]),
false});
}
// Oil in place (in liquid and gas phases)
if (hasFRBKeyword(summaryConfig, "OIP")) {
ADB::V oip = sd.fip[Model::SimulatorData::FIP_LIQUID] +
sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL];
output.emplace_back(data::CellData{
"OIP",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(oip),
false});
}
}
if (vapour_active) {
// Gas in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "GIPG")) {
output.emplace_back(data::CellData{
"GIPG",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPOUR]),
false});
}
// Gas in place (liquid phase only)
if (hasFRBKeyword(summaryConfig, "GIPL")) {
output.emplace_back(data::CellData{
"GIPL",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS]),
false});
}
// Gas in place (in both liquid and gas phases)
if (hasFRBKeyword(summaryConfig, "GIP")) {
ADB::V gip = sd.fip[Model::SimulatorData::FIP_VAPOUR] +
sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS];
output.emplace_back(data::CellData{
"GIP",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(gip),
false});
}
}
// Cell pore volume in reservoir conditions
if (hasFRBKeyword(summaryConfig, "RPV")) {
output.emplace_back(data::CellData{
"RPV",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_PV]),
false});
}
// Pressure averaged value (hydrocarbon pore volume weighted)
if (summaryConfig.hasKeyword("FPRH") || summaryConfig.hasKeyword("RPRH")) {
output.emplace_back(data::CellData{
"PRH",
Opm::UnitSystem::measure::pressure,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_WEIGHTED_PRESSURE]),
false});
}
}
}
@ -635,8 +796,15 @@ namespace Opm
bool substep)
{
const RestartConfig& restartConfig = eclipseState_->getRestartConfig();
const SummaryConfig& summaryConfig = eclipseState_->getSummaryConfig();
const int reportStepNum = timer.reportStepNum();
std::vector<data::CellData> cellData = detail::getCellData( phaseUsage_, physicalModel, restartConfig, reportStepNum );
bool logMessages = output_ && parallelOutput_->isIORank();
std::vector<data::CellData> cellData;
detail::getRestartData( cellData, phaseUsage_, physicalModel,
restartConfig, reportStepNum, logMessages );
detail::getSummaryData( cellData, phaseUsage_, physicalModel, summaryConfig );
writeTimeStepWithCellProperties(timer, localState, localWellState, cellData, substep);
}
}

View File

@ -22,6 +22,8 @@
#ifndef OPM_STANDARDWELLS_HEADER_INCLUDED
#define OPM_STANDARDWELLS_HEADER_INCLUDED
#include <dune/common/parallel/mpihelper.hh>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
@ -40,7 +42,7 @@
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/autodiff/VFPProperties.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/simulators/WellSwitchingLogger.hpp>
namespace Opm {
@ -58,7 +60,9 @@ namespace Opm {
// --------- Types ---------
using ADB = AutoDiffBlock<double>;
using Vector = ADB::V;
using V = ADB::V;
using Communication =
Dune::CollectiveCommunication<typename Dune::MPIHelper
::MPICommunicator>;
// copied from BlackoilModelBase
// should put to somewhere better

View File

@ -72,7 +72,7 @@ namespace Opm
StandardWells::StandardWells(const Wells* wells_arg)
StandardWells::StandardWells(const Wells* wells_arg)
: wells_active_(wells_arg!=nullptr)
, wells_(wells_arg)
, wops_(wells_arg)
@ -708,9 +708,10 @@ namespace Opm
StandardWells::
updateWellControls(WellState& xw) const
{
wellhelpers::WellSwitchingLogger logger;
if( !localWellsActive() ) return ;
std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = wells().number_of_phases;
@ -745,11 +746,9 @@ namespace Opm
// Constraint number ctrl_index was broken, switch to it.
// We disregard terminal_ouput here as with it only messages
// for wells on one process will be printed.
std::ostringstream ss;
ss << "Switching control mode for well " << wells().name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
OpmLog::info(ss.str());
logger.wellSwitched(wells().name[w],
well_controls_iget_type(wc, current),
well_controls_iget_type(wc, ctrl_index));
xw.currentControls()[w] = ctrl_index;
current = xw.currentControls()[w];
}

View File

@ -35,6 +35,8 @@ namespace Opm {
namespace wellhelpers
{
inline
double rateToCompare(const std::vector<double>& well_phase_flow_rate,
const int well,

View File

@ -33,6 +33,37 @@ namespace Opm
/// One solvent fraction per well connection
std::vector<double>& solventFraction() { return solvent_fraction_; }
const std::vector<double>& solventFraction() const { return solvent_fraction_; }
data::Wells report(const PhaseUsage &pu) const override {
data::Wells res = WellStateFullyImplicitBlackoil::report(pu);
const int nw = WellState::numWells();
// If there are now wells numPhases throws a floating point
// exception.
if (nw == 0) {
return res;
}
const int np = BaseType::numPhases();
assert( np == 3 ); // the solvent model assumes 3 phases in the base model
// completions aren't supported yet
for( auto w = 0; w < nw; ++w ) {
using rt = data::Rates::opt;
double solvent_well_rate = 0.0;
for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
auto solvent_rate_this = BaseType::perfPhaseRates()[np*perf + pu.phase_pos[BlackoilPhases::Vapour]] * solventFraction()[perf];
solvent_well_rate += solvent_rate_this;
}
res.at( wells_->name[ w ]).rates.set( rt::solvent, solvent_well_rate );
}
return res;
}
private:
std::vector<double> solvent_fraction_;
};

View File

@ -33,6 +33,7 @@
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/polymer/PolymerBlackoilState.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/well_controls.h>
@ -555,7 +556,7 @@ namespace {
const std::vector<ADB>& sat = state.saturation;
const std::vector<PhasePresence> cond = phaseCondition();
std::vector<ADB> pressure = computePressures(state);
std::vector<ADB> pressure = computePressures(state);
const ADB pv_mult = poroMult(press);
const V& pv = geo_.poreVolume();
@ -596,6 +597,18 @@ namespace {
}
}
// Fluid in place is not implemented in this class.
// See BlackoilModelBase::computeFluidInPlace(...) for how it's implemented there
// FIXME: This following code has not been tested.
if (sd_.fip[0].size() == 0) {
OpmLog::warning("NOT_COMPUTING_FIP",
"Computing fluid in place is not implemented for summary files.");
for (int i = 0; i < 7; ++i) {
sd_.fip[i] = V::Zero(nc);
}
}
return values;
}

View File

@ -80,11 +80,23 @@ namespace Opm {
: rq(num_phases)
, rsSat(ADB::null())
, rvSat(ADB::null())
, fip()
{
}
enum FipId {
FIP_AQUA = BlackoilPropsAdInterface::Water,
FIP_LIQUID = BlackoilPropsAdInterface::Oil,
FIP_VAPOUR = BlackoilPropsAdInterface::Gas,
FIP_DISSOLVED_GAS = 3,
FIP_VAPORIZED_OIL = 4,
FIP_PV = 5, //< Pore volume
FIP_WEIGHTED_PRESSURE = 6
};
std::vector<ReservoirResidualQuant> rq;
ADB rsSat;
ADB rvSat;
std::array<V, 7> fip;
};
/// Construct a solver. It will retain references to the

View File

@ -0,0 +1,132 @@
/*
Copyright 2016 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2016 STATOIL AS.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PARALLELFILEMERGER_HEADER_INCLUDED
#define OPM_PARALLELFILEMERGER_HEADER_INCLUDED
#include <memory>
#include <boost/filesystem.hpp>
#include <boost/filesystem/fstream.hpp>
#include <boost/regex.hpp>
namespace Opm
{
namespace detail
{
namespace fs = boost::filesystem;
/// \brief A functor that merges multiple files of a parallel run to one file.
///
/// Without care multiple processes might log messages in a parallel run.
/// Non-root processes will do that to seperate files
/// <basename>.<rank>.<extension. This functor will append those file
/// to usual ones and delete the other files.
class ParallelFileMerger
{
public:
/// \brief Constructor
/// \param output_dir The output directory to use for reading/Writing.
/// \param deckanme The name of the deck.
ParallelFileMerger(const fs::path& output_dir,
const std::string& deckname)
: debugFileRegex_("\\."+deckname+"\\.\\d+\\.DEBUG"),
logFileRegex_(deckname+"\\.\\d+\\.PRT")
{
auto debugPath = output_dir;
debugPath /= (std::string(".") + deckname + ".DEBUG");
debugStream_.reset(new fs::ofstream(debugPath,
std::ofstream::app));
auto logPath = output_dir;
logPath /= ( deckname + ".PRT");
logStream_.reset(new fs::ofstream(logPath,
std::ofstream::app));
}
void operator()(const fs::path& file)
{
const static boost::regex regex(".+\\.(\\d+)\\..+");
boost::smatch matches;
std::string filename = file.filename().native();
if ( boost::regex_match(filename, matches, regex) )
{
std::string rank = boost::regex_replace(filename, regex, "\\1");
if( boost::regex_match(filename, logFileRegex_) )
{
appendFile(*logStream_, file, rank);
}
else
{
if (boost::regex_match(filename, debugFileRegex_) )
{
appendFile(*debugStream_, file, rank);
}
else
{
OPM_THROW(std::runtime_error,
"Unrecognized file with name "
<< filename
<< " from parallel run.");
}
}
}
}
private:
/// \brief Append contents of a file to a stream
/// \brief of The output stream to use.
/// \brief file The file whose content to append.
/// \brief rank The rank that wrote the file.
void appendFile(fs::ofstream& of, const fs::path& file, const std::string& rank)
{
if( fs::file_size(file) )
{
std::cerr << "WARNING: There has been logging to file "
<< file.string() <<" by process "
<< rank << std::endl;
fs::ifstream in(file);
of<<std::endl<< std::endl;
of<<"=======================================================";
of<<std::endl<<std::endl;
of << " Output written by rank " << rank << " to file " << file.string();
of << ":" << std::endl << std::endl;
of << in.rdbuf() << std::endl << std::endl;
of << "======================== end output =====================";
of << std::endl;
in.close();
}
fs::remove(file);
}
/// \brief Regex to capture .*.DEBUG
boost::regex debugFileRegex_;
/// \brief Regex to capture *.PRT
boost::regex logFileRegex_;
/// \brief Stream to *.DEBUG file
std::unique_ptr<fs::ofstream> debugStream_;
/// \brief Stream to *.PRT file
std::unique_ptr<fs::ofstream> logStream_;
};
} // end namespace detail
} // end namespace OPM
#endif // end header guard

View File

@ -0,0 +1,205 @@
/*
Copyright 2016 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2016 Statoil AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/simulators/WellSwitchingLogger.hpp>
#include <numeric>
namespace Opm
{
namespace wellhelpers
{
#if HAVE_MPI
int WellSwitchingLogger::calculateMessageSize(std::vector<int>& well_name_lengths)
{
// Each process will send a message to the root process with
// the following data:
// total number of switches, for each switch the length of the
// well name, for each switch the well name and the two controls.
well_name_lengths.reserve(switchMap_.size());
for(const auto& switchEntry : switchMap_)
{
int length = switchEntry.first.size() +1; //we write an additional \0
well_name_lengths.push_back(length);
}
// compute the message size
int message_size = 0;
int increment = 0;
// number of switches
MPI_Pack_size(1, MPI_INT, MPI_COMM_WORLD, &message_size);
// const char* length include delimiter for each switch
MPI_Pack_size(switchMap_.size(), MPI_INT, MPI_COMM_WORLD, &increment);
message_size += increment;
// for each well the name + two controls in one write
for(const auto& length : well_name_lengths)
{
// well name
MPI_Pack_size(length, MPI_CHAR, MPI_COMM_WORLD, &increment);
message_size += increment;
// controls
MPI_Pack_size(2, MPI_CHAR, MPI_COMM_WORLD, &increment);
message_size += increment;
}
return message_size;
}
void WellSwitchingLogger::packData(std::vector<int>& well_name_lengths,
std::vector<char>& buffer)
{
// Pack the data
// number of switches
int offset = 0;
int no_switches = switchMap_.size();
MPI_Pack(&no_switches, 1, MPI_INT, buffer.data(), buffer.size(),
&offset, MPI_COMM_WORLD);
MPI_Pack(well_name_lengths.data(), well_name_lengths.size(),
MPI_INT, buffer.data(), buffer.size(),
&offset, MPI_COMM_WORLD);
for(const auto& switchEntry : switchMap_)
{
// well name
auto& well_name = switchEntry.first;
MPI_Pack(const_cast<char*>(well_name.c_str()), well_name.size()+1,
MPI_CHAR, buffer.data(), buffer.size(),
&offset, MPI_COMM_WORLD);
// controls
MPI_Pack(const_cast<char*>(switchEntry.second.data()), 2 , MPI_CHAR,
buffer.data(), buffer.size(), &offset, MPI_COMM_WORLD);
}
}
void WellSwitchingLogger::unpackDataAndLog(std::vector<char>& recv_buffer,
const std::vector<int>& displ)
{
for(int p=1; p < cc_.size(); ++p)
{
int offset = displ[p];
int no_switches = 0;
MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset,
&no_switches, 1, MPI_INT, MPI_COMM_WORLD);
if ( no_switches == 0 )
{
continue;
}
std::vector<int> well_name_lengths(no_switches);
MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset,
well_name_lengths.data(), well_name_lengths.size(),
MPI_INT, MPI_COMM_WORLD);
std::vector<char> well_name;
for ( int i = 0; i < no_switches; ++i )
{
well_name.resize(well_name_lengths[i]);
MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset,
well_name.data(), well_name_lengths[i], MPI_CHAR,
MPI_COMM_WORLD);
std::array<char,2> fromto{{}};
MPI_Unpack(recv_buffer.data(), recv_buffer.size(), &offset,
fromto.data(), 2, MPI_CHAR, MPI_COMM_WORLD);
logSwitch(well_name.data(), fromto, p);
}
}
}
void WellSwitchingLogger::logSwitch(const char* name, std::array<char,2> fromto,
int rank)
{
std::ostringstream ss;
ss << "Switching control mode for well " << name
<< " from " << modestring[WellControlType(fromto[0])]
<< " to " << modestring[WellControlType(fromto[1])]
<< " on rank " << rank << std::endl;
OpmLog::info(ss.str());
}
#endif
void WellSwitchingLogger::gatherDataAndLog()
{
#if HAVE_MPI
if(cc_.size() == 1)
{
return;
}
std::vector<int> message_sizes;
std::vector<int> well_name_lengths;
int message_size = calculateMessageSize(well_name_lengths);
if ( cc_.rank() == 0 ){
for(const auto& entry : switchMap_)
{
logSwitch(entry.first.c_str(), entry.second,0);
}
message_sizes.resize(cc_.size());
}
MPI_Gather(&message_size, 1, MPI_INT, message_sizes.data(),
1, MPI_INT, 0, MPI_COMM_WORLD);
std::vector<char> buffer(message_size);
packData(well_name_lengths, buffer);
std::vector<int> displ;
if ( cc_.rank() == 0){
// last entry will be total size of
displ.resize(cc_.size() + 1, 0);
std::partial_sum(message_sizes.begin(), message_sizes.end(),
displ.begin()+1);
}
std::vector<char> recv_buffer;
if ( cc_.rank() == 0 ){
recv_buffer.resize(displ[cc_.size()]);
}
MPI_Gatherv(buffer.data(), buffer.size(), MPI_PACKED,
recv_buffer.data(), message_sizes.data(),
displ.data(), MPI_PACKED, 0, MPI_COMM_WORLD);
if ( cc_.rank() == 0 )
{
unpackDataAndLog(recv_buffer, displ);
}
#endif
}
WellSwitchingLogger::~WellSwitchingLogger()
{
gatherDataAndLog();
}
} // end namespace wellhelpers
} // end namespace Opm

View File

@ -0,0 +1,114 @@
/*
Copyright 2016 Dr. Blatt - HPC-Simulation-Software & Services
Copyright 2016 Statoil AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLSWITCHINGLOGGER_HEADER_INCLUDED
#define OPM_WELLSWITCHINGLOGGER_HEADER_INCLUDED
#include <array>
#include <map>
#include <string>
#include <vector>
#include <dune/common/parallel/mpihelper.hh>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/core/well_controls.h>
namespace Opm
{
namespace wellhelpers
{
/// \brief Utility class to handle the log messages about well switching.
///
/// In parallel all the messages will be send to a root processor
/// and logged there.
class WellSwitchingLogger
{
typedef std::map<std::string, std::array<char,2> > SwitchMap;
public:
/// \brief The type of the collective communication used.
typedef Dune::CollectiveCommunication<typename Dune::MPIHelper::MPICommunicator>
Communication;
/// \brief Constructor.
///
/// \param cc The collective communication to use.
explicit WellSwitchingLogger(const Communication& cc =
Dune::MPIHelper::getCollectiveCommunication())
: cc_(cc)
{}
/// \brief Log that a well switched.
/// \param name The name of the well.
/// \param from The control of the well before the switch.
/// \param to The control of the well after the switch.
void wellSwitched(std::string name,
WellControlType from,
WellControlType to)
{
if( cc_.size() > 1 )
{
using Pair = typename SwitchMap::value_type;
switchMap_.insert(Pair(name, {char(from), char(to)}));
}
else
{
std::ostringstream ss;
ss << "Switching control mode for well " << name
<< " from " << modestring[from]
<< " to " << modestring[to] << std::endl;
OpmLog::info(ss.str());
}
}
/// \brief Destructor send does the actual logging.
~WellSwitchingLogger();
private:
#if HAVE_MPI
void unpackDataAndLog(std::vector<char>& recv_buffer,
const std::vector<int>& displ);
void packData(std::vector<int>& well_name_length,
std::vector<char>& buffer);
int calculateMessageSize(std::vector<int>& well_name_length);
void logSwitch(const char* name, std::array<char,2> fromto,
int rank);
#endif // HAVE_MPI
void gatherDataAndLog();
/// \brief A map containing the local switches
SwitchMap switchMap_;
/// \brief Collective communication object.
Communication cc_;
/// \brief The strings for printing.
const std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
};
} // end namespace wellhelpers
} // end namespace Opm
#endif

View File

@ -187,7 +187,7 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
int pvtRegionIdx = pvtRegion[cellIdx];
double T = initialState.temperature()[cellIdx];
double p = initialState.pressure()[cellIdx];
double p = phasePressure[wpos][cellIdx];
double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b;
@ -201,21 +201,22 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
int pvtRegionIdx = pvtRegion[cellIdx];
double T = initialState.temperature()[cellIdx];
double p = initialState.pressure()[cellIdx];
double p = phasePressure[opos][cellIdx];
double Rs = initialState.gasoilratio()[cellIdx];
double RsSat = pvto.saturatedGasDissolutionFactor(pvtRegionIdx, T, p);
double b;
if (Rs >= RsSat) {
double b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
}
else {
double b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
if (pu.phase_used[BlackoilPhases::Vapour]) {
int gpos = pu.phase_pos[BlackoilPhases::Vapour];
rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b;
}
b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
}
rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
if (pu.phase_used[BlackoilPhases::Vapour]) {
int gpos = pu.phase_pos[BlackoilPhases::Vapour];
rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b;
}
}
}
@ -227,22 +228,21 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
int pvtRegionIdx = pvtRegion[cellIdx];
double T = initialState.temperature()[cellIdx];
double p = initialState.pressure()[cellIdx];
double p = phasePressure[gpos][cellIdx];
double Rv = initialState.rv()[cellIdx];
double RvSat = pvtg.saturatedOilVaporizationFactor(pvtRegionIdx, T, p);
double b;
if (Rv >= RvSat) {
double b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
}
else {
double b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
if (pu.phase_used[BlackoilPhases::Liquid]) {
int opos = pu.phase_pos[BlackoilPhases::Liquid];
rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b;
}
b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
}
rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
if (pu.phase_used[BlackoilPhases::Liquid]) {
int opos = pu.phase_pos[BlackoilPhases::Liquid];
rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b;
}
}
}
@ -270,7 +270,7 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
// update the maximum pressure potential difference between the two
// regions
const auto barrierId = std::make_pair(eq1, eq2);
const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2));
if (maxDp.count(barrierId) == 0) {
maxDp[barrierId] = 0.0;
}
@ -278,7 +278,6 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const double z1 = UgGridHelpers::cellCenterDepth(grid, c1);
const double z2 = UgGridHelpers::cellCenterDepth(grid, c2);
const double zAvg = (z1 + z2)/2; // average depth
const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2;
@ -289,8 +288,8 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
const double sResid2 = minSat[numPhases*c2 + phaseIdx];
// compute gravity corrected pressure potentials at the average depth
const double p1 = phasePressure[phaseIdx][c1] + rhoAvg*gravity*(zAvg - z1);
const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(zAvg - z2);
const double p1 = phasePressure[phaseIdx][c1];
const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(z1 - z2);
if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2))
maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));
@ -356,8 +355,11 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
// set the threshold pressure for faces of PVT regions where the third item
// has been defaulted to the maximum pressure potential difference between
// these regions
const auto barrierId = std::make_pair(eq1, eq2);
thpres_vals[face] = maxDp.at(barrierId);
const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2));
if (maxDp.count(barrierId) > 0)
thpres_vals[face] = maxDp.at(barrierId);
else
thpres_vals[face] = 0.0;
}
}

View File

@ -16,7 +16,7 @@ BuildRequires: blas-devel lapack-devel dune-common-devel opm-output-devel
BuildRequires: git suitesparse-devel doxygen bc
BuildRequires: opm-parser-devel opm-core-devel opm-grid-devel
BuildRequires: tinyxml-devel dune-istl-devel eigen3-devel ert.ecl-devel
%{?el6:BuildRequires: cmake28 devtoolset-2 boost148-devel}
%{?el6:BuildRequires: cmake28 devtoolset-3-toolchain boost148-devel}
%{!?el6:BuildRequires: cmake gcc gcc-c++ boost-devel}
BuildRoot: %{_tmppath}/%{name}-%{version}-build
Requires: libopm-simulators1 = %{version}
@ -59,8 +59,8 @@ This package contains the applications for opm-simulators
%setup -q -n %{name}-release-%{version}-%{tag}
%build
%{?el6:scl enable devtoolset-2 bash}
%{?el6:cmake28} %{?!el6:cmake} -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix} -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF %{?el6:-DCMAKE_CXX_COMPILER=/opt/rh/devtoolset-2/root/usr/bin/g++ -DCMAKE_C_COMPILER=/opt/rh/devtoolset-2/root/usr/bin/gcc -DBOOST_LIBRARYDIR=%{_libdir}/boost148 -DBOOST_INCLUDEDIR=%{_includedir}/boost148}
%{?el6:scl enable devtoolset-3 bash}
%{?el6:cmake28} %{?!el6:cmake} -DBUILD_SHARED_LIBS=1 -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=%{_prefix} -DCMAKE_INSTALL_DOCDIR=share/doc/%{name}-%{version} -DUSE_RUNPATH=OFF %{?el6:-DCMAKE_CXX_COMPILER=/opt/rh/devtoolset-3/root/usr/bin/g++ -DCMAKE_C_COMPILER=/opt/rh/devtoolset-3/root/usr/bin/gcc -DBOOST_LIBRARYDIR=%{_libdir}/boost148 -DBOOST_INCLUDEDIR=%{_includedir}/boost148}
make
%install

View File

@ -26,6 +26,7 @@
#endif
#define BOOST_TEST_MODULE MultisegmentWellsTest
#define BOOST_TEST_NO_MAIN
#include <vector>
#include <unordered_set>
@ -167,3 +168,16 @@ BOOST_AUTO_TEST_CASE(testStructure)
BOOST_CHECK_EQUAL(0, ms_wells->topWellSegments()[0]);
BOOST_CHECK_EQUAL(1, ms_wells->topWellSegments()[1]);
}
bool
init_unit_test_func()
{
return true;
}
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
boost::unit_test::unit_test_main(&init_unit_test_func,
argc, argv);
}

View File

@ -0,0 +1,65 @@
#include <config.h>
#include <dune/common/version.hh>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define BOOST_TEST_MODULE DistributedCpGridTests
#define BOOST_TEST_NO_MAIN
#include <boost/test/unit_test.hpp>
#include <opm/simulators/WellSwitchingLogger.hpp>
#if HAVE_MPI
class MPIError {
public:
/** @brief Constructor. */
MPIError(std::string s, int e) : errorstring(s), errorcode(e){}
/** @brief The error string. */
std::string errorstring;
/** @brief The mpi error code. */
int errorcode;
};
void MPI_err_handler(MPI_Comm *, int *err_code, ...){
char *err_string=new char[MPI_MAX_ERROR_STRING];
int err_length;
MPI_Error_string(*err_code, err_string, &err_length);
std::string s(err_string, err_length);
std::cerr << "An MPI Error ocurred:"<<std::endl<<s<<std::endl;
delete[] err_string;
throw MPIError(s, *err_code);
}
#endif
bool
init_unit_test_func()
{
return true;
}
BOOST_AUTO_TEST_CASE(wellswitchlog)
{
auto cc = Dune::MPIHelper::getCollectiveCommunication();
Opm::wellhelpers::WellSwitchingLogger logger(cc);
std::ostringstream name;
name <<"Well on rank "<<cc.rank()<<std::flush;
logger.wellSwitched(name.str(), BHP, THP);
}
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
#if HAVE_MPI
// register a throwing error handler to allow for
// debugging with "catch throw" in gdb
MPI_Errhandler handler;
MPI_Comm_create_errhandler(MPI_err_handler, &handler);
MPI_Comm_set_errhandler(MPI_COMM_WORLD, handler);
#endif
boost::unit_test::unit_test_main(&init_unit_test_func,
argc, argv);
}