OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-24 10:10:18 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #881 from andlaus/frankenstein_merge_master_v7

Browse Source 
Frankenstein merge master v7
This commit is contained in:
Andreas Lauser 2016-10-28 09:32:06 +02:00 committed by GitHub
commit 54bb57cce0
11 changed files with 191 additions and 72 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -1744,16 +1744,28 @@ typedef Eigen::Array<double,
if (std::isnan(mass_balance_residual[idx])
|| std::isnan(CNV[idx])
|| (idx < np && std::isnan(well_flux_residual[idx]))) {
OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx));
const auto msg = std::string("NaN residual for phase ") + materialName(idx);
if (terminal_output_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
}
if (mass_balance_residual[idx] > maxResidualAllowed()
|| CNV[idx] > maxResidualAllowed()
|| (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx));
const auto msg = std::string("Too large residual for phase ") + materialName(idx);
if (terminal_output_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
}
}
if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) {
OPM_THROW(Opm::NumericalProblem, "NaN or too large residual for well control equation");
const auto msg = std::string("NaN or too large residual for well control equation");
if (terminal_output_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
}
return converged;
@ -1809,10 +1821,18 @@ typedef Eigen::Array<double,
// if one of the residuals is NaN, throw exception, so that the solver can be restarted
for (int idx = 0; idx < np; ++idx) {
if (std::isnan(well_flux_residual[idx])) {
OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx));
const auto msg = std::string("NaN residual for phase ") + materialName(idx);
if (terminal_output_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
}
if (well_flux_residual[idx] > maxResidualAllowed()) {
OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx));
const auto msg = std::string("Too large residual for phase ") + materialName(idx);
if (terminal_output_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
}
}
@ -2198,10 +2218,12 @@ typedef Eigen::Array<double,
{
//Accumulate phases for each region
for (int phase = 0; phase < maxnp; ++phase) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
values[region][phase] += sd_.fip[phase][c];
if (active_[ phase ]) {
for (int c = 0; c < nc; ++c) {
const int region = fipnum[c] - 1;
if (region != -1) {
values[region][phase] += sd_.fip[phase][c];
}
}
}
}

4
opm/autodiff/BlackoilSequentialModel.hpp
View File

@ -247,7 +247,7 @@ namespace Opm {
/// Return the well model
const WellModel& wellModel() const
{
return pressure_model_->wellModel();
return pressure_solver_.model().wellModel();
}
@ -265,7 +265,7 @@ namespace Opm {
/// Return reservoir simulation data (for output functionality)
const SimulatorData& getSimulatorData() const {
return transport_model_->getSimulatorData();
return transport_solver_.model().getSimulatorData();
}

71
opm/autodiff/Compat.hpp
View File

@ -21,10 +21,13 @@
#ifndef OPM_SIMULATORS_COMPAT_HPP
#define OPM_SIMULATORS_COMPAT_HPP
#include <algorithm>
#include <cassert>
#include <opm/common/data/SimulationDataContainer.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilSolventState.hpp>
#include <opm/output/data/Cells.hpp>
#include <opm/output/data/Solution.hpp>
@ -167,14 +170,68 @@ inline void solutionToSim( const data::Solution& sol,
inline void wellsToState( const data::Wells& wells,
PhaseUsage phases,
WellStateFullyImplicitBlackoil& state ) {
using rt = data::Rates::opt;
inline void wellsToState( const data::Wells& wells, WellState& state ) {
state.bhp() = wells.bhp;
state.temperature() = wells.temperature;
state.wellRates() = wells.well_rate;
state.perfPress() = wells.perf_pressure;
state.perfRates() = wells.perf_rate;
const auto np = phases.num_phases;
std::vector< rt > phs( np );
if( phases.phase_used[BlackoilPhases::Aqua] ) {
phs.at( phases.phase_pos[BlackoilPhases::Aqua] ) = rt::wat;
}
if( phases.phase_used[BlackoilPhases::Liquid] ) {
phs.at( phases.phase_pos[BlackoilPhases::Liquid] ) = rt::oil;
}
if( phases.phase_used[BlackoilPhases::Vapour] ) {
phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas;
}
for( const auto& wm : state.wellMap() ) {
const auto well_index = wm.second[ 0 ];
const auto& well = wells.at( wm.first );
state.bhp()[ well_index ] = well.bhp;
state.temperature()[ well_index ] = well.temperature;
state.currentControls()[ well_index ] = well.control;
const auto wellrate_index = well_index * np;
for( size_t i = 0; i < phs.size(); ++i ) {
assert( well.rates.has( phs[ i ] ) );
state.wellRates()[ wellrate_index + i ] = well.rates.get( phs[ i ] );
}
const auto perforation_pressure = []( const data::Completion& comp ) {
return comp.pressure;
};
const auto perforation_reservoir_rate = []( const data::Completion& comp ) {
return comp.reservoir_rate;
};
std::transform( well.completions.begin(),
well.completions.end(),
state.perfPress().begin() + wm.second[ 1 ],
perforation_pressure );
std::transform( well.completions.begin(),
well.completions.end(),
state.perfRates().begin() + wm.second[ 1 ],
perforation_reservoir_rate );
int local_comp_index = 0;
for (const data::Completion& comp : well.completions) {
const int global_comp_index = wm.second[1] + local_comp_index;
for (int phase_index = 0; phase_index < np; ++phase_index) {
state.perfPhaseRates()[global_comp_index*np + phase_index] = comp.rates.get(phs[phase_index]);
}
++local_comp_index;
}
}
}
}

16
opm/autodiff/FlowMain.hpp
View File

@ -141,10 +141,10 @@ namespace Opm
return EXIT_FAILURE;
}
asImpl().setupOutput();
asImpl().setupLogging();
asImpl().readDeckInput();
asImpl().setupGridAndProps();
asImpl().setupLogging();
asImpl().extractMessages();
asImpl().setupGridAndProps();
asImpl().runDiagnostics();
asImpl().setupState();
asImpl().writeInit();
@ -424,8 +424,16 @@ namespace Opm
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 = eclipse_state_->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));
streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10, limits));
streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true));
// Read parameters.
@ -541,7 +549,7 @@ namespace Opm
fluidprops_.reset(new BlackoilPropsAdFromDeck(*deck_, *eclipse_state_, material_law_manager_, grid));
// Rock compressibility.
rock_comp_.reset(new RockCompressibility(*deck_, *eclipse_state_));
rock_comp_.reset(new RockCompressibility(*deck_, *eclipse_state_, output_cout_));
// Gravity.
assert(UgGridHelpers::dimensions(grid) == 3);

12
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@ -174,6 +174,7 @@ namespace Opm
const boost::any& parallelInformation_arg=boost::any())
: iterations_( 0 ),
parallelInformation_(parallelInformation_arg),
isIORank_(isIORank(parallelInformation_arg)),
parameters_( param )
{
}
@ -487,7 +488,11 @@ namespace Opm
// Check for failure of linear solver.
if (!parameters_.ignoreConvergenceFailure_ && !result.converged) {
OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver.");
const std::string msg("Convergence failure for linear solver.");
if (isIORank_) {
OpmLog::problem(msg);
}
OPM_THROW_NOLOG(LinearSolverProblem, msg);
}
// Copy solver output to dx.
@ -509,6 +514,7 @@ namespace Opm
protected:
mutable int iterations_;
boost::any parallelInformation_;
bool isIORank_;
NewtonIterationBlackoilInterleavedParameters parameters_;
}; // end NewtonIterationBlackoilInterleavedImpl
@ -657,7 +663,9 @@ namespace Opm
// Check for failure of linear solver.
if (!result.converged) {
OPM_THROW(LinearSolverProblem, "Convergence failure for linear solver in computePressureIncrement().");
const std::string msg("Convergence failure for linear solver in computePressureIncrement().");
OpmLog::problem(msg);
OPM_THROW_NOLOG(LinearSolverProblem, msg);
}
// Copy solver output to dx.

20
opm/autodiff/NewtonIterationUtilities.cpp
View File

@ -22,6 +22,7 @@
#include <opm/autodiff/NewtonIterationUtilities.hpp>
#include <opm/autodiff/AutoDiffHelpers.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/utility/platform_dependent/disable_warnings.h>
@ -288,5 +289,24 @@ namespace Opm
/// Return true if this is a serial run, or rank zero on an MPI run.
bool isIORank(const boost::any& parallel_info)
{
#if HAVE_MPI
if (parallel_info.type() == typeid(ParallelISTLInformation)) {
const ParallelISTLInformation& info =
boost::any_cast<const ParallelISTLInformation&>(parallel_info);
return info.communicator().rank() == 0;
} else {
return true;
}
#else
static_cast<void>(parallel_info); // Suppress unused argument warning.
return true;
#endif
}
} // namespace Opm

3
opm/autodiff/NewtonIterationUtilities.hpp
View File

@ -21,6 +21,7 @@
#define OPM_NEWTONITERATIONUTILITIES_HEADER_INCLUDED
#include <opm/autodiff/AutoDiffBlock.hpp>
#include <boost/any.hpp>
#include <vector>
namespace Opm
@ -58,6 +59,8 @@ namespace Opm
Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
AutoDiffBlock<double>::V& b);
/// Return true if this is a serial run, or rank zero on an MPI run.
bool isIORank(const boost::any& parallel_info);
} // namespace Opm

4
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
View File

@ -374,10 +374,12 @@ namespace Opm
if (initConfig.restartRequested() && ((initConfig.getRestartStep()) == (timer.currentStepNum()))) {
std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl;
} else {
data::Solution combined_sol = simToSolution(state, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...)
combined_sol.insert(simProps.begin(), simProps.end()); // ... insert "extra" data (KR, VISC, ...)
eclWriter_->writeTimeStep(timer.reportStepNum(),
substep,
timer.simulationTimeElapsed(),
simToSolution( state, phaseUsage_ ),
combined_sol,
wellState.report(phaseUsage_));
}
}

24
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@ -277,7 +277,7 @@ namespace Opm
/** \brief Whether this process does write to disk */
bool isIORank () const
{
parallelOutput_->isIORank();
return parallelOutput_->isIORank();
}
void restore(SimulatorTimerInterface& timer,
@ -420,7 +420,7 @@ namespace Opm
Opm::UgGridHelpers::numCells(grid) );
solutionToSim( restarted.first, phaseusage, simulatorstate );
wellsToState( restarted.second, wellstate );
wellsToState( restarted.second, phaseusage, wellstate );
}
@ -691,24 +691,26 @@ namespace Opm
data::TargetType::SUMMARY );
}
if (liquid_active) {
const ADB::V& oipl = sd.fip[Model::SimulatorData::FIP_LIQUID];
const ADB::V& oipg = vapour_active ? sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL] : ADB::V();
const ADB::V& oip = vapour_active ? oipl + oipg : oipl;
//Oil in place (liquid phase only)
if (hasFRBKeyword(summaryConfig, "OIPL")) {
output.insert("OIPL",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_LIQUID]),
adbVToDoubleVector(oipl),
data::TargetType::SUMMARY );
}
//Oil in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "OIPG")) {
output.insert("OIPG",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPORIZED_OIL]),
adbVToDoubleVector(oipg),
data::TargetType::SUMMARY );
}
// 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.insert("OIP",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(oip),
@ -716,24 +718,26 @@ namespace Opm
}
}
if (vapour_active) {
const ADB::V& gipg = sd.fip[Model::SimulatorData::FIP_VAPOUR];
const ADB::V& gipl = liquid_active ? sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS] : ADB::V();
const ADB::V& gip = liquid_active ? gipg + gipl : gipg;
// Gas in place (gas phase only)
if (hasFRBKeyword(summaryConfig, "GIPG")) {
output.insert("GIPG",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_VAPOUR]),
adbVToDoubleVector(gipg),
data::TargetType::SUMMARY );
}
// Gas in place (liquid phase only)
if (hasFRBKeyword(summaryConfig, "GIPL")) {
output.insert("GIPL",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(sd.fip[Model::SimulatorData::FIP_DISSOLVED_GAS]),
adbVToDoubleVector(gipl),
data::TargetType::SUMMARY );
}
// 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.insert("GIP",
Opm::UnitSystem::measure::volume,
adbVToDoubleVector(gip),

2
opm/autodiff/SimulatorFullyImplicitBlackoilOutputEbos.hpp
View File

@ -252,7 +252,7 @@ namespace Opm
Opm::UgGridHelpers::numCells(grid) );
solutionToSim( restarted.first, phaseusage, simulatorstate );
wellsToState( restarted.second, wellstate );
wellsToState( restarted.second, phaseusage, wellstate );
}

67
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@ -280,9 +280,23 @@ namespace Opm
data::Wells res = WellState::report(pu);
const int nw = this->numWells();
// If there are now wells numPhases throws a floating point
// exception.
const int np = nw ? this->numPhases() : -1;
if( nw == 0 ) return res;
const int np = pu.num_phases;
using rt = data::Rates::opt;
std::vector< rt > phs( np );
if( pu.phase_used[BlackoilPhases::Aqua] ) {
phs.at( pu.phase_pos[BlackoilPhases::Aqua] ) = rt::wat;
}
if( pu.phase_used[BlackoilPhases::Liquid] ) {
phs.at( pu.phase_pos[BlackoilPhases::Liquid] ) = rt::oil;
}
if( pu.phase_used[BlackoilPhases::Vapour] ) {
phs.at( pu.phase_pos[BlackoilPhases::Vapour] ) = rt::gas;
}
/* this is a reference or example on **how** to convert from
* WellState to something understood by opm-output. it is intended
@ -291,42 +305,23 @@ namespace Opm
* representations.
*/
for( auto w = 0; w < nw; ++w ) {
using rt = data::Rates::opt;
std::map< size_t, data::Completion > completions;
for( const auto& wt : this->wellMap() ) {
const auto w = wt.second[ 0 ];
auto& well = res.at( wt.first );
well.control = this->currentControls()[ w ];
// completions aren't supported yet
//const auto* begin = wells_->well_connpos + w;
//const auto* end = wells_->well_connpos + w + 1;
//for( auto* i = begin; i != end; ++i ) {
// const auto perfrate = this->perfPhaseRates().begin() + *i;
// data::Rates perfrates;
// perfrates.set( rt::wat, *(perfrate + 0) );
// perfrates.set( rt::oil, *(perfrate + 1) );
// perfrates.set( rt::gas, *(perfrate + 2) );
int local_comp_index = 0;
for( auto& comp : well.completions ) {
const auto rates = this->perfPhaseRates().begin()
+ (np * wt.second[ 1 ])
+ (np * local_comp_index);
++local_comp_index;
// const size_t active_index = wells_->well_cells[ *i ];
// completions.emplace( active_index,
// data::Completion{ active_index, perfrates } );
//}
const auto wellrate_index = np * w;
const auto& wv = this->wellRates();
data::Rates wellrates;
if( np == 3 ) {
/* only write if 3-phase solution */
wellrates.set( rt::wat, wv[ wellrate_index + 0 ] );
wellrates.set( rt::oil, wv[ wellrate_index + 1 ] );
wellrates.set( rt::gas, wv[ wellrate_index + 2 ] );
for( int i = 0; i < np; ++i ) {
comp.rates.set( phs[ i ], *(rates + i) );
}
}
const double bhp = this->bhp()[ w ];
const double thp = this->thp()[ w ];
res.emplace( wells_->name[ w ],
data::Well { wellrates, bhp, thp, std::move( completions ) } );
assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]);
}
return res;