OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge from upstream.

Browse Source
This commit is contained in:
Bård Skaflestad 2012-05-09 18:56:09 +02:00
commit 8afc0571f9
15 changed files with 536 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

115
examples/spu_2p.cpp
View File

@ -297,6 +297,8 @@ main(int argc, char** argv)
boost::scoped_ptr<Opm::RockCompressibility> rock_comp;
Opm::SimulatorTimer simtimer;
Opm::TwophaseState state;
bool check_well_controls = false;
int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
@ -309,6 +311,8 @@ main(int argc, char** argv)
props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
// Wells init.
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
check_well_controls = param.getDefault("check_well_controls", false);
max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Timer init.
if (deck.hasField("TSTEP")) {
simtimer.init(deck);
@ -512,10 +516,13 @@ main(int argc, char** argv)
Opm::WellReport wellreport;
std::vector<double> well_bhp;
std::vector<double> well_perfrates;
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
if (wells->c_wells()) {
const int nw = wells->c_wells()->number_of_wells;
well_bhp.resize(nw, 0.0);
well_perfrates.resize(wells->c_wells()->well_connpos[nw], 0.0);
well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0);
wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates);
}
for (; !simtimer.done(); ++simtimer) {
@ -535,42 +542,82 @@ main(int argc, char** argv)
if (wells->c_wells()) {
Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), gravity[2], true, wdp);
}
pressure_timer.start();
if (rock_comp->isActive()) {
rc.resize(num_cells);
std::vector<double> initial_pressure = state.pressure();
std::vector<double> prev_pressure;
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
prev_pressure = state.pressure();
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
}
state.pressure() = initial_pressure;
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
state.pressure(), state.faceflux(), well_bhp, well_perfrates);
double max_change = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) {
break;
}
}
computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
} else {
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
if (check_well_controls) {
computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows);
}
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
if (check_well_controls) {
wells->applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
bool well_control_passed = !check_well_controls;
int well_control_iteration = 0;
do {
pressure_timer.start();
if (rock_comp->isActive()) {
rc.resize(num_cells);
std::vector<double> initial_pressure = state.pressure();
std::vector<double> initial_porevolume(num_cells);
computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume);
std::vector<double> pressure_increment(num_cells);
std::vector<double> prev_pressure;
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
}
computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
prev_pressure = state.pressure();
// compute pressure increment
psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc,
prev_pressure, initial_porevolume, simtimer.currentStepLength(),
pressure_increment);
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_perfrates, reorder_src);
if (!use_reorder) {
double max_change = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += pressure_increment[cell];
max_change = std::max(max_change, std::fabs(pressure_increment[cell]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) {
break;
}
}
psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
} else {
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
}
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
if (check_well_controls) {
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
fractional_flows,
well_perfrates,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_perfrates, reorder_src);
if (!use_reorder) {
clear_transport_source(tsrc);
for (int cell = 0; cell < num_cells; ++cell) {
if (reorder_src[cell] > 0.0) {

3
opm/core/InjectionSpecification.cpp
View File

@ -8,7 +8,8 @@ namespace Opm
surface_flow_max_rate_(-1e100),
reservoir_flow_max_rate_(-1e100),
BHP_limit_(-1e100),
reinjection_fraction_target_(-1e100),
reinjection_fraction_target_(1),
voidage_replacment_fraction_(1),
guide_rate_(1.0),
guide_rate_type_(NONE_GRT)
{

1
opm/core/InjectionSpecification.hpp
View File

@ -32,6 +32,7 @@ namespace Opm
double reservoir_flow_max_rate_;
double BHP_limit_;
double reinjection_fraction_target_;
double voidage_replacment_fraction_;
double guide_rate_;
GuideRateType guide_rate_type_;
};

2
opm/core/ProductionSpecification.hpp
View File

@ -11,7 +11,7 @@ namespace Opm
enum ControlMode
{
NONE, ORAT, WRAT, GRAT, LRAT, CRAT, RESV, PRBL, BHP, THP, GRUP, FLD
NONE = 0, ORAT = 1, WRAT=2, GRAT=3, LRAT=4, CRAT=5, RESV=6, PRBL=7, BHP=8, THP=9, GRUP=10, FLD=11
};
enum Procedure

18
opm/core/WellCollection.cpp
View File

@ -122,4 +122,22 @@ namespace Opm
roots_[i]->applyInjGroupControls();
}
}
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
void WellCollection::applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase)
{
for (size_t i = 0; i < roots_.size(); ++i) {
roots_[i]->applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase);
}
}
}

12
opm/core/WellCollection.hpp
View File

@ -88,6 +88,18 @@ namespace Opm
/// Applies all group controls (injection and production)
void applyGroupControls();
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase);
private:
// To account for the possibility of a forest
std::vector<std::tr1::shared_ptr<WellsGroupInterface> > roots_;

129
opm/core/WellsGroup.cpp
View File

@ -453,7 +453,7 @@ namespace Opm
// as that would check if we're under group control, something we're not.
const double children_guide_rate = children_[i]->productionGuideRate(true);
children_[i]->applyProdGroupControl(prod_mode,
(children_guide_rate / my_guide_rate) * getTarget(prod_mode),
(children_guide_rate / my_guide_rate) * getTarget(prod_mode),
false);
}
break;
@ -489,8 +489,9 @@ namespace Opm
}
return;
}
case InjectionSpecification::VREP:
case InjectionSpecification::REIN:
std::cout << "WARNING: Ignoring control type REIN" << std::endl;
std::cout << "Replacement keywords found, remember to call applyExplicitReinjectionControls." << std::endl;
return;
case InjectionSpecification::FLD:
case InjectionSpecification::NONE:
@ -528,7 +529,90 @@ namespace Opm
return sum;
}
/// Gets the total production flow of the given phase.
/// \param[in] phase_flows A vector containing rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] phase The phase for which to sum up.
double WellsGroup::getTotalProductionFlow(const std::vector<double>& phase_flows,
const BlackoilPhases::PhaseIndex phase)
{
double sum = 0.0;
for (size_t i = 0; i < children_.size(); ++i) {
sum += children_[i]->getTotalProductionFlow(phase_flows, phase);
}
return sum;
}
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
void WellsGroup::applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase)
{
if (injSpec().control_mode_ == InjectionSpecification::REIN) {
BlackoilPhases::PhaseIndex phase;
switch (injSpec().injector_type_) {
case InjectionSpecification::WATER:
phase = BlackoilPhases::Aqua;
break;
case InjectionSpecification::GAS:
phase = BlackoilPhases::Vapour;
break;
case InjectionSpecification::OIL:
phase = BlackoilPhases::Liquid;
break;
}
const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase);
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
// Note, we do _not_ want to call the applyProdGroupControl in this object,
// as that would check if we're under group control, something we're not.
const double children_guide_rate = children_[i]->injectionGuideRate(true);
#ifdef DIRTY_WELLCTRL_HACK
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
false);
#else
children_[i]->applyInjGroupControl(InjectionSpecification::RATE,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
false);
#endif
}
}
else if (injSpec().control_mode_ == InjectionSpecification::VREP) {
double total_produced = 0.0;
if (phaseUsage().phase_used[BlackoilPhases::Aqua]) {
total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Aqua);
}
if (phaseUsage().phase_used[BlackoilPhases::Liquid]) {
total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Liquid);
}
if (phaseUsage().phase_used[BlackoilPhases::Vapour]) {
total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour);
}
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
// Note, we do _not_ want to call the applyProdGroupControl in this object,
// as that would check if we're under group control, something we're not.
const double children_guide_rate = children_[i]->injectionGuideRate(true);
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_,
false);
}
}
}
// ============== WellNode members ============
@ -691,7 +775,7 @@ namespace Opm
{
// Not changing if we're not forced to change
if (!forced
&& (injSpec().control_mode_ != InjectionSpecification::GRUP || injSpec().control_mode_ != InjectionSpecification::NONE)) {
&& (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) {
return;
}
if (!wells_->type[self_index_] == INJECTOR) {
@ -728,13 +812,47 @@ namespace Opm
}
/// Gets the total production flow of the given phase.
/// \param[in] phase_flows A vector containing rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] phase The phase for which to sum up.
double WellNode::getTotalProductionFlow(const std::vector<double>& phase_flows,
const BlackoilPhases::PhaseIndex phase)
{
if (type() == INJECTOR) {
return 0.0;
}
return phase_flows[self_index_*phaseUsage().num_phases + phaseUsage().phase_pos[phase]];
}
WellType WellNode::type() const {
return wells_->type[self_index_];
}
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
void WellNode::applyExplicitReinjectionControls(const std::vector<double>&,
const std::vector<double>&)
{
// Do nothing at well level.
}
void WellNode::applyProdGroupControl(const ProductionSpecification::ControlMode control_mode,
const double target,
const bool forced)
{
// Not changing if we're not forced to change
if (!forced && (prodSpec().control_mode_ != ProductionSpecification::GRUP
|| prodSpec().control_mode_ != ProductionSpecification::NONE)) {
&& prodSpec().control_mode_ != ProductionSpecification::NONE)) {
std::cout << "Returning" << std::endl;
return;
}
if (!wells_->type[self_index_] == PRODUCER) {
@ -771,6 +889,7 @@ namespace Opm
distr[phase_pos[BlackoilPhases::Vapour]] = 1.0;
break;
case ProductionSpecification::LRAT:
std::cout << "applying rate" << std::endl;
wct = SURFACE_RATE;
if (!phase_used[BlackoilPhases::Liquid]) {
THROW("Oil phase not active and LRAT control specified.");
@ -981,6 +1100,8 @@ namespace Opm
injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
injection_specification.reservoir_flow_max_rate_ = line.resv_flow_max_rate_;
injection_specification.reinjection_fraction_target_ = line.reinjection_fraction_target_;
injection_specification.voidage_replacment_fraction_ = line.voidage_replacement_fraction_;
}
}
}

64
opm/core/WellsGroup.hpp
View File

@ -179,6 +179,27 @@ namespace Opm
/// wells under group control
virtual double injectionGuideRate(bool only_group) = 0;
/// Gets the total production flow of the given phase.
/// \param[in] phase_flows A vector containing rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] phase The phase for which to sum up.
virtual double getTotalProductionFlow(const std::vector<double>& phase_flows,
const BlackoilPhases::PhaseIndex phase) = 0;
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
virtual void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase) = 0;
protected:
/// Calculates the correct rate for the given ProductionSpecification::ControlMode
double rateByMode(const double* res_rates,
@ -258,6 +279,26 @@ namespace Opm
/// \param[in] only_group If true, will only accumelate guide rates for
/// wells under group control
virtual double injectionGuideRate(bool only_group);
/// Gets the total production flow of the given phase.
/// \param[in] phase_flows A vector containing rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] phase The phase for which to sum up.
virtual double getTotalProductionFlow(const std::vector<double>& phase_flows,
const BlackoilPhases::PhaseIndex phase);
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
virtual void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase);
private:
std::vector<std::tr1::shared_ptr<WellsGroupInterface> > children_;
@ -327,6 +368,29 @@ namespace Opm
/// \param[in] only_group If true, will only accumelate guide rates for
/// wells under group control
virtual double injectionGuideRate(bool only_group);
/// Gets the total production flow of the given phase.
/// \param[in] phase_flows A vector containing rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] phase The phase for which to sum up.
virtual double getTotalProductionFlow(const std::vector<double>& phase_flows,
const BlackoilPhases::PhaseIndex phase);
/// Returns the type of the well.
WellType type() const;
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
virtual void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase);
private:
Wells* wells_;

16
opm/core/WellsManager.cpp
View File

@ -702,4 +702,20 @@ namespace Opm
well_surfacerates_phase);
}
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
void WellsManager::applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase)
{
well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase);
}
} // namespace Opm

12
opm/core/WellsManager.hpp
View File

@ -88,6 +88,18 @@ namespace Opm
bool conditionsMet(const std::vector<double>& well_bhp,
const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase);
/// Applies explicit reinjection controls. This must be called at each timestep to be correct.
/// \param[in] well_reservoirrates_phase
/// A vector containing reservoir rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
/// \param[in] well_surfacerates_phase
/// A vector containing surface rates by phase for each well.
/// Is assumed to be ordered the same way as the related Wells-struct,
/// with all phase rates of a single well adjacent in the array.
void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& well_surfacerates_phase);
private:
// Disable copying and assignment.

5
opm/core/eclipse/SpecialEclipseFields.hpp
View File

@ -1001,12 +1001,14 @@ struct GconinjeLine
double surface_flow_max_rate_; // Surface flow rate target or upper limit
double resv_flow_max_rate_; // Reservoir flow rate target or upper limit
double reinjection_fraction_target_;
double voidage_replacement_fraction_;
// Default values
GconinjeLine() :
surface_flow_max_rate_(-1.0E20),
resv_flow_max_rate_(-1E20),
reinjection_fraction_target_(-1E20)
reinjection_fraction_target_(1),
voidage_replacement_fraction_(1)
{
}
};
@ -1048,6 +1050,7 @@ struct GCONINJE : public SpecialBase
gconinje_line.surface_flow_max_rate_ = double_data[0];
gconinje_line.resv_flow_max_rate_ = double_data[1];
gconinje_line.reinjection_fraction_target_ = double_data[2];
gconinje_line.voidage_replacement_fraction_ = double_data[3];
// HACK! Ignore any further items
if (num_read == num_to_read) {
ignoreSlashLine(is);

109
opm/core/pressure/IncompTpfa.cpp
View File

@ -26,6 +26,7 @@
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/newwells.h>
#include <string.h>
namespace Opm
{
@ -210,13 +211,16 @@ namespace Opm
F.bc = bcs;
F.totmob = &totmob[0];
F.wdp = &wdp[0];
int ok = true;
if (rock_comp.empty()) {
ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
ok = ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
} else {
ifs_tpfa_assemble_comprock(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
ok = ifs_tpfa_assemble_comprock(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
&porevol[0], &rock_comp[0], dt, &pressure[0], h_);
}
if (!ok) {
THROW("Failed assembling pressure system.");
}
linsolver_.solve(h_->A, h_->b, h_->x);
@ -233,9 +237,108 @@ namespace Opm
soln.well_flux = &well_rate[0];
soln.well_press = &well_bhp[0];
}
ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
}
void IncompTpfa::solveIncrement(const std::vector<double>& totmob,
const std::vector<double>& omega,
const std::vector<double>& src,
const std::vector<double>& wdp,
const FlowBoundaryConditions* bcs,
const std::vector<double>& porevol,
const std::vector<double>& rock_comp,
const std::vector<double>& prev_pressure,
const std::vector<double>& initial_porevol,
const double dt,
std::vector<double>& pressure_increment)
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
if (!omega.empty()) {
if (gpress_.empty()) {
THROW("Nozero omega argument given, but gravity was null in constructor.");
}
mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
&omega[0],
&gpress_[0], &gpress_omegaweighted_[0]);
} else {
if (!gpress_.empty()) {
THROW("Empty omega argument given, but gravity was non-null in constructor.");
}
}
ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
if (! src.empty()) { F.src = &src[0]; }
F.bc = bcs;
F.totmob = &totmob[0];
F.wdp = &wdp[0];
if (rock_comp.empty()) {
ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
} else {
ifs_tpfa_assemble_comprock_increment(gg, &F, &trans_[0], &gpress_omegaweighted_[0],
&porevol[0], &rock_comp[0], dt, &prev_pressure[0],
&initial_porevol[0], h_);
}
linsolver_.solve(h_->A, h_->b, &pressure_increment[0]);
}
void IncompTpfa::computeFaceFlux(const std::vector<double>& totmob,
const std::vector<double>& omega,
const std::vector<double>& src,
const std::vector<double>& wdp,
const FlowBoundaryConditions* bcs,
std::vector<double>& pressure,
std::vector<double>& faceflux,
std::vector<double>& well_bhp,
std::vector<double>& well_rate) {
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
if (!omega.empty()) {
if (gpress_.empty()) {
THROW("Nozero omega argument given, but gravity was null in constructor.");
}
mim_ip_density_update(gg->number_of_cells, gg->cell_facepos,
&omega[0],
&gpress_[0], &gpress_omegaweighted_[0]);
} else {
if (!gpress_.empty()) {
THROW("Empty omega argument given, but gravity was non-null in constructor.");
}
}
ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
if (! src.empty()) { F.src = &src[0]; }
F.bc = bcs;
F.totmob = &totmob[0];
F.wdp = &wdp[0];
ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
faceflux.resize(grid_.number_of_faces);
ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL };
soln.cell_press = &pressure[0];
soln.face_flux = &faceflux[0];
if(wells_ != NULL) {
well_bhp.resize(wells_->number_of_wells);
well_rate.resize(wells_->well_connpos[ wells_->number_of_wells ]);
soln.well_flux = &well_rate[0];
soln.well_press = &well_bhp[0];
}
memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x));
ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
}
} // namespace Opm

22
opm/core/pressure/IncompTpfa.hpp
View File

@ -121,6 +121,28 @@ namespace Opm
std::vector<double>& well_bhp,
std::vector<double>& well_rate);
void solveIncrement(const std::vector<double>& totmob,
const std::vector<double>& omega,
const std::vector<double>& src,
const std::vector<double>& wdp,
const FlowBoundaryConditions* bcs,
const std::vector<double>& porevol,
const std::vector<double>& rock_comp,
const std::vector<double>& prev_pressure,
const std::vector<double>& initial_porevol,
const double dt,
std::vector<double>& pressure_increment);
void computeFaceFlux(const std::vector<double>& totmob,
const std::vector<double>& omega,
const std::vector<double>& src,
const std::vector<double>& wdp,
const FlowBoundaryConditions* bcs,
std::vector<double>& pressure,
std::vector<double>& faceflux,
std::vector<double>& well_bhp,
std::vector<double>& well_rate);
/// Expose read-only reference to internal half-transmissibility.
const ::std::vector<double>& getHalfTrans() const { return htrans_; }

57
opm/core/pressure/tpfa/ifs_tpfa.c
View File

@ -11,6 +11,16 @@
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v) {
size_t i,j;
for (j = 0; j < A->m; ++j) {
v[j] = 0;
for (i = (size_t) (A->ia[j]); i < (size_t) (A->ia[j+1]); ++i) {
v[j] += A->sa[i]*u[A->ja[i]];
}
}
}
struct ifs_tpfa_impl {
double *fgrav; /* Accumulated grav contrib/face */
@ -287,7 +297,7 @@ assemble_shut_well(int nc, int w,
size_t jw;
double trans;
/* The equation added for a shut well w is
/* The equation added for a shut well w is
* \sum_{c~w} T_{w,c} * mt_c p_w = 0.0
* where c~w indicates the cell perforated by w, T are production
* indices, mt is total mobility and p_w is the bottom-hole
@ -371,7 +381,7 @@ assemble_well_contrib(int nc ,
/* We cannot handle this case, since we do
* not have access to formation volume factors,
* needed to convert between reservoir and
* surface rates
* surface rates
*/
fprintf(stderr, "ifs_tpfa cannot handle SURFACE_RATE well controls.\n");
*ok = 0;
@ -725,6 +735,49 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid *G ,
return ok;
}
/* ---------------------------------------------------------------------- */
int
ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G ,
const struct ifs_tpfa_forces *F ,
const double *trans ,
const double *gpress ,
const double *porevol ,
const double *rock_comp,
const double dt ,
const double *prev_pressure,
const double *initial_porevolume,
struct ifs_tpfa_data *h )
/* ---------------------------------------------------------------------- */
{
int c;
size_t j;
int system_singular, ok;
double *v;
assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok);
v = malloc(h->A->m * sizeof *v);
mult_csr_matrix(h->A, prev_pressure, v);
/* We want to solve a Newton step for the residual
* (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible
*
*/
if (ok) {
for (c = 0; c < G->number_of_cells; ++c) {
j = csrmatrix_elm_index(c, c, h->A);
h->A->sa[j] += porevol[c] * rock_comp[c] / dt;
h->b[c] += -(porevol[c] - initial_porevolume[c])/dt - v[c];
}
}
free(v);
return ok;
}
/* ---------------------------------------------------------------------- */
void

17
opm/core/pressure/tpfa/ifs_tpfa.h
View File

@ -61,6 +61,10 @@ struct ifs_tpfa_data *
ifs_tpfa_construct(struct UnstructuredGrid *G,
struct Wells *W);
void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v);
int
ifs_tpfa_assemble(struct UnstructuredGrid *G ,
const struct ifs_tpfa_forces *F ,
@ -78,6 +82,19 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid *G ,
const double dt ,
const double *pressure ,
struct ifs_tpfa_data *h );
int
ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G ,
const struct ifs_tpfa_forces *F ,
const double *trans ,
const double *gpress ,
const double *porevol ,
const double *rock_comp,
const double dt ,
const double *prev_pressure ,
const double *initial_porevolume,
struct ifs_tpfa_data *h );
void
ifs_tpfa_press_flux(struct UnstructuredGrid *G ,
const struct ifs_tpfa_forces *F ,