OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

fixing reviewing comments from PR 1279.

Browse Source
This commit is contained in:
Kai Bao 2017-10-16 16:46:28 +02:00
parent 4efaf64cf7
commit ba8eb708d1
14 changed files with 113 additions and 172 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -391,7 +391,7 @@ namespace Opm {
}
catch ( const Dune::FMatrixError& e )
{
OPM_THROW(Opm::NumericalProblem,"The following FMatrixError got caught during well modeling \n" << e.what());
OPM_THROW(Opm::NumericalProblem,"Error encounted when solving well equations");
}
return report;
@ -425,12 +425,12 @@ namespace Opm {
saturationsNew[FluidSystem::waterPhaseIdx] = priVarsNew[Indices::waterSaturationIdx];
oilSaturationNew -= saturationsNew[FluidSystem::waterPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsNew.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
saturationsNew[FluidSystem::gasPhaseIdx] = priVarsNew[Indices::compositionSwitchIdx];
oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
oilSaturationNew -= saturationsNew[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsNew[FluidSystem::oilPhaseIdx] = oilSaturationNew;
}
@ -446,16 +446,16 @@ namespace Opm {
saturationsOld[FluidSystem::waterPhaseIdx] = priVarsOld[Indices::waterSaturationIdx];
oilSaturationOld -= saturationsOld[FluidSystem::waterPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && priVarsOld.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
saturationsOld[FluidSystem::gasPhaseIdx] = priVarsOld[Indices::compositionSwitchIdx];
oilSaturationOld -= saturationsOld[FluidSystem::gasPhaseIdx];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
saturationsOld[FluidSystem::oilPhaseIdx] = oilSaturationOld;
}
Scalar tmp = pressureNew - pressureOld;
resultDelta += tmp*tmp;
resultDenom += pressureNew*pressureNew;
@ -680,7 +680,7 @@ namespace Opm {
maxVal = std::max(std::abs(dsw),maxVal);
maxVal = std::max(std::abs(dsg),maxVal);
maxVal = std::max(std::abs(dso),maxVal);
maxVal = std::max(std::abs(dss),maxVal);
maxVal = std::max(std::abs(dss),maxVal);
double satScaleFactor = 1.0;
if (maxVal > dsMax()) {

14
opm/autodiff/BlackoilModelParameters.cpp
View File

@ -50,11 +50,13 @@ namespace Opm
tolerance_cnv_ = param.getDefault("tolerance_cnv", tolerance_cnv_);
tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
tolerance_well_control_ = param.getDefault("tolerance_well_control", tolerance_well_control_);
tolerance_pressure_ms_wells_ = param.getDefault("tolerance_pressure_ms_wells", tolerance_pressure_ms_wells_);
max_welleq_iter_ = param.getDefault("max_welleq_iter", max_welleq_iter_);
max_pressure_change_ms_wells_ = param.getDefault("max_pressure_change_ms_wells", max_pressure_change_ms_wells_);
use_inner_iterations_ms_wells_ = param.getDefault("use_inner_iterations_ms_wells", use_inner_iterations_ms_wells_);
max_inner_iter_ms_wells_ = param.getDefault("max_inner_iter_ms_wells", max_inner_iter_ms_wells_);
if (use_multisegment_well_) {
tolerance_pressure_ms_wells_ = param.getDefault("tolerance_pressure_ms_wells", tolerance_pressure_ms_wells_);
max_pressure_change_ms_wells_ = param.getDefault("max_pressure_change_ms_wells", max_pressure_change_ms_wells_);
use_inner_iterations_ms_wells_ = param.getDefault("use_inner_iterations_ms_wells", use_inner_iterations_ms_wells_);
max_inner_iter_ms_wells_ = param.getDefault("max_inner_iter_ms_wells", max_inner_iter_ms_wells_);
}
maxSinglePrecisionTimeStep_ = unit::convert::from(
param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day );
max_strict_iter_ = param.getDefault("max_strict_iter",8);
@ -81,9 +83,9 @@ namespace Opm
tolerance_cnv_ = 1.0e-2;
tolerance_wells_ = 1.0e-4;
tolerance_well_control_ = 1.0e-7;
tolerance_pressure_ms_wells_ = 1000.0; // 0.01 bar
tolerance_pressure_ms_wells_ = unit::convert::from(0.01, unit::barsa); // 0.01 bar
max_welleq_iter_ = 15;
max_pressure_change_ms_wells_ = 200000.; // 2.0 bar
max_pressure_change_ms_wells_ = unit::convert::from(2.0, unit::barsa); // 2.0 bar
use_inner_iterations_ms_wells_ = true;
max_inner_iter_ms_wells_ = 10;
maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day );

12
opm/autodiff/BlackoilModelParameters.hpp
View File

@ -81,14 +81,14 @@ namespace Opm
/// Try to detect oscillation or stagnation.
bool use_update_stabilization_;
/// Whehter to use MultisegmentWell to handle multisegment wells
// it is something temporary before the multisegment well model is considered to be
// well developed and tested.
// if it is false, we will handle multisegment wells as standard wells, which will be
// the default behavoir for the moment. Later, we might it to be true by default if necessary
/// Whether to use MultisegmentWell to handle multisegment wells
/// it is something temporary before the multisegment well model is considered to be
/// well developed and tested.
/// if it is false, we will handle multisegment wells as standard wells, which will be
/// the default behavoir for the moment. Later, we might set it to be true by default if necessary
bool use_multisegment_well_;
// The file name of the deck
/// The file name of the deck
std::string deck_file_name_;
/// Construct from user parameters or defaults.

8
opm/autodiff/BlackoilWellModel.hpp
View File

@ -138,7 +138,7 @@ namespace Opm {
/// return true if wells are available on this process
bool localWellsActive() const;
bool getWellConvergence(Simulator& ebosSimulator,
bool getWellConvergence(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg) const;
/// upate the dynamic lists related to economic limits
@ -163,6 +163,8 @@ namespace Opm {
const int number_of_phases_;
const ModelParameters& param_;
using WellInterfacePtr = std::unique_ptr<WellInterface<TypeTag> >;
// a vector of all the wells.
// eventually, the wells_ above should be gone.
@ -176,12 +178,12 @@ namespace Opm {
static std::vector<WellInterfacePtr > createWellContainer(const Wells* wells,
const std::vector<const Well*>& wells_ecl,
const bool use_multisegment_well,
const int time_step);
const int time_step,
const ModelParameters& param);
// Well collection is used to enforce the group control
WellCollection* well_collection_;
ModelParameters param_;
bool terminal_output_;
bool has_solvent_;
bool has_polymer_;

21
opm/autodiff/BlackoilWellModel_impl.hpp
View File

@ -18,9 +18,9 @@ namespace Opm {
, wells_ecl_(wells_ecl)
, number_of_wells_(wells_arg ? (wells_arg->number_of_wells) : 0)
, number_of_phases_(wells_arg ? (wells_arg->number_of_phases) : 0) // TODO: not sure if it is proper for this way
, well_container_(createWellContainer(wells_arg, wells_ecl, param.use_multisegment_well_, current_timeIdx) )
, well_collection_(well_collection)
, param_(param)
, well_container_(createWellContainer(wells_arg, wells_ecl, param.use_multisegment_well_, current_timeIdx, param) )
, well_collection_(well_collection)
, terminal_output_(terminal_output)
, has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
, has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
@ -118,7 +118,8 @@ namespace Opm {
createWellContainer(const Wells* wells,
const std::vector< const Well* >& wells_ecl,
const bool use_multisegment_well,
const int time_step)
const int time_step,
const ModelParameters& param)
{
std::vector<WellInterfacePtr> well_container;
@ -149,9 +150,9 @@ namespace Opm {
const Well* well_ecl = wells_ecl[index_well];
if ( !well_ecl->isMultiSegment(time_step) || !use_multisegment_well) {
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells) );
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells, param) );
} else {
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells) );
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells, param) );
}
}
}
@ -210,7 +211,7 @@ namespace Opm {
bool only_wells) const
{
for (int w = 0; w < number_of_wells_; ++w) {
well_container_[w]->assembleWellEq(ebosSimulator, param_, dt, well_state, only_wells);
well_container_[w]->assembleWellEq(ebosSimulator, dt, well_state, only_wells);
}
}
@ -289,7 +290,7 @@ namespace Opm {
recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const
{
for (auto& well : well_container_) {
well->recoverWellSolutionAndUpdateWellState(x, param_, well_state);
well->recoverWellSolutionAndUpdateWellState(x, well_state);
}
}
@ -433,7 +434,7 @@ namespace Opm {
if( localWellsActive() )
{
for (auto& well : well_container_) {
well->solveEqAndUpdateWellState(param_, well_state);
well->solveEqAndUpdateWellState(well_state);
}
}
// updateWellControls uses communication
@ -477,13 +478,13 @@ namespace Opm {
template<typename TypeTag>
bool
BlackoilWellModel<TypeTag>::
getWellConvergence(Simulator& ebosSimulator,
getWellConvergence(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg) const
{
ConvergenceReport report;
for (const auto& well : well_container_) {
report += well->getWellConvergence(ebosSimulator, B_avg, param_);
report += well->getWellConvergence(B_avg);
}
// checking NaN residuals

4
opm/autodiff/MSWellHelpers.hpp
View File

@ -111,7 +111,7 @@ namespace mswellhelpers
static double haalandFormular(const double re, const double diameter, const double roughness)
inline double haalandFormular(const double re, const double diameter, const double roughness)
{
const double value = -3.6 * std::log10(6.9 / re + std::pow(roughness / (3.7 * diameter), 10. / 9.) );
@ -125,7 +125,7 @@ namespace mswellhelpers
static double calculateFrictionFactor(const double area, const double diameter,
inline double calculateFrictionFactor(const double area, const double diameter,
const double w, const double roughness, const double mu)
{

19
opm/autodiff/MultisegmentWell.hpp
View File

@ -34,7 +34,6 @@ namespace Opm
public:
typedef WellInterface<TypeTag> Base;
// TODO: the WellState does not have any information related to segments
using typename Base::WellState;
using typename Base::Simulator;
using typename Base::IntensiveQuantities;
@ -50,8 +49,7 @@ namespace Opm
using Base::has_polymer;
// TODO: for now, not considering the polymer, solvent and so on to simplify the development process.
// TODO: should I begin with the old primary variable or the new fraction based variable systems?
// Let us begin with the new one
// TODO: we need to have order for the primary variables and also the order for the well equations.
// sometimes, they are similar, while sometimes, they can have very different forms.
@ -96,7 +94,7 @@ namespace Opm
// TODO: for now, we only use one type to save some implementation efforts, while improve later.
typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
MultisegmentWell(const Well* well, const int time_step, const Wells* wells);
MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
@ -108,7 +106,6 @@ namespace Opm
virtual void initPrimaryVariablesEvaluation() const;
virtual void assembleWellEq(Simulator& ebosSimulator,
const ModelParameters& param,
const double dt,
WellState& well_state,
bool only_wells);
@ -119,9 +116,7 @@ namespace Opm
WellState& well_state) const;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator,
const std::vector<double>& B_avg,
const ModelParameters& param) const;
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const;
@ -130,7 +125,7 @@ namespace Opm
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param,
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state) const;
/// computing the well potentials for group control
@ -140,8 +135,7 @@ namespace Opm
virtual void updatePrimaryVariables(const WellState& well_state) const;
virtual void solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state); // const?
virtual void solveEqAndUpdateWellState(WellState& well_state); // const?
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state); // should be const?
@ -173,6 +167,7 @@ namespace Opm
// TODO: the current implementation really relies on the order of the
// perforation does not change from the parser to Wells structure.
using Base::well_cells_;
using Base::param_;
using Base::well_index_;
using Base::well_type_;
using Base::first_perf_;
@ -264,7 +259,6 @@ namespace Opm
// updating the well_state based on well solution dwells
void updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
const bool inner_iteration,
WellState& well_state) const;
@ -340,7 +334,6 @@ namespace Opm
// TODO: try to make ebosSimulator const, as it should be
void iterateWellEquations(Simulator& ebosSimulator,
const ModelParameters& param,
const double dt,
WellState& well_state);

90
opm/autodiff/MultisegmentWell_impl.hpp
View File

@ -27,8 +27,8 @@ namespace Opm
template <typename TypeTag>
MultisegmentWell<TypeTag>::
MultisegmentWell(const Well* well, const int time_step, const Wells* wells)
: Base(well, time_step, wells)
MultisegmentWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
: Base(well, time_step, wells, param)
, segment_perforations_(numberOfSegments())
, segment_inlets_(numberOfSegments())
, cell_perforation_depth_diffs_(number_of_perforations_, 0.0)
@ -228,15 +228,14 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
assembleWellEq(Simulator& ebosSimulator,
const ModelParameters& param,
const double dt,
WellState& well_state,
bool only_wells)
{
const bool use_inner_iterations = param.use_inner_iterations_ms_wells_;
const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
if (use_inner_iterations) {
iterateWellEquations(ebosSimulator, param, dt, well_state);
iterateWellEquations(ebosSimulator, dt, well_state);
}
assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells);
@ -407,18 +406,16 @@ namespace Opm
template <typename TypeTag>
typename MultisegmentWell<TypeTag>::ConvergenceReport
MultisegmentWell<TypeTag>::
getWellConvergence(const Simulator& /* ebosSimulator */,
const std::vector<double>& B_avg,
const ModelParameters& param) const
getWellConvergence(const std::vector<double>& B_avg) const
{
// assert((int(B_avg.size()) == numComponents()) || has_polymer);
assert( (int(B_avg.size()) == numComponents()) );
// checking if any residual is NaN or too large. The two large one is only handled for the well flux
std::vector<std::vector<double>> residual(numberOfSegments(), std::vector<double>(numWellEq, 0.0));
std::vector<std::vector<double>> abs_residual(numberOfSegments(), std::vector<double>(numWellEq, 0.0));
for (int seg = 0; seg < numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]);
abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]);
}
}
@ -429,14 +426,14 @@ namespace Opm
for (int seg = 0; seg < numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
if (eq_idx < numComponents()) { // phase or component mass equations
const double flux_residual = B_avg[eq_idx] * residual[seg][eq_idx];
const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
// TODO: the report can not handle the segment number yet.
if (std::isnan(flux_residual)) {
report.nan_residual_found = true;
const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
report.nan_residual_wells.push_back(problem_well);
} else if (flux_residual > param.max_residual_allowed_) {
} else if (flux_residual > param_.max_residual_allowed_) {
report.too_large_residual_found = true;
const auto& phase_name = FluidSystem::phaseName(flowPhaseToEbosPhaseIdx(eq_idx));
const typename ConvergenceReport::ProblemWell problem_well = {name(), phase_name};
@ -449,7 +446,7 @@ namespace Opm
} else { // pressure equation
// TODO: we should distinguish the rate control equations, bhp control equations
// and the oridnary pressure equations
const double pressure_residal = residual[seg][eq_idx];
const double pressure_residal = abs_residual[seg][eq_idx];
const std::string eq_name("Pressure");
if (std::isnan(pressure_residal)) {
report.nan_residual_found = true;
@ -472,10 +469,10 @@ namespace Opm
// check convergence for flux residuals
for ( int comp_idx = 0; comp_idx < numComponents(); ++comp_idx)
{
report.converged = report.converged && (maximum_residual[comp_idx] < param.tolerance_wells_);
report.converged = report.converged && (maximum_residual[comp_idx] < param_.tolerance_wells_);
}
report.converged = report.converged && (maximum_residual[SPres] < param.tolerance_pressure_ms_wells_);
report.converged = report.converged && (maximum_residual[SPres] < param_.tolerance_pressure_ms_wells_);
} else { // abnormal values found and no need to check the convergence
report.converged = false;
}
@ -526,12 +523,11 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
recoverWellSolutionAndUpdateWellState(const BVector& x,
const ModelParameters& param,
WellState& well_state) const
{
BVectorWell xw(1);
recoverSolutionWell(x, xw);
updateWellState(xw, param, false, well_state);
updateWellState(xw, false, well_state);
}
@ -557,8 +553,6 @@ namespace Opm
MultisegmentWell<TypeTag>::
updatePrimaryVariables(const WellState& well_state) const
{
// TODO: not handling solvent or polymer for now.
// TODO: to test using rate conversion coefficients to see if it will be better than
// this default one
@ -603,7 +597,6 @@ namespace Opm
if (active()[Gas]) {
if (distr[pu.phase_pos[Gas]] > 0.0) {
// TODO: not handling solvent here yet
primary_variables_[seg][GFrac] = 1.0;
} else {
primary_variables_[seg][GFrac] = 0.0;
@ -645,14 +638,13 @@ namespace Opm
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state)
solveEqAndUpdateWellState(WellState& well_state)
{
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_);
updateWellState(dx_well, param, false, well_state);
updateWellState(dx_well, false, well_state);
}
@ -736,16 +728,15 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
const bool inner_iteration,
WellState& well_state) const
{
const bool use_inner_iterations = param.use_inner_iterations_ms_wells_;
const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_;
const double relaxation_factor = (use_inner_iterations && inner_iteration) ? 0.2 : 1.0;
const double dFLimit = param.dwell_fraction_max_;
const double max_pressure_change = param.max_pressure_change_ms_wells_;
const double dFLimit = param_.dwell_fraction_max_;
const double max_pressure_change = param_.max_pressure_change_ms_wells_;
const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
for (int seg = 0; seg < numberOfSegments(); ++seg) {
@ -776,8 +767,6 @@ namespace Opm
primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal];
}
// TODO: not handling solvent related for now
}
updateWellStateFromPrimaryVariables(well_state);
@ -888,11 +877,6 @@ namespace Opm
return primary_variables_evaluation_[seg][GFrac];
}
// TODO: not handling solvent for now
// if (has_solvent && compIdx == contiSolventEqIdx) {
// return primary_variables_evaluation_[seg][SFrac];
// }
// Oil fraction
EvalWell oil_fraction = 1.0;
if (active()[Water]) {
@ -985,11 +969,6 @@ namespace Opm
b_perfcells[phase] = extendEval(fs.invB(phase_idx_ebos));
}
// TODO: not handling solvent for now
// if (has_solvent) {
// b_perfcells[contiSolventEqIdx] = extendEval(intQuants.solventInverseFormationVolumeFactor());
// }
// pressure difference between the segment and the perforation
const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
@ -1044,11 +1023,6 @@ namespace Opm
volume_ratio += cmix_s[watpos] / b_perfcells[watpos];
}
// TODO: not handling
// if (has_solvent) {
// volumeRatio += cmix_s[contiSolventEqIdx] / b_perfcells_dense[contiSolventEqIdx];
// }
if (active()[Oil] && active()[Gas]) {
const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas];
@ -1256,8 +1230,6 @@ namespace Opm
segment_viscosities_[seg] += visc[p] * phase_fraction;
}
// TODO: not handling solvent for now.
EvalWell density(0.0);
for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) {
density += surf_dens[comp_idx] * mix_s[comp_idx];
@ -1356,24 +1328,7 @@ namespace Opm
mob[phase] = extendEval(relativePerms[ebosPhaseIdx] / intQuants.fluidState().viscosity(ebosPhaseIdx));
}
// this may not work if viscosity and relperms has been modified?
// if (has_solvent) {
// OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent");
// }
}
// modify the water mobility if polymer is present
// if (has_polymer) {
// assume fully mixture for wells.
// EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
// if (well_type_ == INJECTOR) {
// const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex());
// mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) );
// }
// TODO: not sure if we should handle shear calculation with MS well
// }
}
@ -1608,8 +1563,6 @@ namespace Opm
fractions[oil_pos] -= fractions[gas_pos];
}
// TODO: not handling solvent related
if (active()[Water]) {
const int water_pos = pu.phase_pos[Water];
if (fractions[water_pos] < 0.0) {
@ -1777,7 +1730,6 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
iterateWellEquations(Simulator& ebosSimulator,
const ModelParameters& param,
const double dt,
WellState& well_state)
{
@ -1786,7 +1738,7 @@ namespace Opm
// if converged, we can update the well_state.
// the function updateWellState() should have a flag to show
// if we will update the well state.
const int max_iter_number = param.max_inner_iter_ms_wells_;
const int max_iter_number = param_.max_inner_iter_ms_wells_;
int it = 0;
for (; it < max_iter_number; ++it) {
@ -1801,13 +1753,13 @@ namespace Opm
// const std::vector<double> B {0.8, 0.8, 0.008};
const std::vector<double> B {0.5, 0.5, 0.005};
const ConvergenceReport report = getWellConvergence(ebosSimulator, B, param);
const ConvergenceReport report = getWellConvergence(B);
if (report.converged) {
break;
}
updateWellState(dx_well, param, true, well_state);
updateWellState(dx_well, true, well_state);
initPrimaryVariablesEvaluation();
}

14
opm/autodiff/StandardWell.hpp
View File

@ -99,7 +99,7 @@ namespace Opm
static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
StandardWell(const Well* well, const int time_step, const Wells* wells);
StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param);
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
@ -111,7 +111,6 @@ namespace Opm
virtual void initPrimaryVariablesEvaluation() const;
virtual void assembleWellEq(Simulator& ebosSimulator,
const ModelParameters& param,
const double dt,
WellState& well_state,
bool only_wells);
@ -122,9 +121,7 @@ namespace Opm
WellState& xw) const;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator,
const std::vector<double>& B_avg,
const ModelParameters& param) const;
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const;
@ -133,7 +130,7 @@ namespace Opm
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param,
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state) const;
/// computing the well potentials for group control
@ -143,8 +140,7 @@ namespace Opm
virtual void updatePrimaryVariables(const WellState& well_state) const;
virtual void solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state);
virtual void solveEqAndUpdateWellState(WellState& well_state);
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state); // should be const?
@ -165,6 +161,7 @@ namespace Opm
// protected member variables from the Base class
using Base::vfp_properties_;
using Base::gravity_;
using Base::param_;
using Base::well_efficiency_factor_;
using Base::first_perf_;
using Base::ref_depth_;
@ -234,7 +231,6 @@ namespace Opm
// updating the well_state based on well solution dwells
void updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
WellState& well_state) const;
// calculate the properties for the well connections

26
opm/autodiff/StandardWell_impl.hpp
View File

@ -24,8 +24,8 @@ namespace Opm
{
template<typename TypeTag>
StandardWell<TypeTag>::
StandardWell(const Well* well, const int time_step, const Wells* wells)
: Base(well, time_step, wells)
StandardWell(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
: Base(well, time_step, wells, param)
, perf_densities_(number_of_perforations_)
, perf_pressure_diffs_(number_of_perforations_)
, primary_variables_(numWellEq, 0.0)
@ -514,7 +514,6 @@ namespace Opm
void
StandardWell<TypeTag>::
assembleWellEq(Simulator& ebosSimulator,
const ModelParameters& /* param */,
const double dt,
WellState& well_state,
bool only_wells)
@ -765,12 +764,11 @@ namespace Opm
void
StandardWell<TypeTag>::
updateWellState(const BVectorWell& dwells,
const BlackoilModelParameters& param,
WellState& well_state) const
{
const int np = number_of_phases_;
const double dBHPLimit = param.dbhp_max_rel_;
const double dFLimit = param.dwell_fraction_max_;
const double dBHPLimit = param_.dbhp_max_rel_;
const double dFLimit = param_.dwell_fraction_max_;
const auto pu = phaseUsage();
const std::vector<double> xvar_well_old = primary_variables_;
@ -1365,9 +1363,7 @@ namespace Opm
template<typename TypeTag>
typename StandardWell<TypeTag>::ConvergenceReport
StandardWell<TypeTag>::
getWellConvergence(const Simulator& /* ebosSimulator */,
const std::vector<double>& B_avg,
const ModelParameters& param) const
getWellConvergence(const std::vector<double>& B_avg) const
{
typedef double Scalar;
typedef std::vector< Scalar > Vector;
@ -1379,8 +1375,8 @@ namespace Opm
// For the polymer case, there is one more mass balance equations of reservoir than wells
assert((int(B_avg.size()) == numComp) || has_polymer);
const double tol_wells = param.tolerance_wells_;
const double maxResidualAllowed = param.max_residual_allowed_;
const double tol_wells = param_.tolerance_wells_;
const double maxResidualAllowed = param_.max_residual_allowed_;
// TODO: it should be the number of numWellEq
// using numComp here for flow_ebos running 2p case.
@ -1492,15 +1488,14 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state)
solveEqAndUpdateWellState(WellState& well_state)
{
// We assemble the well equations, then we check the convergence,
// which is why we do not put the assembleWellEq here.
BVectorWell dx_well(1);
invDuneD_.mv(resWell_, dx_well);
updateWellState(dx_well, param, well_state);
updateWellState(dx_well, well_state);
}
@ -1597,12 +1592,11 @@ namespace Opm
void
StandardWell<TypeTag>::
recoverWellSolutionAndUpdateWellState(const BVector& x,
const ModelParameters& param,
WellState& well_state) const
{
BVectorWell xw(1);
recoverSolutionWell(x, xw);
updateWellState(xw, param, well_state);
updateWellState(xw, well_state);
}

18
opm/autodiff/WellInterface.hpp
View File

@ -86,7 +86,7 @@ namespace Opm
static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
/// Constructor
WellInterface(const Well* well, const int time_step, const Wells* wells);
WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param);
/// Virutal destructor
virtual ~WellInterface() {}
@ -144,15 +144,11 @@ namespace Opm
}
};
virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator,
const std::vector<double>& B_avg,
const ModelParameters& param) const = 0;
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const = 0;
virtual void solveEqAndUpdateWellState(const ModelParameters& param,
WellState& well_state) = 0;
virtual void solveEqAndUpdateWellState(WellState& well_state) = 0;
virtual void assembleWellEq(Simulator& ebosSimulator,
const ModelParameters& param,
const double dt,
WellState& well_state,
bool only_wells) = 0;
@ -166,7 +162,7 @@ namespace Opm
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x, const ModelParameters& param,
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state) const = 0;
/// Ax = Ax - C D^-1 B x
@ -203,6 +199,9 @@ namespace Opm
// the index of well in Wells struct
int index_of_well_;
// simulation parameters
const ModelParameters& param_;
// well type
// INJECTOR or PRODUCER
enum WellType well_type_;
@ -215,21 +214,18 @@ namespace Opm
std::vector<double> comp_frac_;
// controls for this well
// TODO: later will check whehter to let it stay with pointer
struct WellControls* well_controls_;
// number of the perforations for this well
int number_of_perforations_;
// record the index of the first perforation
// TODO: it might not be needed if we refactor WellState to be a vector
// of states of individual well.
int first_perf_;
// well index for each perforation
std::vector<double> well_index_;
// TODO: it might should go to StandardWell
// depth for each perforation
std::vector<double> perf_depth_;

3
opm/autodiff/WellInterface_impl.hpp
View File

@ -25,9 +25,10 @@ namespace Opm
template<typename TypeTag>
WellInterface<TypeTag>::
WellInterface(const Well* well, const int time_step, const Wells* wells)
WellInterface(const Well* well, const int time_step, const Wells* wells, const ModelParameters& param)
: well_ecl_(well)
, current_step_(time_step)
, param_(param)
{
if (!well) {
OPM_THROW(std::invalid_argument, "Null pointer of Well is used to construct WellInterface");

29
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@ -311,7 +311,7 @@ namespace Opm
segrates_.clear();
segrates_.reserve(nw * numPhases());
nseg_ = 0.;
nseg_ = 0;
// in the init function, the well rates and perforation rates have been initialized or copied from prevState
// what we do here, is to set the segment rates and perforation rates
for (int w = 0; w < nw; ++w) {
@ -340,14 +340,15 @@ namespace Opm
}
} else { // it is a multi-segment well
const SegmentSet& segment_set = well_ecl->getSegmentSet(time_step);
// assuming the oder of the perforations in well_ecl is the same with Wells
// assuming the order of the perforations in well_ecl is the same with Wells
const CompletionSet& completion_set = well_ecl->getCompletions(time_step);
const int nseg = segment_set.numberSegment();
// number of segment for this single well
const int well_nseg = segment_set.numberSegment();
const int nperf = completion_set.size();
nseg_ += nseg;
nseg_ += well_nseg;
// we need to know for each segment, how many perforation it has and how many segments using it as outlet_segment
// that is why I think we should use a well model to initialize the WellState here
std::vector<std::vector<int>> segment_perforations(nseg);
std::vector<std::vector<int>> segment_perforations(well_nseg);
for (int perf = 0; perf < nperf; ++perf) {
const Completion& completion = completion_set.get(perf);
const int segment_number = completion.getSegmentNumber();
@ -355,8 +356,8 @@ namespace Opm
segment_perforations[segment_index].push_back(perf);
}
std::vector<std::vector<int>> segment_inlets(nseg);
for (int seg = 0; seg < nseg; ++seg) {
std::vector<std::vector<int>> segment_inlets(well_nseg);
for (int seg = 0; seg < well_nseg; ++seg) {
const Segment& segment = segment_set[seg];
const int segment_number = segment.segmentNumber();
const int outlet_segment_number = segment.outletSegment();
@ -394,10 +395,9 @@ namespace Opm
std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(segrates_));
}
// for the segment pressure, the segment pressure is the same with the first perforation
// for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment
// if there is no perforation associated with this segment, it uses the pressure from the outlet segment
// which requres the ordering is successful
// TODO: maybe also check the relation with the cell?
// Not sure what is the best way to handle the initialization, hopefully, the bad initialization can be
// improved during the solveWellEq process
{
@ -405,7 +405,7 @@ namespace Opm
segpress_.push_back(bhp()[w]);
const int top_segment = top_segment_index_[w];
const int start_perf = wells->well_connpos[w];
for (int seg = 1; seg < nseg; ++seg) {
for (int seg = 1; seg < well_nseg; ++seg) {
if ( !segment_perforations[seg].empty() ) {
const int first_perf = segment_perforations[seg][0];
segpress_.push_back(perfPress()[start_perf + first_perf]);
@ -439,6 +439,7 @@ namespace Opm
const int old_top_segment_index = prev_well_state.topSegmentIndex(old_index_well);
const int new_top_segmnet_index = topSegmentIndex(new_index_well);
int number_of_segment = 0;
// if it is the last well in list
if (new_index_well == int(top_segment_index_.size()) - 1) {
number_of_segment = nseg_ - new_top_segmnet_index;
} else {
@ -464,9 +465,9 @@ namespace Opm
// the rate of the segment equals to the sum of the contribution from the perforations and inlet segment rates.
// the first segment is always the top segment, its rates should be equal to the well rates.
assert(segment_inlets.size() == segment_perforations.size());
const int nseg = segment_inlets.size();
const int well_nseg = segment_inlets.size();
if (segment == 0) { // beginning the calculation
segment_rates.resize(np * nseg, 0.0);
segment_rates.resize(np * well_nseg, 0.0);
}
// contributions from the perforations belong to this segment
for (const int& perf : segment_perforations[segment]) {
@ -550,13 +551,13 @@ namespace Opm
std::vector<bool> is_new_well_;
// MS well related
// for StandardWell, the segment number will be one
// for StandardWell, the number of segments will be one
std::vector<double> segrates_;
std::vector<double> segpress_;
// the index of the top segments, which is used to locate the
// multisegment well related information in WellState
std::vector<int> top_segment_index_;
int nseg_; // number of the segments
int nseg_; // total number of the segments
};

11
tests/test_wellmodel.cpp
View File

@ -47,6 +47,7 @@
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/createGlobalCellArray.hpp>
#include <opm/autodiff/GridInit.hpp>
@ -122,9 +123,10 @@ BOOST_AUTO_TEST_CASE(TestStandardWellInput) {
const auto& wells_ecl = setup_test.ecl_state->getSchedule().getWells(setup_test.current_timestep);
BOOST_CHECK_EQUAL( wells_ecl.size(), 2);
const Opm::Well* well = wells_ecl[1];
BOOST_CHECK_THROW( StandardWell( well, -1, wells), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr), std::invalid_argument);
const Opm::BlackoilModelParameters param;
BOOST_CHECK_THROW( StandardWell( well, -1, wells, param), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( nullptr, 4, wells, param), std::invalid_argument);
BOOST_CHECK_THROW( StandardWell( well, 4, nullptr, param), std::invalid_argument);
}
@ -137,6 +139,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
{
const int nw = wells_struct ? (wells_struct->number_of_wells) : 0;
const Opm::BlackoilModelParameters param;
for (int w = 0; w < nw; ++w) {
const std::string well_name(wells_struct->name[w]);
@ -150,7 +153,7 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
// we should always be able to find the well in wells_ecl
BOOST_CHECK(index_well != wells_ecl.size());
wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct) );
wells.emplace_back(new StandardWell(wells_ecl[index_well], current_timestep, wells_struct, param) );
}
}