diff --git a/opm/autodiff/BlackoilModelParametersEbos.hpp b/opm/autodiff/BlackoilModelParametersEbos.hpp index ca12f0732..66fd9975c 100644 --- a/opm/autodiff/BlackoilModelParametersEbos.hpp +++ b/opm/autodiff/BlackoilModelParametersEbos.hpp @@ -63,7 +63,7 @@ SET_SCALAR_PROP(FlowModelParameters, ToleranceCnv,1e-2); SET_SCALAR_PROP(FlowModelParameters, ToleranceCnvRelaxed, 1e9); SET_SCALAR_PROP(FlowModelParameters, ToleranceWells, 1e-4); SET_SCALAR_PROP(FlowModelParameters, ToleranceWellControl, 1e-7); -SET_INT_PROP(FlowModelParameters, MaxWelleqIter, 15); +SET_INT_PROP(FlowModelParameters, MaxWelleqIter, 30); SET_BOOL_PROP(FlowModelParameters, UseMultisegmentWell, false); SET_SCALAR_PROP(FlowModelParameters, MaxSinglePrecisionDays, 20.0); SET_INT_PROP(FlowModelParameters, MaxStrictIter, 8); @@ -74,7 +74,7 @@ SET_BOOL_PROP(FlowModelParameters, MatrixAddWellContributions, false); SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01 *1e5); SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 2.0 *1e5); SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true); -SET_INT_PROP(FlowModelParameters, MaxInnerIterMsWells, 10); +SET_INT_PROP(FlowModelParameters, MaxInnerIterMsWells, 100); // if openMP is available, determine the number threads per process automatically. #if _OPENMP diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp index 57411760a..cd506669e 100644 --- a/opm/autodiff/BlackoilWellModel.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -358,7 +358,7 @@ namespace Opm { void computeRepRadiusPerfLength(const Grid& grid, Opm::DeferredLogger& deferred_logger); - void computeAverageFormationFactor(std::vector& B_avg) const; + void computeAverageFormationFactor(std::vector& B_avg) const; void applyVREPGroupControl(); @@ -379,7 +379,7 @@ namespace Opm { /// at the beginning of the time step and no derivatives are included in these quantities void calculateExplicitQuantities(Opm::DeferredLogger& deferred_logger) const; - SimulatorReport solveWellEq(const double dt, Opm::DeferredLogger& deferred_logger); + SimulatorReport solveWellEq(const std::vector& B_avg, const double dt, Opm::DeferredLogger& deferred_logger); void initPrimaryVariablesEvaluation() const; @@ -392,7 +392,7 @@ namespace Opm { void resetWellControlFromState() const; - void assembleWellEq(const double dt, Opm::DeferredLogger& deferred_logger); + void assembleWellEq(const std::vector& B_avg, const double dt, Opm::DeferredLogger& deferred_logger); // some preparation work, mostly related to group control and RESV, // at the beginning of each time step (Not report step) diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 09a637b3e..86a6d0602 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -687,15 +687,19 @@ namespace Opm { // Set the well primary variables based on the value of well solutions initPrimaryVariablesEvaluation(); + std::vector< Scalar > B_avg(numComponents(), Scalar() ); + computeAverageFormationFactor(B_avg); + if (param_.solve_welleq_initially_ && iterationIdx == 0) { // solve the well equations as a pre-processing step - last_report_ = solveWellEq(dt, local_deferredLogger); + last_report_ = solveWellEq(B_avg, dt, local_deferredLogger); + if (initial_step_) { // update the explicit quantities to get the initial fluid distribution in the well correct. calculateExplicitQuantities(local_deferredLogger); prepareTimeStep(local_deferredLogger); - last_report_ = solveWellEq(dt, local_deferredLogger); + last_report_ = solveWellEq(B_avg, dt, local_deferredLogger); initial_step_ = false; } // TODO: should we update the explicit related here again, or even prepareTimeStep(). @@ -703,7 +707,8 @@ namespace Opm { // reservoir state, will tihs be a better place to inialize the explict information? } - assembleWellEq(dt, local_deferredLogger); + assembleWellEq(B_avg, dt, local_deferredLogger); + } catch (std::exception& e) { exception_thrown = 1; } @@ -715,10 +720,10 @@ namespace Opm { template void BlackoilWellModel:: - assembleWellEq(const double dt, Opm::DeferredLogger& deferred_logger) + assembleWellEq(const std::vector& B_avg, const double dt, Opm::DeferredLogger& deferred_logger) { for (auto& well : well_container_) { - well->assembleWellEq(ebosSimulator_, dt, well_state_, deferred_logger); + well->assembleWellEq(ebosSimulator_, B_avg, dt, well_state_, deferred_logger); } } @@ -873,14 +878,10 @@ namespace Opm { template SimulatorReport BlackoilWellModel:: - solveWellEq(const double dt, Opm::DeferredLogger& deferred_logger) + solveWellEq(const std::vector& B_avg, const double dt, Opm::DeferredLogger& deferred_logger) { WellState well_state0 = well_state_; - const int numComp = numComponents(); - std::vector< Scalar > B_avg( numComp, Scalar() ); - computeAverageFormationFactor(B_avg); - const int max_iter = param_.max_welleq_iter_; int it = 0; @@ -888,7 +889,7 @@ namespace Opm { int exception_thrown = 0; do { try { - assembleWellEq(dt, deferred_logger); + assembleWellEq(B_avg, dt, deferred_logger); } catch (std::exception& e) { exception_thrown = 1; } @@ -1443,7 +1444,7 @@ namespace Opm { template void BlackoilWellModel:: - computeAverageFormationFactor(std::vector& B_avg) const + computeAverageFormationFactor(std::vector& B_avg) const { const auto& grid = ebosSimulator_.vanguard().grid(); const auto& gridView = grid.leafGridView(); @@ -1736,7 +1737,7 @@ namespace Opm { const auto& segments = well.segments; // \Note: eventually we need to hanlde the situations that some segments are shut - assert(segment_set.size() == segments.size()); + assert(int(segment_set.size()) == segments.size()); for (const auto& segment : segments) { const int segment_index = segment_set.segmentNumberToIndex(segment.first); diff --git a/opm/autodiff/MultisegmentWell.hpp b/opm/autodiff/MultisegmentWell.hpp index f7186332e..b02863e60 100644 --- a/opm/autodiff/MultisegmentWell.hpp +++ b/opm/autodiff/MultisegmentWell.hpp @@ -113,6 +113,7 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const override; virtual void assembleWellEq(const Simulator& ebosSimulator, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger) override; @@ -243,8 +244,8 @@ namespace Opm // the segment the perforation belongs to std::vector perforation_segment_depth_diffs_; - // the intial component compistion of segments - std::vector > segment_comp_initial_; + // the intial amount of fluids in each segment under surface condition + std::vector > segment_fluid_initial_; // the densities of segment fluids // we should not have this member variable @@ -270,9 +271,10 @@ namespace Opm // updating the well_state based on well solution dwells void updateWellState(const BVectorWell& dwells, - const bool inner_iteration, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const; + Opm::DeferredLogger& deferred_logger, + const double relaxation_factor=1.0) const; + // initialize the segment rates with well rates // when there is no more accurate way to initialize the segment rates, we initialize @@ -280,7 +282,7 @@ namespace Opm void initSegmentRatesWithWellRates(WellState& well_state) const; // computing the accumulation term for later use in well mass equations - void computeInitialComposition(); + void computeInitialSegmentFluids(const Simulator& ebos_simulator); // compute the pressure difference between the perforation and cell center void computePerfCellPressDiffs(const Simulator& ebosSimulator); @@ -316,6 +318,8 @@ namespace Opm EvalWell getSegmentRate(const int seg, const int comp_idx) const; + EvalWell getSegmentRateUpwinding(const int seg, const int comp_idx, const bool upwinding, int& seg_upwind) const; + EvalWell getSegmentGTotal(const int seg) const; // get the mobility for specific perforation @@ -352,6 +356,7 @@ namespace Opm // TODO: try to make ebosSimulator const, as it should be void iterateWellEquations(const Simulator& ebosSimulator, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger); @@ -366,6 +371,16 @@ namespace Opm WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override; virtual void updateWaterThroughput(const double dt, WellState& well_state) const override; + + EvalWell getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const; + + std::vector getWellResiduals(const std::vector& B_avg) const; + + void detectOscillations(const std::vector& measure_history, + const int it, bool& oscillate, bool& stagnate) const; + + double getResidualMeasureValue(const std::vector& residuals) const; + }; } diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index 08501d77f..51cdefcc1 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -39,7 +39,7 @@ namespace Opm , cell_perforation_depth_diffs_(number_of_perforations_, 0.0) , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0) , perforation_segment_depth_diffs_(number_of_perforations_, 0.0) - , segment_comp_initial_(numberOfSegments(), std::vector(num_components_, 0.0)) + , segment_fluid_initial_(numberOfSegments(), std::vector(num_components_, 0.0)) , segment_densities_(numberOfSegments(), 0.0) , segment_viscosities_(numberOfSegments(), 0.0) , segment_mass_rates_(numberOfSegments(), 0.0) @@ -54,18 +54,31 @@ namespace Opm OPM_THROW(std::runtime_error, "polymer is not supported by multisegment well yet"); } - if (Base::has_energy) { + if (Base::has_energy) { OPM_THROW(std::runtime_error, "energy is not supported by multisegment well yet"); - } + } // since we decide to use the WellSegments from the well parser. we can reuse a lot from it. // for other facilities needed but not available from parser, we need to process them here - // initialize the segment_perforations_ + // initialize the segment_perforations_ and update perforation_segment_depth_diffs_ const WellConnections& completion_set = well_ecl_->getConnections(current_step_); - for (int perf = 0; perf < number_of_perforations_; ++perf) { + // index of the perforation within wells struct + // there might be some perforations not active, which causes the number of the perforations in + // well_ecl_ and wells struct different + // the current implementation is a temporary solution for now, it should be corrected from the parser + // side + int i_perf_wells = 0; + perf_depth_.resize(number_of_perforations_, 0.); + for (size_t perf = 0; perf < completion_set.size(); ++perf) { const Connection& connection = completion_set.get(perf); - const int segment_index = segmentNumberToIndex(connection.segment()); - segment_perforations_[segment_index].push_back(perf); + if (connection.state() == WellCompletion::OPEN) { + const int segment_index = segmentNumberToIndex(connection.segment()); + segment_perforations_[segment_index].push_back(i_perf_wells); + perf_depth_[i_perf_wells] = connection.depth(); + const double segment_depth = segmentSet()[segment_index].depth(); + perforation_segment_depth_diffs_[i_perf_wells] = perf_depth_[i_perf_wells] - segment_depth; + i_perf_wells++; + } } // initialize the segment_inlets_ @@ -80,16 +93,6 @@ namespace Opm } } - // callcuate the depth difference between perforations and their segments - perf_depth_.resize(number_of_perforations_, 0.); - for (int seg = 0; seg < numberOfSegments(); ++seg) { - const double segment_depth = segmentSet()[seg].depth(); - for (const int perf : segment_perforations_[seg]) { - perf_depth_[perf] = completion_set.get(perf).depth(); - perforation_segment_depth_diffs_[perf] = perf_depth_[perf] - segment_depth; - } - } - // calculating the depth difference between the segment and its oulet_segments // for the top segment, we will make its zero unless we find other purpose to use this value for (int seg = 1; seg < numberOfSegments(); ++seg) { @@ -234,6 +237,7 @@ namespace Opm void MultisegmentWell:: assembleWellEq(const Simulator& ebosSimulator, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger) @@ -241,7 +245,8 @@ namespace Opm const bool use_inner_iterations = param_.use_inner_iterations_ms_wells_; if (use_inner_iterations) { - iterateWellEquations(ebosSimulator, dt, well_state, deferred_logger); + + iterateWellEquations(ebosSimulator, B_avg, dt, well_state, deferred_logger); } assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, deferred_logger); @@ -462,6 +467,7 @@ namespace Opm report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, name()}); } else if (control_residual > param_.max_residual_allowed_) { report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, name()}); + // TODO: we should distinguish the flux residual or pressure residual here } else if (control_residual > param_.tolerance_wells_) { report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, name()}); } @@ -550,7 +556,7 @@ namespace Opm { BVectorWell xw(1); recoverSolutionWell(x, xw); - updateWellState(xw, false, well_state, deferred_logger); + updateWellState(xw, well_state, deferred_logger); } @@ -584,7 +590,7 @@ namespace Opm template void MultisegmentWell:: - updatePrimaryVariables(const WellState& well_state, Opm::DeferredLogger& deferred_logger) const + updatePrimaryVariables(const WellState& well_state, Opm::DeferredLogger& /* deferred_logger */) const { // TODO: to test using rate conversion coefficients to see if it will be better than // this default one @@ -677,7 +683,7 @@ namespace Opm // which is why we do not put the assembleWellEq here. const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); - updateWellState(dx_well, false, well_state, deferred_logger); + updateWellState(dx_well, well_state, deferred_logger); } @@ -742,13 +748,13 @@ namespace Opm template void MultisegmentWell:: - computeInitialComposition() + computeInitialSegmentFluids(const Simulator& ebos_simulator) { for (int seg = 0; seg < numberOfSegments(); ++seg) { - // TODO: probably it should be numWellEq -1 more accurately, - // while by meaning it should be num_comp + // TODO: trying to reduce the times for the surfaceVolumeFraction calculation + const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value(); for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - segment_comp_initial_[seg][comp_idx] = surfaceVolumeFraction(seg, comp_idx).value(); + segment_fluid_initial_[seg][comp_idx] = surface_volume * surfaceVolumeFraction(seg, comp_idx).value(); } } } @@ -761,14 +767,10 @@ namespace Opm void MultisegmentWell:: updateWellState(const BVectorWell& dwells, - const bool inner_iteration, WellState& well_state, - Opm::DeferredLogger& deferred_logger) const + Opm::DeferredLogger& deferred_logger, + const double relaxation_factor) const { - 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 std::vector > old_primary_variables = primary_variables_; @@ -776,13 +778,15 @@ namespace Opm for (int seg = 0; seg < numberOfSegments(); ++seg) { if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { const int sign = dwells[seg][WFrac] > 0. ? 1 : -1; - const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit); + // const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), relaxation_factor * dFLimit); + const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit); primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited; } if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { const int sign = dwells[seg][GFrac] > 0. ? 1 : -1; - const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit); + // const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), relaxation_factor * dFLimit); + const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit); primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited; } @@ -793,6 +797,7 @@ namespace Opm { const int sign = dwells[seg][SPres] > 0.? 1 : -1; const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change); + // const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); primary_variables_[seg][SPres] = old_primary_variables[seg][SPres] - dx_limited; } @@ -814,11 +819,13 @@ namespace Opm void MultisegmentWell:: calculateExplicitQuantities(const Simulator& ebosSimulator, - const WellState& /* well_state */, + const WellState& well_state, Opm::DeferredLogger& deferred_logger) { + updatePrimaryVariables(well_state, deferred_logger); + initPrimaryVariablesEvaluation(); computePerfCellPressDiffs(ebosSimulator); - computeInitialComposition(); + computeInitialSegmentFluids(ebosSimulator); } @@ -1012,7 +1019,7 @@ namespace Opm // pressure difference between the perforation and the grid cell const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf]; - perf_press = pressure_cell + cell_perf_press_diff; + perf_press = pressure_cell - cell_perf_press_diff; // Pressure drawdown (also used to determine direction of flow) // TODO: not 100% sure about the sign of the seg_perf_press_diff const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff); @@ -1276,6 +1283,10 @@ namespace Opm segment_densities_[seg] = density / volrat; // calculate the mass rates + // TODO: for now, we are not considering the upwinding for this amount + // since how to address the fact that the derivatives is not trivial for now + // and segment_mass_rates_ goes a long way with the frictional pressure loss + // and accelerational pressure loss, which needs some work to handle segment_mass_rates_[seg] = 0.; for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { const EvalWell rate = getSegmentRate(seg, comp_idx); @@ -1313,6 +1324,35 @@ namespace Opm + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + getSegmentRateUpwinding(const int seg, + const int comp_idx, + const bool upwinding, + int& seg_upwind) const + { + // not considering upwinding for the injectors for now + if ((!upwinding) || (well_type_ == INJECTOR) || (primary_variables_evaluation_[seg][GTotal] <= 0.) || (seg == 0)) { + seg_upwind = seg; // using the composition from the seg + return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg, comp_idx); + } + + // assert( seg != 0); // if top segment flowing towards the wrong direction, we are not handling it + + // basically here, it a producer and flow is in the injecting direction + // we will use the compsotion from the outlet segment + const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment()); + seg_upwind = outlet_segment_index; + + // TODO: we can refactor above code to use the following return + return primary_variables_evaluation_[seg][GTotal] * volumeFractionScaled(seg_upwind, comp_idx); + } + + + + + template typename MultisegmentWell::EvalWell MultisegmentWell:: @@ -1758,40 +1798,83 @@ namespace Opm void MultisegmentWell:: iterateWellEquations(const Simulator& ebosSimulator, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger) { - // basically, it only iterate through the equations. - // we update the primary variables - // 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 WellState well_state0 = well_state; + const std::vector residuals0 = getWellResiduals(B_avg); + std::vector > residual_history; + std::vector measure_history; int it = 0; + // relaxation factor + double relaxation_factor = 1.; + const double min_relaxation_factor = 0.2; for (; it < max_iter_number; ++it) { assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, deferred_logger); const BVectorWell dx_well = mswellhelpers::invDXDirect(duneD_, resWell_); - // TODO: use these small values for now, not intend to reach the convergence - // in this stage, but, should we? - // We should try to avoid hard-code values in the code. - // If we want to use the real one, we need to find a way to get them. - // const std::vector B {0.8, 0.8, 0.008}; - const std::vector B {0.5, 0.5, 0.005}; - const auto report = getWellConvergence(B, deferred_logger); + const auto report = getWellConvergence(B_avg, deferred_logger); if (report.converged()) { break; } - updateWellState(dx_well, true, well_state, deferred_logger); + residual_history.push_back(getWellResiduals(B_avg)); + measure_history.push_back(getResidualMeasureValue(residual_history[it])); + + bool is_oscillate = false; + bool is_stagnate = false; + + detectOscillations(measure_history, it, is_oscillate, is_stagnate); + // TODO: maybe we should have more sophiscated strategy to recover the relaxation factor, + // for example, to recover it to be bigger + + if (is_oscillate || is_stagnate) { + // a factor value to reduce the relaxation_factor + const double reduction_mutliplier = 0.9; + relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor); + + // debug output + std::ostringstream sstr; + if (is_stagnate) { + sstr << " well " << name() << " observes stagnation within " << it << "th inner iterations\n"; + } + if (is_oscillate) { + sstr << " well " << name() << " osbserves oscillation within " << it <<"th inner iterations\n"; + } + sstr << " relaxation_factor is " << relaxation_factor << " now\n"; + deferred_logger.debug(sstr.str()); + } + + updateWellState(dx_well, well_state, deferred_logger, relaxation_factor); initPrimaryVariablesEvaluation(); } - // TODO: maybe we should not use these values if they are not converged. + + // TODO: we should decide whether to keep the updated well_state, or recover to use the old well_state + if (it < max_iter_number) { + std::ostringstream sstr; + sstr << " well " << name() << " manage to get converged within " << it << " inner iterations"; + deferred_logger.debug(sstr.str()); + } else { + std::ostringstream sstr; + sstr << " well " << name() << " did not get converged within " << it << " inner iterations \n"; + sstr << " outputting the residual history for well " << name() << " during inner iterations \n"; + for (int i = 0; i < it; ++i) { + const auto& residual = residual_history[i]; + sstr << " residual at " << i << "th iteration "; + for (const auto& res : residual) { + sstr << " " << res; + } + sstr << " " << measure_history[i] << " \n"; + } + deferred_logger.debug(sstr.str()); + } } @@ -1826,14 +1909,14 @@ namespace Opm const int nseg = numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { - // calculating the accumulation term // TODO: without considering the efficiencty factor for now - // volume of the segment + // calculating the accumulation term + // TODO: without considering the efficiencty factor for now { - const double volume = segmentSet()[seg].volume(); + const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg); // for each component for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx]) - + getSegmentRate(seg, comp_idx); + const EvalWell accumulation_term = (segment_surface_volume * surfaceVolumeFraction(seg, comp_idx) + - segment_fluid_initial_[seg][comp_idx]) / dt; resWell_[seg][comp_idx] += accumulation_term.value(); for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { @@ -1841,16 +1924,35 @@ namespace Opm } } } + // considering the contributions due to flowing out from the segment + { + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { + int seg_upwind = -1; + const EvalWell segment_rate = getSegmentRateUpwinding(seg, comp_idx, true, seg_upwind); + + assert(seg_upwind >= 0); + + resWell_[seg][comp_idx] -= segment_rate.value(); + duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + numEq); + duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq); + duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq); + // pressure derivative should be zero + } + } // considering the contributions from the inlet segments { for (const int inlet : segment_inlets_[seg]) { for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { - const EvalWell inlet_rate = getSegmentRate(inlet, comp_idx); - resWell_[seg][comp_idx] -= inlet_rate.value(); - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - duneD_[seg][inlet][comp_idx][pv_idx] -= inlet_rate.derivative(pv_idx + numEq); - } + int seg_upwind = -1; + const EvalWell inlet_rate = getSegmentRateUpwinding(inlet, comp_idx, true, seg_upwind); + assert(seg_upwind >= 0); + + resWell_[seg][comp_idx] += inlet_rate.value(); + duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + numEq); + duneD_[seg][seg_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq); + duneD_[seg][seg_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq); + // pressure derivative should be zero } } } @@ -1880,7 +1982,7 @@ namespace Opm connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective); // subtract sum of phase fluxes in the well equations. - resWell_[seg][comp_idx] -= cq_s_effective.value(); + resWell_[seg][comp_idx] += cq_s_effective.value(); // assemble the jacobians for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { @@ -1889,17 +1991,14 @@ namespace Opm duneC_[seg][cell_idx][pv_idx][comp_idx] -= cq_s_effective.derivative(pv_idx + numEq); // intput in transformed matrix // the index name for the D should be eq_idx / pv_idx - duneD_[seg][seg][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx + numEq); + duneD_[seg][seg][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx + numEq); } for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { // also need to consider the efficiency factor when manipulating the jacobians. - duneB_[seg][cell_idx][comp_idx][pv_idx] -= cq_s_effective.derivative(pv_idx); + duneB_[seg][cell_idx][comp_idx][pv_idx] += cq_s_effective.derivative(pv_idx); } } - // TODO: we should save the perforation pressure and preforation rates? - // we do not use it in the simulation for now, while we might need them if - // we handle the pressure in SEG mode. } // the fourth dequation, the pressure drop equation @@ -1938,4 +2037,239 @@ namespace Opm { } + + + + + template + typename MultisegmentWell::EvalWell + MultisegmentWell:: + getSegmentSurfaceVolume(const Simulator& ebos_simulator, const int seg_idx) const + { + EvalWell temperature; + int pvt_region_index; + { + // using the pvt region of first perforated cell + // TODO: it should be a member of the WellInterface, initialized properly + const int cell_idx = well_cells_[0]; + const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); + const auto& fs = intQuants.fluidState(); + temperature.setValue(fs.temperature(FluidSystem::oilPhaseIdx).value()); + pvt_region_index = fs.pvtRegionIndex(); + } + + const EvalWell seg_pressure = getSegmentPressure(seg_idx); + + std::vector mix_s(num_components_, 0.0); + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { + mix_s[comp_idx] = surfaceVolumeFraction(seg_idx, comp_idx); + } + + std::vector b(num_components_, 0.); + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + b[waterCompIdx] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + } + + EvalWell rv(0.0); + // gas phase + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure); + if (mix_s[oilCompIdx] > 0.0) { + if (mix_s[gasCompIdx] > 0.0) { + rv = mix_s[oilCompIdx] / mix_s[gasCompIdx]; + } + + if (rv > rvmax) { + rv = rvmax; + } + b[gasCompIdx] = + FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv); + } else { // no oil exists + b[gasCompIdx] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + } + } else { // no Liquid phase + // it is the same with zero mix_s[Oil] + b[gasCompIdx] = + FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + } + } + + EvalWell rs(0.0); + // oil phase + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure); + if (mix_s[gasCompIdx] > 0.0) { + if (mix_s[oilCompIdx] > 0.0) { + rs = mix_s[gasCompIdx] / mix_s[oilCompIdx]; + } + // std::cout << " rs " << rs.value() << " rsmax " << rsmax.value() << std::endl; + + if (rs > rsmax) { + rs = rsmax; + } + b[oilCompIdx] = + FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs); + } else { // no oil exists + b[oilCompIdx] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + } + } else { // no gas phase + // it is the same with zero mix_s[Gas] + b[oilCompIdx] = + FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure); + } + } + + std::vector mix(mix_s); + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + + const EvalWell d = 1.0 - rs * rv; + if (d <= 0.0 || d > 1.0) { + OPM_THROW(Opm::NumericalIssue, "Problematic d value " << d << " obtained for well " << name() + << " during convertion to surface volume with rs " << rs + << ", rv " << rv << " and pressure " << seg_pressure + << " obtaining d " << d); + } + + if (rs > 0.0) { // rs > 0.0? + mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; + } + if (rv > 0.0) { // rv > 0.0? + mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; + } + } + + EvalWell vol_ratio(0.0); + for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { + vol_ratio += mix[comp_idx] / b[comp_idx]; + } + + // segment volume + const double volume = segmentSet()[seg_idx].volume(); + + return volume / vol_ratio; + } + + + + + + template + std::vector::Scalar> + MultisegmentWell:: + getWellResiduals(const std::vector& B_avg) const + { + assert(int(B_avg.size() ) == num_components_); + std::vector residuals(numWellEq + 1, 0.0); + + // TODO: maybe we should distinguish the bhp control or rate control equations here + for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { + double residual = 0.; + if (eq_idx < num_components_) { + residual = std::abs(resWell_[seg][eq_idx]) * B_avg[eq_idx]; + } else { + if (seg > 0) { + residual = std::abs(resWell_[seg][eq_idx]); + } + } + if (std::isnan(residual) || std::isinf(residual)) { + OPM_THROW(Opm::NumericalIssue, "nan or inf value for residal get for well " << name() + << " segment " << seg << " eq_idx " << eq_idx); + } + + if (residual > residuals[eq_idx]) { + residuals[eq_idx] = residual; + } + } + } + + // handling the control equation residual + { + const double control_residual = std::abs(resWell_[0][numWellEq - 1]); + if (std::isnan(control_residual) || std::isinf(control_residual)) { + OPM_THROW(Opm::NumericalIssue, "nan or inf value for control residal get for well " << name()); + } + residuals[numWellEq] = control_residual; + } + + return residuals; + } + + + + + + /// Detect oscillation or stagnation based on the residual measure history + template + void + MultisegmentWell:: + detectOscillations(const std::vector& measure_history, + const int it, bool& oscillate, bool& stagnate) const + { + if ( it < 2 ) { + oscillate = false; + stagnate = false; + return; + } + + stagnate = true; + const double F0 = measure_history[it]; + const double F1 = measure_history[it - 1]; + const double F2 = measure_history[it - 2]; + const double d1 = std::abs((F0 - F2) / F0); + const double d2 = std::abs((F0 - F1) / F0); + + const double oscillaton_rel_tol = 0.2; + oscillate = (d1 < oscillaton_rel_tol) && (oscillaton_rel_tol < d2); + + const double stagnation_rel_tol = 1.e-2; + stagnate = std::abs((F1 - F2) / F2) <= stagnation_rel_tol; + } + + + + + + template + double + MultisegmentWell:: + getResidualMeasureValue(const std::vector& residuals) const + { + assert(int(residuals.size()) == numWellEq + 1); + + const double rate_tolerance = param_.tolerance_wells_; + int count = 0; + double sum = 0; + for (int eq_idx = 0; eq_idx < numWellEq - 1; ++eq_idx) { + if (residuals[eq_idx] > rate_tolerance) { + sum += residuals[eq_idx] / rate_tolerance; + ++count; + } + } + + const double pressure_tolerance = param_.tolerance_pressure_ms_wells_; + if (residuals[SPres] > pressure_tolerance) { + sum += residuals[SPres] / pressure_tolerance; + ++count; + } + + // if (count == 0), it should be converged. + assert(count != 0); + + // return sum / double(count); + return sum; + } + } diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 4790261ff..44cb377ca 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -142,6 +142,7 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const override; virtual void assembleWellEq(const Simulator& ebosSimulator, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger) override; diff --git a/opm/autodiff/StandardWellV.hpp b/opm/autodiff/StandardWellV.hpp index 09d918162..0ed849f1a 100644 --- a/opm/autodiff/StandardWellV.hpp +++ b/opm/autodiff/StandardWellV.hpp @@ -143,6 +143,7 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const override; virtual void assembleWellEq(const Simulator& ebosSimulator, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger) override; diff --git a/opm/autodiff/StandardWellV_impl.hpp b/opm/autodiff/StandardWellV_impl.hpp index 427dd379a..b322ae1c7 100644 --- a/opm/autodiff/StandardWellV_impl.hpp +++ b/opm/autodiff/StandardWellV_impl.hpp @@ -480,6 +480,7 @@ namespace Opm void StandardWellV:: assembleWellEq(const Simulator& ebosSimulator, + const std::vector& /* B_avg */, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index be8e4c645..bc391085e 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -445,6 +445,7 @@ namespace Opm void StandardWell:: assembleWellEq(const Simulator& ebosSimulator, + const std::vector& /* B_avg */, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger) diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index a570981aa..186fba2be 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -152,6 +152,7 @@ namespace Opm virtual void solveEqAndUpdateWellState(WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0; virtual void assembleWellEq(const Simulator& ebosSimulator, + const std::vector& B_avg, const double dt, WellState& well_state, Opm::DeferredLogger& deferred_logger diff --git a/opm/autodiff/WellInterface_impl.hpp b/opm/autodiff/WellInterface_impl.hpp index 8b59140b6..d8a53c5c4 100644 --- a/opm/autodiff/WellInterface_impl.hpp +++ b/opm/autodiff/WellInterface_impl.hpp @@ -1182,7 +1182,7 @@ namespace Opm bool converged; WellState well_state0 = well_state; do { - assembleWellEq(ebosSimulator, dt, well_state, deferred_logger); + assembleWellEq(ebosSimulator, B_avg, dt, well_state, deferred_logger); auto report = getWellConvergence(B_avg, deferred_logger); diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index e20325659..d193962c6 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -462,12 +462,11 @@ namespace Opm // we need to create a trival segment related values to avoid there will be some // multi-segment wells added later. nseg_ = nw; - seg_number_.clear(); - top_segment_index_.reserve(nw); - seg_number_.reserve(nw); + top_segment_index_.resize(nw); + seg_number_.resize(nw); for (int w = 0; w < nw; ++w) { - top_segment_index_.push_back(w); - seg_number_.push_back(1); // Top segment is segment #1 + top_segment_index_[w] = w; + seg_number_[w] = 1; // Top segment is segment #1 } segpress_ = bhp(); segrates_ = wellRates(); @@ -657,7 +656,8 @@ namespace Opm const WellConnections& completion_set = well_ecl->getConnections(time_step); // number of segment for this single well const int well_nseg = segment_set.size(); - const int nperf = completion_set.size(); + // const int nperf = completion_set.size(); + int n_activeperf = 0; nseg_ += well_nseg; for (auto segID = 0*well_nseg; segID < well_nseg; ++segID) { this->seg_number_.push_back(segment_set[segID].segmentNumber()); @@ -665,10 +665,13 @@ namespace Opm // 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> segment_perforations(well_nseg); - for (int perf = 0; perf < nperf; ++perf) { + for (size_t perf = 0; perf < completion_set.size(); ++perf) { const Connection& connection = completion_set.get(perf); - const int segment_index = segment_set.segmentNumberToIndex(connection.segment()); - segment_perforations[segment_index].push_back(perf); + if (connection.state() == WellCompletion::OPEN) { + const int segment_index = segment_set.segmentNumberToIndex(connection.segment()); + segment_perforations[segment_index].push_back(n_activeperf); + n_activeperf++; + } } std::vector> segment_inlets(well_nseg); @@ -689,7 +692,10 @@ namespace Opm const int np = numPhases(); const int start_perf = wells->well_connpos[w]; const int start_perf_next_well = wells->well_connpos[w + 1]; - assert(nperf == (start_perf_next_well - start_perf)); // make sure the information from wells_ecl consistent with wells + + // make sure the information from wells_ecl consistent with wells + assert(n_activeperf == (start_perf_next_well - start_perf)); + if (pu.phase_used[Gas]) { const int gaspos = pu.phase_pos[Gas]; // scale the phase rates for Gas to avoid too bad initial guess for gas fraction @@ -697,7 +703,7 @@ namespace Opm // TODO: to see if this strategy can benefit StandardWell too // TODO: it might cause big problem for gas rate control or if there is a gas rate limit // maybe the best way is to initialize the fractions first then get the rates - for (int perf = 0; perf < nperf; perf++) { + for (int perf = 0; perf < n_activeperf; perf++) { const int perf_pos = start_perf + perf; perfPhaseRates()[np * perf_pos + gaspos] *= 100.; }