diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 893bffa12..692905f2a 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -47,7 +47,6 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/VFPProdProperties.cpp opm/autodiff/VFPInjProperties.cpp opm/autodiff/WellMultiSegment.cpp - opm/autodiff/StandardWells.cpp opm/autodiff/BlackoilSolventState.cpp opm/autodiff/ThreadHandle.hpp opm/polymer/PolymerState.cpp @@ -194,7 +193,9 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/VFPInjProperties.hpp opm/autodiff/WellStateMultiSegment.hpp opm/autodiff/WellMultiSegment.hpp + opm/autodiff/WellHelpers.hpp opm/autodiff/StandardWells.hpp + opm/autodiff/StandardWellsSolvent.hpp opm/polymer/CompressibleTpfaPolymer.hpp opm/polymer/GravityColumnSolverPolymer.hpp opm/polymer/GravityColumnSolverPolymer_impl.hpp diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index b8ceae236..41af2dc3c 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -383,39 +383,28 @@ namespace Opm { void computeWellConnectionPressures(const SolutionState& state, const WellState& xw); - void computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf); - void assembleMassBalanceEq(const SolutionState& state); + // TODO: only kept for now due to flow_multisegment + // will be removed soon void extractWellPerfProperties(const SolutionState& state, std::vector& mob_perfcells, std::vector& b_perfcells) const; + // TODO: only kept for now due to flow_multisegment + // will be removed soon + void updateWellState(const V& dwells, + WellState& well_state); + + bool solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, SolutionState& state, WellState& well_state); - void - computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s) const; - - void - updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - const SolutionState& state, - WellState& xw) const; - void addWellFluxEq(const std::vector& cq_s, const SolutionState& state); @@ -430,11 +419,6 @@ namespace Opm { const WellState& xw, const V& aliveWells); - void updateWellControls(WellState& xw) const; - - void updateWellState(const V& dwells, - WellState& well_state); - bool getWellConvergence(const int iteration); bool isVFPActive() const; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 8711544b1..7dc6dc2b8 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -758,81 +759,8 @@ namespace detail { } } - template - void BlackoilModelBase::computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf) - { - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; - - // Compute the average pressure in each well block - const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); - V avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - avg_press[perf] = p_avg; - } - } - - const std::vector& well_cells = stdWells().wellOps().well_cells; - - // Use cell values for the temperature as the wells don't knows its temperature yet. - const ADB perf_temp = subset(state.temperature, well_cells); - - // Compute b, rsmax, rvmax values for perforations. - // Evaluate the properties using average well block pressures - // and cell values for rs, rv, phase condition and temperature. - const ADB avg_press_ad = ADB::constant(avg_press); - std::vector perf_cond(nperf); - const std::vector& pc = phaseCondition(); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = pc[well_cells[perf]]; - } - const PhaseUsage& pu = fluid_.phaseUsage(); - DataBlock b(nperf, pu.num_phases); - if (pu.phase_used[BlackoilPhases::Aqua]) { - const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; - } - assert(active_[Oil]); - const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); - if (pu.phase_used[BlackoilPhases::Liquid]) { - const ADB perf_rs = subset(state.rs, well_cells); - const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; - const V rssat = fluidRsSat(avg_press, perf_so, well_cells); - rsmax_perf.assign(rssat.data(), rssat.data() + nperf); - } else { - rsmax_perf.assign(nperf, 0.0); - } - if (pu.phase_used[BlackoilPhases::Vapour]) { - const ADB perf_rv = subset(state.rv, well_cells); - const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; - const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); - rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); - } else { - rvmax_perf.assign(nperf, 0.0); - } - // b is row major, so can just copy data. - b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); - // Surface density. - // The compute density segment wants the surface densities as - // an np * number of wells cells array - V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); - for (int phase = 1; phase < pu.num_phases; ++phase) { - rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); - } template @@ -849,33 +777,19 @@ namespace detail { std::vector rsmax_perf; std::vector rvmax_perf; std::vector surf_dens_perf; - asImpl().computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - // 2. Compute densities - std::vector cd = - WellDensitySegmented::computeConnectionDensities( - wells(), xw, fluid_.phaseUsage(), - b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - const int nperf = wells().well_connpos[wells().number_of_wells]; - const std::vector& well_cells = stdWells().wellOps().well_cells; + asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // Extract well connection depths. - const V depth = cellCentroidsZToEigen(grid_); - const V pdepth = subset(depth, well_cells); - std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); + const StandardWells::Vector depth = cellCentroidsZToEigen(grid_); + const StandardWells::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); // Gravity - double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); + const double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); - // 3. Compute pressure deltas - std::vector cdp = - WellDensitySegmented::computeConnectionPressureDelta( - wells(), perf_depth, cd, grav); + asImpl().stdWells().computeWellConnectionDensitesPressures(xw, fluid_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, grav); - // 4. Store the results - stdWells().wellPerforationDensities() = Eigen::Map(cd.data(), nperf); - stdWells().wellPerforationPressureDiffs() = Eigen::Map(cdp.data(), nperf); } @@ -903,7 +817,10 @@ namespace detail { // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells - asImpl().updateWellControls(well_state); + // asImpl().updateWellControls(well_state); + // asImpl().stdWells().updateWellControls(well_state); + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + asImpl().stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); // Create the primary variables. SolutionState state = asImpl().variableState(reservoir_state, well_state); @@ -938,15 +855,15 @@ namespace detail { std::vector mob_perfcells; std::vector b_perfcells; - asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells); + asImpl().stdWells().extractWellPerfProperties(state, rq_, fluid_.numPhases(), fluid_, active_, mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); } V aliveWells; std::vector cq_s; - asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); - asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); + asImpl().stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); + asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); asImpl().addWellFluxEq(cq_s, state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); asImpl().addWellControlEq(state, well_state, aliveWells); @@ -1111,190 +1028,6 @@ namespace detail { - - template - void - BlackoilModelBase::computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s) const - { - if( ! localWellsActive() ) return ; - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - V Tw = Eigen::Map(wells().WI, nperf); - const std::vector& well_cells = stdWells().wellOps().well_cells; - - // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = stdWells().wellPerforationPressureDiffs(); - // Extract needed quantities for the perforation cells - const ADB& p_perfcells = subset(state.pressure, well_cells); - const ADB& rv_perfcells = subset(state.rv, well_cells); - const ADB& rs_perfcells = subset(state.rs, well_cells); - - // Perforation pressure - const ADB perfpressure = (stdWells().wellOps().w2p * state.bhp) + cdp; - - // Pressure drawdown (also used to determine direction of flow) - const ADB drawdown = p_perfcells - perfpressure; - - // Compute vectors with zero and ones that - // selects the wanted quantities. - - // selects injection perforations - V selectInjectingPerforations = V::Zero(nperf); - // selects producing perforations - V selectProducingPerforations = V::Zero(nperf); - for (int c = 0; c < nperf; ++c){ - if (drawdown.value()[c] < 0) - selectInjectingPerforations[c] = 1; - else - selectProducingPerforations[c] = 1; - } - - // Handle cross flow - const V numInjectingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); - const V numProducingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); - for (int w = 0; w < nw; ++w) { - if (!wells().allow_cf[w]) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - // Crossflow is not allowed; reverse flow is prevented. - // At least one of the perforation must be open in order to have a meeningful - // equation to solve. For the special case where all perforations have reverse flow, - // and the target rate is non-zero all of the perforations are keept open. - if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { - selectProducingPerforations[perf] = 0.0; - } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ - selectInjectingPerforations[perf] = 0.0; - } - } - } - } - - // HANDLE FLOW INTO WELLBORE - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); - cq_ps[phase] = b_perfcells[phase] * cq_p; - } - if (active_[Oil] && active_[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const ADB cq_psOil = cq_ps[oilpos]; - const ADB cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs_perfcells * cq_psOil; - cq_ps[oilpos] += rv_perfcells * cq_psGas; - } - - // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities - ADB total_mob = mob_perfcells[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob += mob_perfcells[phase]; - } - // injection perforations total volume rates - const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); - - // compute wellbore mixture for injecting perforations - // The wellbore mixture depends on the inflow from the reservoar - // and the well injection rates. - - // compute avg. and total wellbore phase volumetric rates at standard conds - const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); - std::vector wbq(np, ADB::null()); - ADB wbqt = ADB::constant(V::Zero(nw)); - for (int phase = 0; phase < np; ++phase) { - const ADB& q_ps = stdWells().wellOps().p2w * cq_ps[phase]; - const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); - Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); - const int pos = pu.phase_pos[phase]; - wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps; - wbqt += wbq[phase]; - } - // compute wellbore mixture at standard conditions. - Selector notDeadWells_selector(wbqt.value(), Selector::Zero); - std::vector cmix_s(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const int pos = pu.phase_pos[phase]; - cmix_s[phase] = stdWells().wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); - } - - // compute volume ratio between connection at standard conditions - ADB volumeRatio = ADB::constant(V::Zero(nperf)); - const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; - for (int phase = 0; phase < np; ++phase) { - ADB tmp = cmix_s[phase]; - - if (phase == Oil && active_[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - tmp -= rv_perfcells * cmix_s[gaspos] / d; - } - if (phase == Gas && active_[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - tmp -= rs_perfcells * cmix_s[oilpos] / d; - } - volumeRatio += tmp / b_perfcells[phase]; - } - - // injecting connections total volumerates at standard conditions - ADB cqt_is = cqt_i/volumeRatio; - - // connection phase volumerates at standard conditions - cq_s.resize(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; - } - - // check for dead wells (used in the well controll equations) - aliveWells = V::Constant(nw, 1.0); - for (int w = 0; w < nw; ++w) { - if (wbqt.value()[w] == 0) { - aliveWells[w] = 0.0; - } - } - } - - - - - - template - void BlackoilModelBase::updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - const SolutionState& state, - WellState& xw) const - { - if ( !asImpl().localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); - for (int phase = 1; phase < np; ++phase) { - cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); - } - xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); - - // Update the perforation pressures. - const V& cdp = stdWells().wellPerforationPressureDiffs(); - const V perfpressure = (stdWells().wellOps().w2p * state.bhp.value().matrix()).array() + cdp; - xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); - } - - - - - template void BlackoilModelBase::addWellFluxEq(const std::vector& cq_s, const SolutionState& state) @@ -1321,138 +1054,6 @@ namespace detail { - namespace detail - { - inline - double rateToCompare(const std::vector& well_phase_flow_rate, - const int well, - const int num_phases, - const double* distr) - { - double rate = 0.0; - for (int phase = 0; phase < num_phases; ++phase) { - // Important: well_phase_flow_rate is ordered with all phase rates for first - // well first, then all phase rates for second well etc. - rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase]; - } - return rate; - } - - inline - bool constraintBroken(const std::vector& bhp, - const std::vector& thp, - const std::vector& well_phase_flow_rate, - const int well, - const int num_phases, - const WellType& well_type, - const WellControls* wc, - const int ctrl_index) - { - const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index); - const double target = well_controls_iget_target(wc, ctrl_index); - const double* distr = well_controls_iget_distr(wc, ctrl_index); - - bool broken = false; - - switch (well_type) { - case INJECTOR: - { - switch (ctrl_type) { - case BHP: - broken = bhp[well] > target; - break; - - case THP: - broken = thp[well] > target; - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - broken = rateToCompare(well_phase_flow_rate, - well, num_phases, distr) > target; - break; - } - } - break; - - case PRODUCER: - { - switch (ctrl_type) { - case BHP: - broken = bhp[well] < target; - break; - - case THP: - broken = thp[well] < target; - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - // Note that the rates compared below are negative, - // so breaking the constraints means: too high flow rate - // (as for injection). - broken = rateToCompare(well_phase_flow_rate, - well, num_phases, distr) < target; - break; - } - } - break; - - default: - OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells."); - } - - return broken; - } - } // namespace detail - - - namespace detail { - /** - * Simple hydrostatic correction for VFP table - * @param wells - wells struct - * @param w Well number - * @param vfp_table VFP table - * @param well_perforation_densities Densities at well perforations - * @param gravity Gravitational constant (e.g., 9.81...) - */ - inline - double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth, - const ADB::V& well_perforation_densities, const double gravity) { - if ( wells.well_connpos[w] == wells.well_connpos[w+1] ) - { - // This is a well with no perforations. - // If this is the last well we would subscript over the - // bounds below. - // we assume well_perforation_densities to be 0 - return 0; - } - const double well_ref_depth = wells.depth_ref[w]; - const double dh = vfp_ref_depth - well_ref_depth; - const int perf = wells.well_connpos[w]; - const double rho = well_perforation_densities[perf]; - const double dp = rho*gravity*dh; - - return dp; - } - - inline - ADB::V computeHydrostaticCorrection(const Wells& wells, const ADB::V vfp_ref_depth, - const ADB::V& well_perforation_densities, const double gravity) { - const int nw = wells.number_of_wells; - ADB::V retval = ADB::V::Zero(nw); - -#pragma omp parallel for schedule(static) - for (int i=0; i bool BlackoilModelBase::isVFPActive() const { @@ -1485,156 +1086,6 @@ namespace detail { } - template - void BlackoilModelBase::updateWellControls(WellState& xw) const - { - if( ! localWellsActive() ) return ; - - std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; - // Find, for each well, if any constraints are broken. If so, - // switch control to first broken constraint. - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); -#pragma omp parallel for schedule(dynamic) - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; - // The current control in the well state overrides - // the current control set in the Wells struct, which - // is instead treated as a default. - int current = xw.currentControls()[w]; - // Loop over all controls except the current one, and also - // skip any RESERVOIR_RATE controls, since we cannot - // handle those. - const int nwc = well_controls_get_num(wc); - int ctrl_index = 0; - for (; ctrl_index < nwc; ++ctrl_index) { - if (ctrl_index == current) { - // This is the currently used control, so it is - // used as an equation. So this is not used as an - // inequality constraint, and therefore skipped. - continue; - } - if (detail::constraintBroken( - xw.bhp(), xw.thp(), xw.wellRates(), - w, np, wells().type[w], wc, ctrl_index)) { - // ctrl_index will be the index of the broken constraint after the loop. - break; - } - } - if (ctrl_index != nwc) { - // Constraint number ctrl_index was broken, switch to it. - if (terminal_output_) - { - std::cout << "Switching control mode for well " << wells().name[w] - << " from " << modestring[well_controls_iget_type(wc, current)] - << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; - } - xw.currentControls()[w] = ctrl_index; - current = xw.currentControls()[w]; - } - - //Get gravity for THP hydrostatic corrrection - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - - // Updating well state and primary variables. - // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE - const double target = well_controls_iget_target(wc, current); - const double* distr = well_controls_iget_distr(wc, current); - switch (well_controls_iget_type(wc, current)) { - case BHP: - xw.bhp()[w] = target; - break; - - case THP: { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - if (active_[ Water ]) { - aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if (active_[ Oil ]) { - liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if (active_[ Gas ]) { - vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; - } - - const int vfp = well_controls_iget_vfp(wc, current); - const double& thp = well_controls_iget_target(wc, current); - const double& alq = well_controls_iget_alq(wc, current); - - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; - - if (well_type == INJECTOR) { - double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); - - xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - } - else if (well_type == PRODUCER) { - double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); - - xw.bhp()[w] = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - } - else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - break; - } - - case RESERVOIR_RATE: - // No direct change to any observable quantity at - // surface condition. In this case, use existing - // flow rates as initial conditions as reservoir - // rate acts only in aggregate. - break; - - case SURFACE_RATE: - // assign target value as initial guess for injectors and - // single phase producers (orat, grat, wrat) - const WellType& well_type = wells().type[w]; - if (well_type == INJECTOR) { - for (int phase = 0; phase < np; ++phase) { - const double& compi = wells().comp_frac[np * w + phase]; - if (compi > 0.0) { - xw.wellRates()[np*w + phase] = target * compi; - } - } - } else if (well_type == PRODUCER) { - - // only set target as initial rates for single phase - // producers. (orat, grat and wrat, and not lrat) - // lrat will result in numPhasesWithTargetsUnderThisControl == 2 - int numPhasesWithTargetsUnderThisControl = 0; - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - numPhasesWithTargetsUnderThisControl += 1; - } - } - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { - xw.wellRates()[np*w + phase] = target * distr[phase]; - } - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - - - break; - } - - } - } - - - template @@ -1674,8 +1125,8 @@ namespace detail { SolutionState wellSolutionState = state0; asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState); - asImpl().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); - asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); + asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); + asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); asImpl().addWellFluxEq(cq_s, wellSolutionState); asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); converged = getWellConvergence(it); @@ -1700,8 +1151,10 @@ namespace detail { ADB::V total_residual_v = total_residual.value(); const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); assert(dx.size() == total_residual_v.size()); - asImpl().updateWellState(dx.array(), well_state); - asImpl().updateWellControls(well_state); + // asImpl().updateWellState(dx.array(), well_state); + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + asImpl().stdWells().updateWellState(dx.array(), gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state); + asImpl().stdWells(). updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); } } while (it < 15); @@ -1810,7 +1263,7 @@ namespace detail { case THP: { const int perf = wells().well_connpos[w]; - rho_v[w] = stdWells().wellPerforationDensities()[perf]; + rho_v[w] = asImpl().stdWells().wellPerforationDensities()[perf]; const int table_id = well_controls_iget_vfp(wc, current); const double target = well_controls_iget_target(wc, current); @@ -1873,7 +1326,7 @@ namespace detail { //Perform hydrostatic correction to computed targets double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - const ADB::V dp_v = detail::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, stdWells().wellPerforationDensities(), gravity); + const ADB::V dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, asImpl().stdWells().wellPerforationDensities(), gravity); const ADB dp = ADB::constant(dp_v); const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw); @@ -2231,6 +1684,9 @@ namespace detail { } + // TODO: gravity should be stored as a member + // const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + // asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state); asImpl().updateWellState(dwells,well_state); // Update phase conditions used for property calculations. @@ -2241,106 +1697,6 @@ namespace detail { - template - void - BlackoilModelBase::updateWellState(const V& dwells, - WellState& well_state) - { - - if( localWellsActive() ) - { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - // Extract parts of dwells corresponding to each part. - int varstart = 0; - const V dqs = subset(dwells, Span(np*nw, 1, varstart)); - varstart += dqs.size(); - const V dbhp = subset(dwells, Span(nw, 1, varstart)); - varstart += dbhp.size(); - assert(varstart == dwells.size()); - const double dpmaxrel = dpMaxRel(); - - - // Qs update. - // Since we need to update the wellrates, that are ordered by wells, - // from dqs which are ordered by phase, the simplest is to compute - // dwr, which is the data from dqs but ordered by wells. - const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); - const V dwr = Eigen::Map(wwr.data(), nw*np); - const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); - const V wr = wr_old - dwr; - std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); - - // Bhp update. - const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); - const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); - const V bhp = bhp_old - dbhp_limited; - std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); - - //Get gravity for THP hydrostatic correction - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - - // Thp update - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - //Loop over all wells -#pragma omp parallel for schedule(static) - for (int w=0; wgetTable(table_id)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); - - well_state.thp()[w] = vfp_properties_.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp); - } - else if (well_type == PRODUCER) { - double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); - - well_state.thp()[w] = vfp_properties_.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - - //Assume only one THP control specified for each well - break; - } - } - } - } - } - - - - - template std::vector BlackoilModelBase::computeRelPerm(const SolutionState& state) const @@ -3171,6 +2527,24 @@ namespace detail { + + // TODO: only kept for now due to flow_multisegment + // will be removed soon + template + void + BlackoilModelBase::updateWellState(const V& dwells, + WellState& well_state) + { + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), + active_, vfp_properties_, well_state); + + } + + + + + } // namespace Opm #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index b59e709c2..1e2c68c31 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -222,6 +222,7 @@ namespace Opm { using Base::convergenceReduction; using Base::maxResidualAllowed; using Base::variableState; + using Base::variableWellStateIndices; using Base::asImpl; const std::vector& wellsMultiSegment() const { return wells_multisegment_; } @@ -239,6 +240,13 @@ namespace Opm { void computeWellConnectionPressures(const SolutionState& state, const WellState& xw); + /// added to fixing the flow_multisegment running + bool + baseSolveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state); + bool solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index 84d63fa8d..c87c2c7ac 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -964,7 +964,7 @@ namespace Opm { // inequality constraint, and therefore skipped. continue; } - if (detail::constraintBroken( + if (wellhelpers::constraintBroken( xw.bhp(), xw.thp(), xw.wellRates(), w, np, wellsMultiSegment()[w]->wellType(), wc, ctrl_index)) { // ctrl_index will be the index of the broken constraint after the loop. @@ -1033,7 +1033,7 @@ namespace Opm { SolutionState& state, WellState& well_state) { - const bool converged = Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); + const bool converged = baseSolveWellEq(mob_perfcells, b_perfcells, state, well_state); if (converged) { // We must now update the state.segp and state.segqs members, @@ -1547,6 +1547,112 @@ namespace Opm { } + + + /// added to fixing the flow_multisegment running + template + bool + BlackoilMultiSegmentModel::baseSolveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state) { + V aliveWells; + const int np = wells().number_of_phases; + std::vector cq_s(np, ADB::null()); + std::vector indices = variableWellStateIndices(); + SolutionState state0 = state; + WellState well_state0 = well_state; + makeConstantState(state0); + + std::vector mob_perfcells_const(np, ADB::null()); + std::vector b_perfcells_const(np, ADB::null()); + + if ( Base::localWellsActive() ){ + // If there are non well in the sudomain of the process + // thene mob_perfcells_const and b_perfcells_const would be empty + for (int phase = 0; phase < np; ++phase) { + mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); + b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); + } + } + + int it = 0; + bool converged; + do { + // bhp and Q for the wells + std::vector vars0; + vars0.reserve(2); + variableWellStateInitials(well_state, vars0); + std::vector vars = ADB::variables(vars0); + + SolutionState wellSolutionState = state0; + variableStateExtractWellsVars(indices, vars, wellSolutionState); + computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); + updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); + addWellFluxEq(cq_s, wellSolutionState); + addWellControlEq(wellSolutionState, well_state, aliveWells); + converged = Base::getWellConvergence(it); + + if (converged) { + break; + } + + ++it; + if( Base::localWellsActive() ) + { + std::vector eqs; + eqs.reserve(2); + eqs.push_back(residual_.well_flux_eq); + eqs.push_back(residual_.well_eq); + ADB total_residual = vertcatCollapseJacs(eqs); + const std::vector& Jn = total_residual.derivative(); + typedef Eigen::SparseMatrix Sp; + Sp Jn0; + Jn[0].toSparse(Jn0); + const Eigen::SparseLU< Sp > solver(Jn0); + ADB::V total_residual_v = total_residual.value(); + const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); + assert(dx.size() == total_residual_v.size()); + // asImpl().updateWellState(dx.array(), well_state); + updateWellState(dx.array(), well_state); + updateWellControls(well_state); + } + } while (it < 15); + + if (converged) { + if ( terminal_output_ ) { + std::cout << "well converged iter: " << it << std::endl; + } + const int nw = wells().number_of_wells; + { + // We will set the bhp primary variable to the new ones, + // but we do not change the derivatives here. + ADB::V new_bhp = Eigen::Map(well_state.bhp().data(), nw); + // Avoiding the copy below would require a value setter method + // in AutoDiffBlock. + std::vector old_derivs = state.bhp.derivative(); + state.bhp = ADB::function(std::move(new_bhp), std::move(old_derivs)); + } + { + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + // The transpose() below switches the ordering. + const DataBlock wrates = Eigen::Map(well_state.wellRates().data(), nw, np).transpose(); + ADB::V new_qs = Eigen::Map(wrates.data(), nw*np); + std::vector old_derivs = state.qs.derivative(); + state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); + } + computeWellConnectionPressures(state, well_state); + } + + if (!converged) { + well_state = well_state0; + } + + return converged; + } + + } // namespace Opm #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 505f7a3a6..c7881a3a4 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace Opm { @@ -103,6 +104,9 @@ namespace Opm { const bool is_miscible_; std::vector mu_eff_; std::vector b_eff_; + StandardWellsSolvent std_wells_; + const StandardWellsSolvent& stdWells() const { return std_wells_; } + StandardWellsSolvent& stdWells() { return std_wells_; } // Need to declare Base members we want to use here. @@ -130,7 +134,7 @@ namespace Opm { // --------- Protected methods --------- // Need to declare Base members we want to use here. - using Base::stdWells; + // using Base::stdWells; using Base::wells; using Base::variableState; using Base::computeGasPressure; @@ -145,10 +149,10 @@ namespace Opm { using Base::dsMax; using Base::drMaxRel; using Base::maxResidualAllowed; - using Base::updateWellControls; + // using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; - using Base::computePropertiesForWellConnectionPressures; + // using Base::computePropertiesForWellConnectionPressures; std::vector computeRelPerm(const SolutionState& state) const; @@ -202,13 +206,6 @@ namespace Opm { const SolutionState& state, WellState& xw); - void computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf); - void updateEquationsScaling(); void diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index c292c0626..bb86bc578 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -89,7 +89,8 @@ namespace Opm { has_solvent_(has_solvent), solvent_pos_(detail::solventPos(fluid.phaseUsage())), solvent_props_(solvent_props), - is_miscible_(is_miscible) + is_miscible_(is_miscible), + std_wells_(wells_arg, solvent_props, solvent_pos_) { if (has_solvent_) { @@ -381,132 +382,6 @@ namespace Opm { } } - template - void BlackoilSolventModel::computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf) - { - using namespace Opm::AutoDiffGrid; - // 1. Compute properties required by computeConnectionPressureDelta(). - // Note that some of the complexity of this part is due to the function - // taking std::vector arguments, and not Eigen objects. - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - - // Compute the average pressure in each well block - const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); - V avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - avg_press[perf] = p_avg; - } - } - - // Use cell values for the temperature as the wells don't knows its temperature yet. - const ADB perf_temp = subset(state.temperature, well_cells); - - // Compute b, rsmax, rvmax values for perforations. - // Evaluate the properties using average well block pressures - // and cell values for rs, rv, phase condition and temperature. - const ADB avg_press_ad = ADB::constant(avg_press); - std::vector perf_cond(nperf); - const std::vector& pc = phaseCondition(); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = pc[well_cells[perf]]; - } - - const PhaseUsage& pu = fluid_.phaseUsage(); - DataBlock b(nperf, pu.num_phases); - - const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); - if (pu.phase_used[BlackoilPhases::Aqua]) { - b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; - } - - assert(active_[Oil]); - assert(active_[Gas]); - const ADB perf_rv = subset(state.rv, well_cells); - const ADB perf_rs = subset(state.rs, well_cells); - const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); - if (pu.phase_used[BlackoilPhases::Liquid]) { - const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); - //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; - const V rssat = fluidRsSat(avg_press, perf_so, well_cells); - rsmax_perf.assign(rssat.data(), rssat.data() + nperf); - } else { - rsmax_perf.assign(0.0, nperf); - } - V surf_dens_copy = superset(fluid_.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); - for (int phase = 1; phase < pu.num_phases; ++phase) { - if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) { - continue; // the gas surface density is added after the solvent is accounted for. - } - surf_dens_copy += superset(fluid_.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - - if (pu.phase_used[BlackoilPhases::Vapour]) { - // Unclear wether the effective or the pure values should be used for the wells - // the current usage of unmodified properties values gives best match. - //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value(); - V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); - V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells); - if (has_solvent_) { - - const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value(); - //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value(); - - // A weighted sum of the b-factors of gas and solvent are used. - const int nc = Opm::AutoDiffGrid::numCells(grid_); - - const ADB zero = ADB::constant(V::Zero(nc)); - const ADB& ss = state.solvent_saturation; - const ADB& sg = (active_[ Gas ] - ? state.saturation[ pu.phase_pos[ Gas ] ] - : zero); - - Selector zero_selector(ss.value() + sg.value(), Selector::Zero); - V F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value(); - - V injectedSolventFraction = Eigen::Map(&xw.solventFraction()[0], nperf); - - V isProducer = V::Zero(nperf); - V ones = V::Constant(nperf,1.0); - for (int w = 0; w < nw; ++w) { - if(wells().type[w] == PRODUCER) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - isProducer[perf] = 1; - } - } - } - - F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction; - - bg = bg * (ones - F_solvent); - bg = bg + F_solvent * bs; - - const V& rhos = solvent_props_.solventSurfaceDensity(well_cells); - rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos); - } - b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; - surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases); - - const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); - rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); - } else { - rvmax_perf.assign(0.0, nperf); - } - - // b and surf_dens_perf is row major, so can just copy data. - b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); - surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases); - } diff --git a/opm/autodiff/StandardWells.cpp b/opm/autodiff/StandardWells.cpp deleted file mode 100644 index 6b49be6ab..000000000 --- a/opm/autodiff/StandardWells.cpp +++ /dev/null @@ -1,159 +0,0 @@ -/* - Copyright 2016 SINTEF ICT, Applied Mathematics. - Copyright 2016 Statoil ASA. - - 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 . -*/ - - -#include - - - -namespace Opm -{ - - - StandardWells:: - WellOps::WellOps(const Wells* wells) - : w2p(), - p2w(), - well_cells() - { - if( wells ) - { - w2p = Eigen::SparseMatrix(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); - p2w = Eigen::SparseMatrix(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); - - const int nw = wells->number_of_wells; - const int* const wpos = wells->well_connpos; - - typedef Eigen::Triplet Tri; - - std::vector scatter, gather; - scatter.reserve(wpos[nw]); - gather .reserve(wpos[nw]); - - for (int w = 0, i = 0; w < nw; ++w) { - for (; i < wpos[ w + 1 ]; ++i) { - scatter.push_back(Tri(i, w, 1.0)); - gather .push_back(Tri(w, i, 1.0)); - } - } - - w2p.setFromTriplets(scatter.begin(), scatter.end()); - p2w.setFromTriplets(gather .begin(), gather .end()); - - well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]); - } - } - - - - - - StandardWells::StandardWells(const Wells* wells_arg) - : wells_(wells_arg) - , wops_(wells_arg) - , well_perforation_densities_(Vector()) - , well_perforation_pressure_diffs_(Vector()) - { - } - - - - - - const Wells& StandardWells::wells() const - { - assert(wells_ != 0); - return *(wells_); - } - - - - - - bool StandardWells::wellsActive() const - { - return wells_active_; - } - - - - - - void StandardWells::setWellsActive(const bool wells_active) - { - wells_active_ = wells_active; - } - - - - - - bool StandardWells::localWellsActive() const - { - return wells_ ? (wells_->number_of_wells > 0 ) : false; - } - - - - - - const StandardWells::WellOps& - StandardWells::wellOps() const - { - return wops_; - } - - - - - - Vector& StandardWells::wellPerforationDensities() - { - return well_perforation_densities_; - } - - - - - - const Vector& StandardWells::wellPerforationDensities() const - { - return well_perforation_densities_; - } - - - - - - Vector& StandardWells::wellPerforationPressureDiffs() - { - return well_perforation_pressure_diffs_; - } - - - - - - const Vector& StandardWells::wellPerforationPressureDiffs() const - { - return well_perforation_pressure_diffs_; - } - -} diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 1c3d556b7..33a15812f 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -31,13 +31,12 @@ #include #include +#include +#include namespace Opm { - // --------- Types --------- - typedef AutoDiffBlock ADB; - typedef ADB::V Vector; /// Class for handling the standard well model. class StandardWells { @@ -50,6 +49,16 @@ namespace Opm { }; public: + // --------- Types --------- + using ADB = AutoDiffBlock; + using Vector = ADB::V; + + // copied from BlackoilModelBase + // should put to somewhere better + using DataBlock = Eigen::Array; // --------- Public methods --------- explicit StandardWells(const Wells* wells); @@ -64,13 +73,75 @@ namespace Opm { const WellOps& wellOps() const; /// Density of each well perforation - Vector& wellPerforationDensities(); + Vector& wellPerforationDensities(); // mutable version kept for BlackoilMultisegmentModel const Vector& wellPerforationDensities() const; /// Diff to bhp for each well perforation. - Vector& wellPerforationPressureDiffs(); + Vector& wellPerforationPressureDiffs(); // mutable version kept for BlackoilMultisegmentModel const Vector& wellPerforationPressureDiffs() const; + template + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); + + template + void computeWellConnectionDensitesPressures(const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav); + + template + void extractWellPerfProperties(const SolutionState& state, + const std::vector& rq, + const int np, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + std::vector& mob_perfcells, + std::vector& b_perfcells) const; + + template + void computeWellFlux(const SolutionState& state, + const Opm::PhaseUsage& phase_usage, + const std::vector& active, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + std::vector& cq_s) const; + + template + void updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const; + + template + void updateWellState(const Vector& dwells, + const double gravity, + const double dpmaxrel, + const Opm::PhaseUsage& pu, + const std::vector& active, + const VFPProperties& vfp_properties, + WellState& well_state); + + template + void updateWellControls(const Opm::PhaseUsage& pu, + const double gravity, + const VFPProperties& vfp_properties, + const bool terminal_output, + const std::vector& active, + WellState& xw) const; + + protected: bool wells_active_; const Wells* wells_; @@ -81,4 +152,7 @@ namespace Opm { } // namespace Opm + +#include "StandardWells_impl.hpp" + #endif diff --git a/opm/autodiff/StandardWellsSolvent.hpp b/opm/autodiff/StandardWellsSolvent.hpp new file mode 100644 index 000000000..083b30ab7 --- /dev/null +++ b/opm/autodiff/StandardWellsSolvent.hpp @@ -0,0 +1,72 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + 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 . +*/ + + +#ifndef OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED +#define OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED + +#include +#include + +namespace Opm { + + + /// Class for handling the standard well model for solvent model + class StandardWellsSolvent : public StandardWells + { + public: + + using Base = StandardWells; + + // --------- Public methods --------- + explicit StandardWellsSolvent(const Wells* wells, const SolventPropsAdFromDeck& solvent_props, const int solvent_pos); + + template + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); + + // TODO: fluid and active may be can put in the member list + template + void extractWellPerfProperties(const SolutionState& state, + const std::vector& rq, + const int np, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + std::vector& mob_perfcells, + std::vector& b_perfcells) const; + protected: + const SolventPropsAdFromDeck& solvent_props_; + const int solvent_pos_; + + }; + + +} // namespace Opm + +#include "StandardWellsSolvent_impl.hpp" + +#endif diff --git a/opm/autodiff/StandardWellsSolvent_impl.hpp b/opm/autodiff/StandardWellsSolvent_impl.hpp new file mode 100644 index 000000000..ee69a342d --- /dev/null +++ b/opm/autodiff/StandardWellsSolvent_impl.hpp @@ -0,0 +1,231 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + 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 . +*/ + + +#include + + + + +namespace Opm +{ + + + + + + + StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, + const SolventPropsAdFromDeck& solvent_props, + const int solvent_pos) + : Base(wells_arg) + , solvent_props_(solvent_props) + , solvent_pos_(solvent_pos) + { + } + + + + + + template + void + StandardWellsSolvent:: + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + + // Compute the average pressure in each well block + const Vector perf_press = Eigen::Map(xw.perfPress().data(), nperf); + Vector avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + } + + const std::vector& well_cells = wellOps().well_cells; + + // Use cell values for the temperature as the wells don't knows its temperature yet. + const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); + std::vector perf_cond(nperf); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = pc[well_cells[perf]]; + } + + const PhaseUsage& pu = fluid.phaseUsage(); + DataBlock b(nperf, pu.num_phases); + + const Vector bw = fluid.bWat(avg_press_ad, perf_temp, well_cells).value(); + if (pu.phase_used[BlackoilPhases::Aqua]) { + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; + } + + assert(active[Oil]); + assert(active[Gas]); + const ADB perf_rv = subset(state.rv, well_cells); + const ADB perf_rs = subset(state.rs, well_cells); + const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); + if (pu.phase_used[BlackoilPhases::Liquid]) { + const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + // const Vector rssat = fluidRsSat(avg_press, perf_so, well_cells); + const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + } else { + rsmax_perf.assign(0.0, nperf); + } + V surf_dens_copy = superset(fluid.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) { + continue; // the gas surface density is added after the solvent is accounted for. + } + surf_dens_copy += superset(fluid.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + // Unclear wether the effective or the pure values should be used for the wells + // the current usage of unmodified properties values gives best match. + //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value(); + Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + Vector rhog = fluid.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells); + // to handle solvent related + { + + const Vector bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value(); + //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value(); + + // number of cells + const int nc = state.pressure.size(); + + const ADB zero = ADB::constant(Vector::Zero(nc)); + const ADB& ss = state.solvent_saturation; + const ADB& sg = (active[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : zero); + + Selector zero_selector(ss.value() + sg.value(), Selector::Zero); + Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value(); + + Vector injectedSolventFraction = Eigen::Map(&xw.solventFraction()[0], nperf); + + Vector isProducer = Vector::Zero(nperf); + Vector ones = Vector::Constant(nperf,1.0); + for (int w = 0; w < nw; ++w) { + if(wells().type[w] == PRODUCER) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + isProducer[perf] = 1; + } + } + } + + F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction; + + bg = bg * (ones - F_solvent); + bg = bg + F_solvent * bs; + + const Vector& rhos = solvent_props_.solventSurfaceDensity(well_cells); + rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos); + } + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases); + + // const Vector rvsat = fluidRvSat(avg_press, perf_so, well_cells); + const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); + } else { + rvmax_perf.assign(0.0, nperf); + } + + // b and surf_dens_perf is row major, so can just copy data. + b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases); + } + + + + + + + template + void + StandardWellsSolvent:: + extractWellPerfProperties(const SolutionState& state, + const std::vector& rq, + const int np, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + std::vector& mob_perfcells, + std::vector& b_perfcells) const + { + Base::extractWellPerfProperties(state, rq, np, fluid, active, mob_perfcells, b_perfcells); + // handle the solvent related + { + int gas_pos = fluid.phaseUsage().phase_pos[Gas]; + const std::vector& well_cells = wellOps().well_cells; + const int nperf = well_cells.size(); + // Gas and solvent is combinded and solved together + // The input in the well equation is then the + // total gas phase = hydro carbon gas + solvent gas + + // The total mobility is the sum of the solvent and gas mobiliy + mob_perfcells[gas_pos] += subset(rq[solvent_pos_].mob, well_cells); + + // A weighted sum of the b-factors of gas and solvent are used. + const int nc = rq[solvent_pos_].mob.size(); + + const Opm::PhaseUsage& pu = fluid.phaseUsage(); + const ADB zero = ADB::constant(Vector::Zero(nc)); + const ADB& ss = state.solvent_saturation; + const ADB& sg = (active[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : zero); + + Selector zero_selector(ss.value() + sg.value(), Selector::Zero); + ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells); + Vector ones = Vector::Constant(nperf,1.0); + + b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos]; + b_perfcells[gas_pos] += (F_solvent * subset(rq[solvent_pos_].b, well_cells)); + } + } + + +} diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp new file mode 100644 index 000000000..b086984f5 --- /dev/null +++ b/opm/autodiff/StandardWells_impl.hpp @@ -0,0 +1,763 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + 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 . +*/ + + +#include +#include + +#include +#include +#include + + + + +namespace Opm +{ + + + StandardWells:: + WellOps::WellOps(const Wells* wells) + : w2p(), + p2w(), + well_cells() + { + if( wells ) + { + w2p = Eigen::SparseMatrix(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); + p2w = Eigen::SparseMatrix(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); + + const int nw = wells->number_of_wells; + const int* const wpos = wells->well_connpos; + + typedef Eigen::Triplet Tri; + + std::vector scatter, gather; + scatter.reserve(wpos[nw]); + gather .reserve(wpos[nw]); + + for (int w = 0, i = 0; w < nw; ++w) { + for (; i < wpos[ w + 1 ]; ++i) { + scatter.push_back(Tri(i, w, 1.0)); + gather .push_back(Tri(w, i, 1.0)); + } + } + + w2p.setFromTriplets(scatter.begin(), scatter.end()); + p2w.setFromTriplets(gather .begin(), gather .end()); + + well_cells.assign(wells->well_cells, wells->well_cells + wells->well_connpos[wells->number_of_wells]); + } + } + + + + + + StandardWells::StandardWells(const Wells* wells_arg) + : wells_(wells_arg) + , wops_(wells_arg) + , well_perforation_densities_(Vector()) + , well_perforation_pressure_diffs_(Vector()) + { + } + + + + + + const Wells& StandardWells::wells() const + { + assert(wells_ != 0); + return *(wells_); + } + + + + + + bool StandardWells::wellsActive() const + { + return wells_active_; + } + + + + + + void StandardWells::setWellsActive(const bool wells_active) + { + wells_active_ = wells_active; + } + + + + + + bool StandardWells::localWellsActive() const + { + return wells_ ? (wells_->number_of_wells > 0 ) : false; + } + + + + + + const StandardWells::WellOps& + StandardWells::wellOps() const + { + return wops_; + } + + + + + + StandardWells::Vector& StandardWells::wellPerforationDensities() + { + return well_perforation_densities_; + } + + + + + + const StandardWells::Vector& + StandardWells::wellPerforationDensities() const + { + return well_perforation_densities_; + } + + + + + + StandardWells::Vector& + StandardWells::wellPerforationPressureDiffs() + { + return well_perforation_pressure_diffs_; + } + + + + + + const StandardWells::Vector& + StandardWells::wellPerforationPressureDiffs() const + { + return well_perforation_pressure_diffs_; + } + + + + + template + void + StandardWells:: + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + + // Compute the average pressure in each well block + const Vector perf_press = Eigen::Map(xw.perfPress().data(), nperf); + Vector avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + } + + const std::vector& well_cells = wellOps().well_cells; + + // Use cell values for the temperature as the wells don't knows its temperature yet. + const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); + std::vector perf_cond(nperf); + // const std::vector& pc = phaseCondition(); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = pc[well_cells[perf]]; + } + const PhaseUsage& pu = fluid.phaseUsage(); + DataBlock b(nperf, pu.num_phases); + if (pu.phase_used[BlackoilPhases::Aqua]) { + const Vector bw = fluid.bWat(avg_press_ad, perf_temp, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; + } + assert(active[Oil]); + const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); + if (pu.phase_used[BlackoilPhases::Liquid]) { + const ADB perf_rs = subset(state.rs, well_cells); + const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + // const V rssat = fluidRsSat(avg_press, perf_so, well_cells); + const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + } else { + rsmax_perf.assign(nperf, 0.0); + } + if (pu.phase_used[BlackoilPhases::Vapour]) { + const ADB perf_rv = subset(state.rv, well_cells); + const Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + // const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); + const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); + } else { + rvmax_perf.assign(nperf, 0.0); + } + // b is row major, so can just copy data. + b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + + // Surface density. + // The compute density segment wants the surface densities as + // an np * number of wells cells array + Vector rho = superset(fluid.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + rho += superset(fluid.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); + + } + + + + + + template + void + StandardWells:: + computeWellConnectionDensitesPressures(const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav) + { + // Compute densities + std::vector cd = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, fluid.phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const int nperf = wells().well_connpos[wells().number_of_wells]; + + // Compute pressure deltas + std::vector cdp = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), depth_perf, cd, grav); + + // Store the results + well_perforation_densities_ = Eigen::Map(cd.data(), nperf); + well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + } + + + + + + template + void + StandardWells:: + extractWellPerfProperties(const SolutionState& /* state */, + const std::vector& rq, + const int np, + const BlackoilPropsAdInterface& /* fluid */, + const std::vector& /* active */, + std::vector& mob_perfcells, + std::vector& b_perfcells) const + { + // If we have wells, extract the mobilities and b-factors for + // the well-perforated cells. + if ( !localWellsActive() ) { + mob_perfcells.clear(); + b_perfcells.clear(); + return; + } else { + const std::vector& well_cells = wellOps().well_cells; + mob_perfcells.resize(np, ADB::null()); + b_perfcells.resize(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + mob_perfcells[phase] = subset(rq[phase].mob, well_cells); + b_perfcells[phase] = subset(rq[phase].b, well_cells); + } + } + } + + + + + + + template + void + StandardWells:: + computeWellFlux(const SolutionState& state, + const Opm::PhaseUsage& pu, + const std::vector& active, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + std::vector& cq_s) const + { + if( ! localWellsActive() ) return ; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + Vector Tw = Eigen::Map(wells().WI, nperf); + const std::vector& well_cells = wellOps().well_cells; + + // pressure diffs computed already (once per step, not changing per iteration) + const Vector& cdp = wellPerforationPressureDiffs(); + // Extract needed quantities for the perforation cells + const ADB& p_perfcells = subset(state.pressure, well_cells); + const ADB& rv_perfcells = subset(state.rv, well_cells); + const ADB& rs_perfcells = subset(state.rs, well_cells); + + // Perforation pressure + const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp; + + // Pressure drawdown (also used to determine direction of flow) + const ADB drawdown = p_perfcells - perfpressure; + + // Compute vectors with zero and ones that + // selects the wanted quantities. + + // selects injection perforations + Vector selectInjectingPerforations = Vector::Zero(nperf); + // selects producing perforations + Vector selectProducingPerforations = Vector::Zero(nperf); + for (int c = 0; c < nperf; ++c){ + if (drawdown.value()[c] < 0) + selectInjectingPerforations[c] = 1; + else + selectProducingPerforations[c] = 1; + } + + // Handle cross flow + const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); + const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); + for (int w = 0; w < nw; ++w) { + if (!wells().allow_cf[w]) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + // Crossflow is not allowed; reverse flow is prevented. + // At least one of the perforation must be open in order to have a meeningful + // equation to solve. For the special case where all perforations have reverse flow, + // and the target rate is non-zero all of the perforations are keept open. + if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { + selectProducingPerforations[perf] = 0.0; + } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ + selectInjectingPerforations[perf] = 0.0; + } + } + } + } + + // HANDLE FLOW INTO WELLBORE + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); + cq_ps[phase] = b_perfcells[phase] * cq_p; + } + if (active[Oil] && active[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const ADB cq_psOil = cq_ps[oilpos]; + const ADB cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs_perfcells * cq_psOil; + cq_ps[oilpos] += rv_perfcells * cq_psGas; + } + + // HANDLE FLOW OUT FROM WELLBORE + // Using total mobilities + ADB total_mob = mob_perfcells[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob += mob_perfcells[phase]; + } + // injection perforations total volume rates + const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); + + // compute wellbore mixture for injecting perforations + // The wellbore mixture depends on the inflow from the reservoar + // and the well injection rates. + + // compute avg. and total wellbore phase volumetric rates at standard conds + const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); + std::vector wbq(np, ADB::null()); + ADB wbqt = ADB::constant(Vector::Zero(nw)); + for (int phase = 0; phase < np; ++phase) { + const ADB& q_ps = wellOps().p2w * cq_ps[phase]; + const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); + Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); + const int pos = pu.phase_pos[phase]; + wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps; + wbqt += wbq[phase]; + } + // compute wellbore mixture at standard conditions. + Selector notDeadWells_selector(wbqt.value(), Selector::Zero); + std::vector cmix_s(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const int pos = pu.phase_pos[phase]; + cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); + } + + // compute volume ratio between connection at standard conditions + ADB volumeRatio = ADB::constant(Vector::Zero(nperf)); + const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; + for (int phase = 0; phase < np; ++phase) { + ADB tmp = cmix_s[phase]; + + if (phase == Oil && active[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + tmp -= rv_perfcells * cmix_s[gaspos] / d; + } + if (phase == Gas && active[Oil]) { + const int oilpos = pu.phase_pos[Oil]; + tmp -= rs_perfcells * cmix_s[oilpos] / d; + } + volumeRatio += tmp / b_perfcells[phase]; + } + + // injecting connections total volumerates at standard conditions + ADB cqt_is = cqt_i/volumeRatio; + + // connection phase volumerates at standard conditions + cq_s.resize(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; + } + + // check for dead wells (used in the well controll equations) + aliveWells = Vector::Constant(nw, 1.0); + for (int w = 0; w < nw; ++w) { + if (wbqt.value()[w] == 0) { + aliveWells[w] = 0.0; + } + } + } + + + + + + template + void + StandardWells:: + updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const + { + if ( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); + } + xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); + + // Update the perforation pressures. + const Vector& cdp = wellPerforationPressureDiffs(); + const Vector perfpressure = (wellOps().w2p * state.bhp.value().matrix()).array() + cdp; + xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); + } + + + + + + template + void + StandardWells:: + updateWellState(const Vector& dwells, + const double gravity, + const double dpmaxrel, + const Opm::PhaseUsage& pu, + const std::vector& active, + const VFPProperties& vfp_properties, + WellState& well_state) + { + if( localWellsActive() ) + { + // TODO: these parameter should be stored in the StandardWells class + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + // Extract parts of dwells corresponding to each part. + int varstart = 0; + const Vector dqs = subset(dwells, Span(np*nw, 1, varstart)); + varstart += dqs.size(); + const Vector dbhp = subset(dwells, Span(nw, 1, varstart)); + varstart += dbhp.size(); + assert(varstart == dwells.size()); + + + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const Vector dwr = Eigen::Map(wwr.data(), nw*np); + const Vector wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const Vector wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + + // Bhp update. + const Vector bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); + const Vector dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); + const Vector bhp = bhp_old - dbhp_limited; + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + + //Loop over all wells +#pragma omp parallel for schedule(static) + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + const int nwc = well_controls_get_num(wc); + //Loop over all controls until we find a THP control + //that specifies what we need... + //Will only update THP for wells with THP control + for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(wc, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + if (active[ Water ]) { + aqua = wr[w*np + pu.phase_pos[ Water ] ]; + } + if (active[ Oil ]) { + liquid = wr[w*np + pu.phase_pos[ Oil ] ]; + } + if (active[ Gas ]) { + vapour = wr[w*np + pu.phase_pos[ Gas ] ]; + } + + double alq = well_controls_iget_alq(wc, ctrl_index); + int table_id = well_controls_iget_vfp(wc, ctrl_index); + + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getInj()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity); + + well_state.thp()[w] = vfp_properties.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp); + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getProd()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity); + + well_state.thp()[w] = vfp_properties.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + + //Assume only one THP control specified for each well + break; + } + } + } + } + } + + + + + + template + void + StandardWells:: + updateWellControls(const Opm::PhaseUsage& pu, + const double gravity, + const VFPProperties& vfp_properties, + const bool terminal_output, + const std::vector& active, + WellState& xw) const + { + if( !localWellsActive() ) return ; + + std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; +#pragma omp parallel for schedule(dynamic) + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + int current = xw.currentControls()[w]; + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (wellhelpers::constraintBroken( + xw.bhp(), xw.thp(), xw.wellRates(), + w, np, wells().type[w], wc, ctrl_index)) { + // ctrl_index will be the index of the broken constraint after the loop. + break; + } + } + if (ctrl_index != nwc) { + // Constraint number ctrl_index was broken, switch to it. + if (terminal_output) + { + std::cout << "Switching control mode for well " << wells().name[w] + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + } + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + } + + // Updating well state and primary variables. + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[w] = target; + break; + + case THP: { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + if (active[ Water ]) { + aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if (active[ Oil ]) { + liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if (active[ Gas ]) { + vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + const int vfp = well_controls_iget_vfp(wc, current); + const double& thp = well_controls_iget_target(wc, current); + const double& alq = well_controls_iget_alq(wc, current); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity); + + xw.bhp()[w] = vfp_properties.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity); + + xw.bhp()[w] = vfp_properties.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + break; + } + + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + break; + + case SURFACE_RATE: + // assign target value as initial guess for injectors and + // single phase producers (orat, grat, wrat) + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + for (int phase = 0; phase < np; ++phase) { + const double& compi = wells().comp_frac[np * w + phase]; + if (compi > 0.0) { + xw.wellRates()[np*w + phase] = target * compi; + } + } + } else if (well_type == PRODUCER) { + + // only set target as initial rates for single phase + // producers. (orat, grat and wrat, and not lrat) + // lrat will result in numPhasesWithTargetsUnderThisControl == 2 + int numPhasesWithTargetsUnderThisControl = 0; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + + + break; + } + } + + } + +} diff --git a/opm/autodiff/WellHelpers.hpp b/opm/autodiff/WellHelpers.hpp new file mode 100644 index 000000000..b92bb0682 --- /dev/null +++ b/opm/autodiff/WellHelpers.hpp @@ -0,0 +1,170 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + 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 . +*/ + + +#ifndef OPM_WELLHELPERS_HEADER_INCLUDED +#define OPM_WELLHELPERS_HEADER_INCLUDED + +#include +#include +// #include + +#include + +namespace Opm { + + + + + namespace wellhelpers + { + inline + double rateToCompare(const std::vector& well_phase_flow_rate, + const int well, + const int num_phases, + const double* distr) + { + double rate = 0.0; + for (int phase = 0; phase < num_phases; ++phase) { + // Important: well_phase_flow_rate is ordered with all phase rates for first + // well first, then all phase rates for second well etc. + rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase]; + } + return rate; + } + + inline + bool constraintBroken(const std::vector& bhp, + const std::vector& thp, + const std::vector& well_phase_flow_rate, + const int well, + const int num_phases, + const WellType& well_type, + const WellControls* wc, + const int ctrl_index) + { + const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index); + const double target = well_controls_iget_target(wc, ctrl_index); + const double* distr = well_controls_iget_distr(wc, ctrl_index); + + bool broken = false; + + switch (well_type) { + case INJECTOR: + { + switch (ctrl_type) { + case BHP: + broken = bhp[well] > target; + break; + + case THP: + broken = thp[well] > target; + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + broken = rateToCompare(well_phase_flow_rate, + well, num_phases, distr) > target; + break; + } + } + break; + + case PRODUCER: + { + switch (ctrl_type) { + case BHP: + broken = bhp[well] < target; + break; + + case THP: + broken = thp[well] < target; + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + // Note that the rates compared below are negative, + // so breaking the constraints means: too high flow rate + // (as for injection). + broken = rateToCompare(well_phase_flow_rate, + well, num_phases, distr) < target; + break; + } + } + break; + + default: + OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells."); + } + + return broken; + } + + + // --------- Types --------- + using Vector = AutoDiffBlock::V; + + /** + * Simple hydrostatic correction for VFP table + * @param wells - wells struct + * @param w Well number + * @param vfp_table VFP table + * @param well_perforation_densities Densities at well perforations + * @param gravity Gravitational constant (e.g., 9.81...) + */ + inline + double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth, + const Vector& well_perforation_densities, const double gravity) { + if ( wells.well_connpos[w] == wells.well_connpos[w+1] ) + { + // This is a well with no perforations. + // If this is the last well we would subscript over the + // bounds below. + // we assume well_perforation_densities to be 0 + return 0; + } + const double well_ref_depth = wells.depth_ref[w]; + const double dh = vfp_ref_depth - well_ref_depth; + const int perf = wells.well_connpos[w]; + const double rho = well_perforation_densities[perf]; + const double dp = rho*gravity*dh; + + return dp; + } + + inline + Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth, + const Vector& well_perforation_densities, const double gravity) { + const int nw = wells.number_of_wells; + Vector retval = Vector::Zero(nw); + +#pragma omp parallel for schedule(static) + for (int i=0; i water_vel_wells; @@ -569,8 +571,8 @@ namespace Opm { mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; } - Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); - Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); + stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); + stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); Base::addWellFluxEq(cq_s, state); addWellContributionToMassBalanceEq(cq_s, state, well_state); addWellControlEq(state, well_state, aliveWells);