Merge branch 'master' into new_well_formulation

Conflicts:
	opm/autodiff/FullyImplicitBlackoilSolver.cpp

To resolve conflicts, WellState was changed to WellStateFullyImplicitBlackoil
in multiple places, and perfRate() changed to perfPhaseRate() in
WellDensitySegmented.
This commit is contained in:
Atgeirr Flø Rasmussen 2014-03-20 10:04:45 +01:00
commit fd4f5f5a26
13 changed files with 206 additions and 84 deletions

View File

@ -106,4 +106,5 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/SimulatorIncompTwophaseAd.hpp opm/autodiff/SimulatorIncompTwophaseAd.hpp
opm/autodiff/TransportSolverTwophaseAd.hpp opm/autodiff/TransportSolverTwophaseAd.hpp
opm/autodiff/WellDensitySegmented.hpp opm/autodiff/WellDensitySegmented.hpp
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
) )

View File

@ -39,7 +39,7 @@
#include <opm/core/linalg/LinearSolverFactory.hpp> #include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp> #include <opm/autodiff/SimulatorFullyImplicitBlackoil.hpp>
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp> #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
@ -200,46 +200,66 @@ try
param.writeParam(output_dir + "/simulation.param"); param.writeParam(output_dir + "/simulation.param");
} }
#warning "TODO: convert the well handling code to the new parser"
#if USE_NEW_PARSER #if USE_NEW_PARSER
std::cout << "\n\n================ Starting main simulation loop ===============\n" std::cout << "\n\n================ Starting main simulation loop ===============\n"
<< std::flush; << std::flush;
WellState well_state; WellStateFullyImplicitBlackoil well_state;
Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck)); Opm::TimeMapPtr timeMap(new Opm::TimeMap(newParserDeck));
SimulatorTimer simtimer; SimulatorTimer simtimer;
std::shared_ptr<EclipseState> eclipseState(new EclipseState(newParserDeck));
// Create new wells, well_state // initialize variables
WellsManager wells(*deck, *grid->c_grid(), props->permeability()); simtimer.init(timeMap, /*beginReportStepIdx=*/0, /*endReportStepIdx=*/0);
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
// since number of wells may change etc.
well_state.init(wells.c_wells(), state);
// Create and run simulator. SimulatorReport fullReport;
SimulatorFullyImplicitBlackoil simulator(param, for (size_t episodeIdx = 0; episodeIdx < timeMap->numTimesteps(); ++episodeIdx) {
*grid->c_grid(), WellsManager wells(eclipseState,
*new_props, episodeIdx,
rock_comp->isActive() ? rock_comp.get() : 0, *grid->c_grid(),
wells, props->permeability());
linsolver,
grav,
outputWriter);
simtimer.init(timeMap);
SimulatorReport report = simulator.run(simtimer, state, well_state);
if (output) { if (episodeIdx == 0) {
report.reportParam(outStream); // @@@ HACK: we should really make a new well state and
warnIfUnusedParams(param); // properly transfer old well state to it every epoch,
// since number of wells may change etc.
well_state.init(wells.c_wells(), state);
}
simtimer.init(timeMap,
/*beginReportStepIdx=*/episodeIdx,
/*endReportStepIdx=*/episodeIdx + 1);
if (episodeIdx == 0)
outputWriter.writeInit(simtimer, state, well_state.basicWellState());
// Create and run simulator.
SimulatorFullyImplicitBlackoil simulator(param,
*grid->c_grid(),
*new_props,
rock_comp->isActive() ? rock_comp.get() : 0,
wells,
linsolver,
grav,
outputWriter);
SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
fullReport += episodeReport;
if (output) {
episodeReport.reportParam(outStream);
}
} }
std::cout << "\n\n================ End of simulation ===============\n\n"; std::cout << "\n\n================ End of simulation ===============\n\n";
report.report(std::cout); fullReport.report(std::cout);
if (output) { if (output) {
std::string filename = output_dir + "/walltime.param"; std::string filename = output_dir + "/walltime.param";
std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
report.reportParam(tot_os); fullReport.reportParam(tot_os);
warnIfUnusedParams(param);
} }
#else #else
std::cout << "\n\n================ Starting main simulation loop ===============\n" std::cout << "\n\n================ Starting main simulation loop ===============\n"
@ -248,12 +268,13 @@ try
SimulatorReport rep; SimulatorReport rep;
// With a deck, we may have more epochs etc. // With a deck, we may have more epochs etc.
WellState well_state; WellStateFullyImplicitBlackoil well_state;
int step = 0; int step = 0;
SimulatorTimer simtimer; SimulatorTimer simtimer;
// Use timer for last epoch to obtain total time. // Use timer for last epoch to obtain total time.
deck->setCurrentEpoch(deck->numberOfEpochs() - 1); deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
simtimer.init(*deck); simtimer.init(*deck);
const double total_time = simtimer.totalTime(); const double total_time = simtimer.totalTime();
for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) { for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
// Set epoch index. // Set epoch index.
@ -285,6 +306,9 @@ try
well_state.init(wells.c_wells(), state); well_state.init(wells.c_wells(), state);
} }
if (epoch == 0)
outputWriter.writeInit(simtimer, state, well_state.basicWellState());
// Create and run simulator. // Create and run simulator.
SimulatorFullyImplicitBlackoil simulator(param, SimulatorFullyImplicitBlackoil simulator(param,
*grid->c_grid(), *grid->c_grid(),
@ -294,6 +318,8 @@ try
linsolver, linsolver,
grav, grav,
outputWriter); outputWriter);
outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
if (epoch == 0) { if (epoch == 0) {
warnIfUnusedParams(param); warnIfUnusedParams(param);
} }

View File

@ -23,6 +23,7 @@
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp> #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/GeoProps.hpp> #include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/BlackoilPropsAd.hpp> #include <opm/autodiff/BlackoilPropsAd.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/wells.h> #include <opm/core/wells.h>
@ -37,7 +38,6 @@
#include <opm/core/utility/Units.hpp> #include <opm/core/utility/Units.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/initState.hpp> #include <opm/core/simulator/initState.hpp>
#include <boost/shared_ptr.hpp> #include <boost/shared_ptr.hpp>
@ -110,7 +110,7 @@ try
initStateBasic(*g, props0, param, 0.0, state); initStateBasic(*g, props0, param, 0.0, state);
initBlackoilSurfvol(*g, props0, state); initBlackoilSurfvol(*g, props0, state);
Opm::WellState well_state; Opm::WellStateFullyImplicitBlackoil well_state;
well_state.init(wells.get(), state); well_state.init(wells.get(), state);
solver.step(1.0, state, well_state); solver.step(1.0, state, well_state);

View File

@ -24,12 +24,12 @@
#include <opm/autodiff/BlackoilPropsAdInterface.hpp> #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/GeoProps.hpp> #include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/WellDensitySegmented.hpp> #include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp> #include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp> #include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/well_controls.h> #include <opm/core/well_controls.h>
@ -223,7 +223,7 @@ namespace {
FullyImplicitBlackoilSolver:: FullyImplicitBlackoilSolver::
step(const double dt, step(const double dt,
BlackoilState& x , BlackoilState& x ,
WellState& xw) WellStateFullyImplicitBlackoil& xw)
{ {
const V pvdt = geo_.poreVolume() / dt; const V pvdt = geo_.poreVolume() / dt;
@ -332,7 +332,7 @@ namespace {
FullyImplicitBlackoilSolver::SolutionState FullyImplicitBlackoilSolver::SolutionState
FullyImplicitBlackoilSolver::constantState(const BlackoilState& x, FullyImplicitBlackoilSolver::constantState(const BlackoilState& x,
const WellState& xw) const WellStateFullyImplicitBlackoil& xw)
{ {
const int nc = grid_.number_of_cells; const int nc = grid_.number_of_cells;
const int np = x.numPhases(); const int np = x.numPhases();
@ -429,7 +429,7 @@ namespace {
FullyImplicitBlackoilSolver::SolutionState FullyImplicitBlackoilSolver::SolutionState
FullyImplicitBlackoilSolver::variableState(const BlackoilState& x, FullyImplicitBlackoilSolver::variableState(const BlackoilState& x,
const WellState& xw) const WellStateFullyImplicitBlackoil& xw)
{ {
const int nc = grid_.number_of_cells; const int nc = grid_.number_of_cells;
const int np = x.numPhases(); const int np = x.numPhases();
@ -615,7 +615,7 @@ namespace {
void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state, void FullyImplicitBlackoilSolver::computeWellConnectionPressures(const SolutionState& state,
const WellState& xw) const WellStateFullyImplicitBlackoil& xw)
{ {
// 1. Compute properties required by computeConnectionPressureDelta(). // 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function // Note that some of the complexity of this part is due to the function
@ -687,7 +687,7 @@ namespace {
FullyImplicitBlackoilSolver:: FullyImplicitBlackoilSolver::
assemble(const V& pvdt, assemble(const V& pvdt,
const BlackoilState& x , const BlackoilState& x ,
const WellState& xw ) const WellStateFullyImplicitBlackoil& xw )
{ {
// Create the primary variables. // Create the primary variables.
const SolutionState state = variableState(x, xw); const SolutionState state = variableState(x, xw);
@ -1148,7 +1148,7 @@ namespace {
void FullyImplicitBlackoilSolver::updateState(const V& dx, void FullyImplicitBlackoilSolver::updateState(const V& dx,
BlackoilState& state, BlackoilState& state,
WellState& well_state) WellStateFullyImplicitBlackoil& well_state)
{ {
const int np = fluid_.numPhases(); const int np = fluid_.numPhases();
const int nc = grid_.number_of_cells; const int nc = grid_.number_of_cells;

View File

@ -33,7 +33,7 @@ namespace Opm {
class RockCompressibility; class RockCompressibility;
class LinearSolverInterface; class LinearSolverInterface;
class BlackoilState; class BlackoilState;
class WellState; class WellStateFullyImplicitBlackoil;
/// A fully implicit solver for the black-oil problem. /// A fully implicit solver for the black-oil problem.
@ -76,7 +76,7 @@ namespace Opm {
void void
step(const double dt , step(const double dt ,
BlackoilState& state , BlackoilState& state ,
WellState& wstate); WellStateFullyImplicitBlackoil& wstate);
private: private:
// Types and enums // Types and enums
@ -149,18 +149,18 @@ namespace Opm {
// Private methods. // Private methods.
SolutionState SolutionState
constantState(const BlackoilState& x, constantState(const BlackoilState& x,
const WellState& xw); const WellStateFullyImplicitBlackoil& xw);
SolutionState SolutionState
variableState(const BlackoilState& x, variableState(const BlackoilState& x,
const WellState& xw); const WellStateFullyImplicitBlackoil& xw);
void void
computeAccum(const SolutionState& state, computeAccum(const SolutionState& state,
const int aix ); const int aix );
void computeWellConnectionPressures(const SolutionState& state, void computeWellConnectionPressures(const SolutionState& state,
const WellState& xw); const WellStateFullyImplicitBlackoil& xw);
void void
addOldWellEq(const SolutionState& state); addOldWellEq(const SolutionState& state);
@ -173,14 +173,14 @@ namespace Opm {
void void
assemble(const V& dtpv, assemble(const V& dtpv,
const BlackoilState& x , const BlackoilState& x,
const WellState& xw ); const WellStateFullyImplicitBlackoil& xw);
V solveJacobianSystem() const; V solveJacobianSystem() const;
void updateState(const V& dx, void updateState(const V& dx,
BlackoilState& state, BlackoilState& state,
WellState& well_state); WellStateFullyImplicitBlackoil& well_state);
std::vector<ADB> std::vector<ADB>
computePressures(const SolutionState& state) const; computePressures(const SolutionState& state) const;

View File

@ -496,13 +496,13 @@ namespace Opm
<< std::endl; << std::endl;
std::cout.precision(8); std::cout.precision(8);
watercut.push(timer.currentTime() + timer.currentStepLength(), watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]), produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init); tot_produced[0]/tot_porevol_init);
if (wells_) { if (wells_) {
wellreport.push(props_, *wells_, wellreport.push(props_, *wells_,
state.pressure(), state.surfacevol(), state.saturation(), state.pressure(), state.surfacevol(), state.saturation(),
timer.currentTime() + timer.currentStepLength(), timer.simulationTimeElapsed() + timer.currentStepLength(),
well_state.bhp(), well_state.perfRates()); well_state.bhp(), well_state.perfRates());
} }
sreport.total_time = step_timer.secsSinceStart(); sreport.total_time = step_timer.secsSinceStart();

View File

@ -29,6 +29,7 @@
#include <opm/autodiff/GeoProps.hpp> #include <opm/autodiff/GeoProps.hpp>
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp> #include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp> #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/wells.h> #include <opm/core/wells.h>
@ -48,7 +49,6 @@
#include <opm/core/grid/ColumnExtract.hpp> #include <opm/core/grid/ColumnExtract.hpp>
#include <opm/core/simulator/BlackoilState.hpp> #include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp> #include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <boost/filesystem.hpp> #include <boost/filesystem.hpp>
@ -77,7 +77,7 @@ namespace Opm
SimulatorReport run(SimulatorTimer& timer, SimulatorReport run(SimulatorTimer& timer,
BlackoilState& state, BlackoilState& state,
WellState& well_state); WellStateFullyImplicitBlackoil& well_state);
private: private:
// Data. // Data.
@ -127,7 +127,7 @@ namespace Opm
SimulatorReport SimulatorFullyImplicitBlackoil::run(SimulatorTimer& timer, SimulatorReport SimulatorFullyImplicitBlackoil::run(SimulatorTimer& timer,
BlackoilState& state, BlackoilState& state,
WellState& well_state) WellStateFullyImplicitBlackoil& well_state)
{ {
return pimpl_->run(timer, state, well_state); return pimpl_->run(timer, state, well_state);
} }
@ -198,7 +198,7 @@ namespace Opm
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n")); std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
} }
} }
static void outputWellStateMatlab(const Opm::WellState& well_state, static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state,
const int step, const int step,
const std::string& output_dir) const std::string& output_dir)
{ {
@ -312,11 +312,8 @@ namespace Opm
SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer, SimulatorReport SimulatorFullyImplicitBlackoil::Impl::run(SimulatorTimer& timer,
BlackoilState& state, BlackoilState& state,
WellState& well_state) WellStateFullyImplicitBlackoil& well_state)
{ {
eclipseWriter_.writeInit(timer, state, well_state);
eclipseWriter_.writeTimeStep(timer, state, well_state);
// Initialisation. // Initialisation.
std::vector<double> porevol; std::vector<double> porevol;
if (rock_comp_props_ && rock_comp_props_->isActive()) { if (rock_comp_props_ && rock_comp_props_->isActive()) {
@ -384,13 +381,9 @@ namespace Opm
// Optionally, check if well controls are satisfied. // Optionally, check if well controls are satisfied.
if (check_well_controls_) { if (check_well_controls_) {
Opm::computePhaseFlowRatesPerWell(*wells_,
well_state.perfRates(),
fractional_flows,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl; std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir // For testing we set surface := reservoir
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase); well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_state.wellRates(), well_state.wellRates());
++well_control_iteration; ++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) { if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries."); OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
@ -424,13 +417,13 @@ namespace Opm
tstep_os.close(); tstep_os.close();
} }
// write an output file for later inspection
if (output_) {
eclipseWriter_.writeTimeStep(timer, state, well_state);
}
// advance to next timestep before reporting at this location // advance to next timestep before reporting at this location
++timer; ++timer;
// write an output file for later inspection
if (output_) {
eclipseWriter_.writeTimeStep(timer, state, well_state.basicWellState());
}
} }
total_timer.stop(); total_timer.stop();

View File

@ -37,7 +37,7 @@ namespace Opm
class LinearSolverInterface; class LinearSolverInterface;
class SimulatorTimer; class SimulatorTimer;
class BlackoilState; class BlackoilState;
class WellState; class WellStateFullyImplicitBlackoil;
struct SimulatorReport; struct SimulatorReport;
/// Class collecting all necessary components for a two-phase simulation. /// Class collecting all necessary components for a two-phase simulation.
@ -84,7 +84,7 @@ namespace Opm
/// \return simulation report, with timing data /// \return simulation report, with timing data
SimulatorReport run(SimulatorTimer& timer, SimulatorReport run(SimulatorTimer& timer,
BlackoilState& state, BlackoilState& state,
WellState& well_state); WellStateFullyImplicitBlackoil& well_state);
private: private:
class Impl; class Impl;

View File

@ -579,12 +579,12 @@ namespace Opm
dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_) dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
.solveGravity(&initial_porevol[0], stepsize, state); .solveGravity(&initial_porevol[0], stepsize, state);
} }
watercut.push(timer.currentTime() + timer.currentStepLength(), watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]), produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init); tot_produced[0]/tot_porevol_init);
if (wells_) { if (wells_) {
wellreport.push(props_, *wells_, state.saturation(), wellreport.push(props_, *wells_, state.saturation(),
timer.currentTime() + timer.currentStepLength(), timer.simulationTimeElapsed() + timer.currentStepLength(),
well_state.bhp(), well_state.perfRates()); well_state.bhp(), well_state.perfRates());
} }
} }

View File

@ -19,22 +19,23 @@
#include <opm/autodiff/WellDensitySegmented.hpp> #include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/core/wells.h> #include <opm/core/wells.h>
#include <opm/core/simulator/WellState.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/props/BlackoilPhases.hpp> #include <opm/core/props/BlackoilPhases.hpp>
#include <numeric> #include <numeric>
#include <cmath> #include <cmath>
std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, std::vector<double>
const WellState& wstate, Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells,
const PhaseUsage& phase_usage, const WellStateFullyImplicitBlackoil& wstate,
const std::vector<double>& b_perf, const PhaseUsage& phase_usage,
const std::vector<double>& rsmax_perf, const std::vector<double>& b_perf,
const std::vector<double>& rvmax_perf, const std::vector<double>& rsmax_perf,
const std::vector<double>& z_perf, const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens, const std::vector<double>& z_perf,
const double gravity) const std::vector<double>& surf_dens,
const double gravity)
{ {
// Verify that we have consistent input. // Verify that we have consistent input.
const int np = wells.number_of_phases; const int np = wells.number_of_phases;
@ -46,7 +47,7 @@ std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(co
if (surf_dens.size() != size_t(wells.number_of_phases)) { if (surf_dens.size() != size_t(wells.number_of_phases)) {
OPM_THROW(std::logic_error, "Inconsistent input: surf_dens vs. phase_usage."); OPM_THROW(std::logic_error, "Inconsistent input: surf_dens vs. phase_usage.");
} }
if (nperf*np != int(wstate.perfRates().size())) { if (nperf*np != int(wstate.perfPhaseRates().size())) {
OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate."); OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate.");
} }
if (nperf*np != int(b_perf.size())) { if (nperf*np != int(b_perf.size())) {
@ -90,7 +91,7 @@ std::vector<double> Opm::WellDensitySegmented::computeConnectionPressureDelta(co
q_out_perf[perf*np + phase] = q_out_perf[(perf+1)*np + phase]; q_out_perf[perf*np + phase] = q_out_perf[(perf+1)*np + phase];
} }
// Subtract outflow through perforation. // Subtract outflow through perforation.
q_out_perf[perf*np + phase] -= wstate.perfRates()[perf*np + phase]; q_out_perf[perf*np + phase] -= wstate.perfPhaseRates()[perf*np + phase];
} }
} }
} }

View File

@ -27,7 +27,7 @@ struct Wells;
namespace Opm namespace Opm
{ {
class WellState; class WellStateFullyImplicitBlackoil;
struct PhaseUsage; struct PhaseUsage;
@ -51,7 +51,7 @@ namespace Opm
/// \param[in] surf_dens surface densities for active components, size P /// \param[in] surf_dens surface densities for active components, size P
/// \param[in] gravity gravity acceleration constant /// \param[in] gravity gravity acceleration constant
static std::vector<double> computeConnectionPressureDelta(const Wells& wells, static std::vector<double> computeConnectionPressureDelta(const Wells& wells,
const WellState& wstate, const WellStateFullyImplicitBlackoil& wstate,
const PhaseUsage& phase_usage, const PhaseUsage& phase_usage,
const std::vector<double>& b_perf, const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf, const std::vector<double>& rsmax_perf,

View File

@ -0,0 +1,101 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/simulator/WellState.hpp>
#include <vector>
#include <cassert>
namespace Opm
{
/// The state of a set of wells, tailored for use by the fully
/// implicit blackoil simulator.
class WellStateFullyImplicitBlackoil
{
public:
/// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls
template <class State>
void init(const Wells* wells, const State& state)
{
if (wells == 0) {
return;
}
// We use the WellState::init() function to do bhp and well rates init.
// The alternative would be to copy that function wholesale.
basic_well_state_.init(wells, state);
// Initialize perfphaserates_, which must be done here.
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
const int nperf = wells->well_connpos[nw];
perfphaserates_.resize(nperf * np, 0.0);
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
const WellControls* ctrl = wells->ctrls[w];
if (well_controls_well_is_shut(ctrl)) {
// Shut well: perfphaserates_ are all zero.
} else {
// Open well: Initialize perfphaserates_ to well
// rates divided by the number of perforations.
const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
for (int p = 0; p < np; ++p) {
perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
}
}
}
}
}
/// One bhp pressure per well.
std::vector<double>& bhp() { return basic_well_state_.bhp(); }
const std::vector<double>& bhp() const { return basic_well_state_.bhp(); }
/// One rate per well and phase.
std::vector<double>& wellRates() { return basic_well_state_.wellRates(); }
const std::vector<double>& wellRates() const { return basic_well_state_.wellRates(); }
/// One rate per phase and well connection.
std::vector<double>& perfPhaseRates() { return perfphaserates_; }
const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
/// For interfacing with functions that take a WellState.
const WellState& basicWellState() const
{
return basic_well_state_;
}
private:
WellState basic_well_state_;
std::vector<double> perfphaserates_;
};
} // namespace Opm
#endif // OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED

View File

@ -28,7 +28,7 @@
#include <opm/autodiff/WellDensitySegmented.hpp> #include <opm/autodiff/WellDensitySegmented.hpp>
#include <opm/core/wells.h> #include <opm/core/wells.h>
#include <opm/core/simulator/WellState.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/core/props/BlackoilPhases.hpp> #include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/Units.hpp> #include <opm/core/utility/Units.hpp>
@ -65,8 +65,8 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas)
1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
1.0, 0.0, 0.0 }; 1.0, 0.0, 0.0 };
WellState wellstate; WellStateFullyImplicitBlackoil wellstate;
wellstate.perfRates() = rates; wellstate.perfPhaseRates() = rates;
PhaseUsage pu; PhaseUsage pu;
pu.num_phases = 3; pu.num_phases = 3;
pu.phase_used[0] = true; pu.phase_used[0] = true;