diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index fd0947a28..9aff8340f 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -53,6 +53,7 @@ namespace Opm tolerance_pressure_ms_wells_ = param.getDefault("tolerance_pressure_ms_wells", tolerance_pressure_ms_wells_); max_welleq_iter_ = param.getDefault("max_welleq_iter", max_welleq_iter_); max_pressure_change_ms_wells_ = param.getDefault("max_pressure_change_ms_wells", max_pressure_change_ms_wells_); + use_inner_iterations_ms_wells_ = param.getDefault("use_inner_iterations_ms_wells", use_inner_iterations_ms_wells_); maxSinglePrecisionTimeStep_ = unit::convert::from( param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day ); max_strict_iter_ = param.getDefault("max_strict_iter",8); @@ -81,7 +82,8 @@ namespace Opm tolerance_well_control_ = 1.0e-7; tolerance_pressure_ms_wells_ = 100.0; max_welleq_iter_ = 15; - max_pressure_change_ms_wells_ = 100000.; // 1.0 bar + max_pressure_change_ms_wells_ = 200000.; // 2.0 bar + use_inner_iterations_ms_wells_ = true; maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day ); solve_welleq_initially_ = true; update_equations_scaling_ = false; diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index af366c87d..7d577b8c5 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -56,6 +56,9 @@ namespace Opm /// Maximum pressure change over an iteratio for ms wells double max_pressure_change_ms_wells_; + /// Whether to use inner iterations for ms wells + bool use_inner_iterations_ms_wells_; + /// Maximum iteration number of the well equation solution int max_welleq_iter_; diff --git a/opm/autodiff/BlackoilWellModel.hpp b/opm/autodiff/BlackoilWellModel.hpp index 1610aaf07..3e7c591b0 100644 --- a/opm/autodiff/BlackoilWellModel.hpp +++ b/opm/autodiff/BlackoilWellModel.hpp @@ -212,7 +212,7 @@ namespace Opm { void computeRepRadiusPerfLength(const Grid& grid); - void computeAverageFormationFactor(Simulator& ebosSimulator, + void computeAverageFormationFactor(const Simulator& ebosSimulator, std::vector& B_avg) const; void applyVREPGroupControl(WellState& well_state) const; diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp index 0d658656c..4e232f39c 100644 --- a/opm/autodiff/BlackoilWellModel_impl.hpp +++ b/opm/autodiff/BlackoilWellModel_impl.hpp @@ -210,7 +210,7 @@ namespace Opm { bool only_wells) const { for (int w = 0; w < number_of_wells_; ++w) { - well_container_[w]->assembleWellEq(ebosSimulator, dt, well_state, only_wells); + well_container_[w]->assembleWellEq(ebosSimulator, param_, dt, well_state, only_wells); } } @@ -929,7 +929,7 @@ namespace Opm { template void BlackoilWellModel:: - computeAverageFormationFactor(Simulator& ebosSimulator, + computeAverageFormationFactor(const Simulator& ebosSimulator, std::vector& B_avg) const { const int np = numPhases(); diff --git a/opm/autodiff/MultisegmentWell.hpp b/opm/autodiff/MultisegmentWell.hpp index 0a2dcf3a5..48a9d5770 100644 --- a/opm/autodiff/MultisegmentWell.hpp +++ b/opm/autodiff/MultisegmentWell.hpp @@ -105,6 +105,7 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const; virtual void assembleWellEq(Simulator& ebosSimulator, + const ModelParameters& param, const double dt, WellState& well_state, bool only_wells); @@ -115,7 +116,7 @@ namespace Opm WellState& well_state) const; /// check whether the well equations get converged for this well - virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, + virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator, const std::vector& B_avg, const ModelParameters& param) const; @@ -327,6 +328,17 @@ namespace Opm bool frictionalPressureLossConsidered() const; bool accelerationalPressureLossConsidered() const; + + // TODO: try to make ebosSimulator const, as it should be + void iterateWellEquations(Simulator& ebosSimulator, + const ModelParameters& param, + const double dt, + WellState& well_state); + + void assembleWellEqWithoutIteration(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells); }; } diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index e12396ff7..7f69bb647 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -220,120 +220,18 @@ namespace Opm void MultisegmentWell:: assembleWellEq(Simulator& ebosSimulator, + const ModelParameters& param, const double dt, WellState& well_state, bool only_wells) { - // clear all entries - if (!only_wells) { - duneB_ = 0.0; - duneC_ = 0.0; + + const bool use_inner_iterations = param.use_inner_iterations_ms_wells_; + if (use_inner_iterations) { + iterateWellEquations(ebosSimulator, param, dt, well_state); } - duneD_ = 0.0; - resWell_ = 0.0; - - // for the black oil cases, there will be four equations, - // the first three of them are the mass balance equations, the last one is the pressure equations. - // - // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same. - - auto& ebosJac = ebosSimulator.model().linearizer().matrix(); - auto& ebosResid = ebosSimulator.model().linearizer().residual(); - - const bool allow_cf = getAllowCrossFlow(); - - const int nseg = numberOfSegments(); - const int num_comp = numComponents(); - - // TODO: finding better place to put it - computeSegmentFluidProperties(ebosSimulator); - - for (int seg = 0; seg < nseg; ++seg) { - // calculating the accumulation term // TODO: without considering the efficiencty factor for now - // volume of the segment - { - const double volume = segmentSet()[seg].volume(); - // for each component - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { - EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx]) - + getSegmentRate(seg, comp_idx); - - resWell_[seg][comp_idx] += accumulation_term.value(); - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq); - } - } - } - - // considering the contributions from the inlet segments - { - for (const int inlet : segment_inlets_[seg]) { - for (int comp_idx = 0; comp_idx < num_comp; ++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); - } - } - } - } - - // calculating the perforation rate for each perforation that belongs to this segment - const EvalWell seg_pressure = getSegmentPressure(seg); - for (const int perf : segment_perforations_[seg]) { - const int cell_idx = well_cells_[perf]; - const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); - std::vector mob(num_comp, 0.0); - getMobility(ebosSimulator, perf, mob); - std::vector cq_s(num_comp, 0.0); - computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s); - - for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { - // the cq_s entering mass balance equations need to consider the efficiency factors. - const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_; - - if (!only_wells) { - // subtract sum of component fluxes in the reservoir equation. - // need to consider the efficiency factor - // TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input - // is a component index :D - ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value(); - } - - // subtract sum of phase fluxes in the well equations. - resWell_[seg][comp_idx] -= cq_s_effective.value(); - - // assemble the jacobians - for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(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); - } - - for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { - if (!only_wells) { - // also need to consider the efficiency factor when manipulating the jacobians. - ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(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 - if (seg == 0) { // top segment, pressure equation is the control equation - assembleControlEq(); - } else { - assemblePressureEq(seg); - } - } + assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, only_wells); } @@ -491,7 +389,7 @@ namespace Opm template typename MultisegmentWell::ConvergenceReport MultisegmentWell:: - getWellConvergence(Simulator& ebosSimulator, + getWellConvergence(const Simulator& /* ebosSimulator */, const std::vector& B_avg, const ModelParameters& param) const { @@ -823,6 +721,10 @@ namespace Opm const BlackoilModelParameters& param, WellState& well_state) const { + const bool use_inner_iterations = param.use_inner_iterations_ms_wells_; + + const double relaxation_factor = use_inner_iterations ? 0.4 : 1.0; + // I guess the following can also be applied to the segmnet pressure // maybe better to give it a different name const double dBHPLimit = param.dbhp_max_rel_; @@ -833,13 +735,13 @@ namespace Opm for (int seg = 0; seg < numberOfSegments(); ++seg) { if (active()[ Water ]) { const int sign = dwells[seg][WFrac] > 0. ? 1 : -1; - const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]), 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 (active()[ Gas ]) { const int sign = dwells[seg][GFrac] > 0. ? 1 : -1; - const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]), 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; } @@ -856,7 +758,7 @@ namespace Opm // update the total rate // TODO: should we have a limitation of the total rate change? { - primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - dwells[seg][GTotal]; + primary_variables_[seg][GTotal] = old_primary_variables[seg][GTotal] - relaxation_factor * dwells[seg][GTotal]; } // TODO: not handling solvent related for now @@ -1860,4 +1762,177 @@ namespace Opm return (segmentSet().compPressureDrop() == WellSegment::HFA); } + + + + + template + void + MultisegmentWell:: + iterateWellEquations(Simulator& ebosSimulator, + const ModelParameters& param, + const double dt, + WellState& well_state) + { + std::cout << " beginning iterateWellEquations " << std::endl; + // 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. + // assembleWellEq( + const int max_iter_number = 3; + int it = 0; + for (; it < max_iter_number; ++it) { + + assembleWellEqWithoutIteration(ebosSimulator, dt, well_state, true); + + const BVectorWell dx_well = mswellhelpers::invDX(duneD_, resWell_); + + // TODO: use these small values for now, not intend to reach the convergence + // in this stage, but, should we? + // If we want to use the real one, we need to find a way to get them. + const std::vector B {0.1, 0.1, 0.001}; + + const ConvergenceReport report = getWellConvergence(ebosSimulator, B, param); + + if (report.converged) { + std::cout << " converged in iterateWellEquations " << std::endl; + break; + } + + updateWellState(dx_well, param, well_state); + + initPrimaryVariablesEvaluation(); + } + + if (it >= max_iter_number) { + std::cout << " the iterateWellEquations did not converged " << std::endl; + } + } + + + + + + template + void + MultisegmentWell:: + assembleWellEqWithoutIteration(Simulator& ebosSimulator, + const double dt, + WellState& well_state, + bool only_wells) + { + // calculate the fluid properties needed. + computeSegmentFluidProperties(ebosSimulator); + + // clear all entries + if (!only_wells) { + duneB_ = 0.0; + duneC_ = 0.0; + } + + duneD_ = 0.0; + resWell_ = 0.0; + + // for the black oil cases, there will be four equations, + // the first three of them are the mass balance equations, the last one is the pressure equations. + // + // but for the top segment, the pressure equation will be the well control equation, and the other three will be the same. + + auto& ebosJac = ebosSimulator.model().linearizer().matrix(); + auto& ebosResid = ebosSimulator.model().linearizer().residual(); + + const bool allow_cf = getAllowCrossFlow(); + + const int nseg = numberOfSegments(); + const int num_comp = numComponents(); + + for (int seg = 0; seg < nseg; ++seg) { + // calculating the accumulation term // TODO: without considering the efficiencty factor for now + // volume of the segment + { + const double volume = segmentSet()[seg].volume(); + // for each component + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + EvalWell accumulation_term = volume / dt * (surfaceVolumeFraction(seg, comp_idx) - segment_comp_initial_[seg][comp_idx]) + + getSegmentRate(seg, comp_idx); + + resWell_[seg][comp_idx] += accumulation_term.value(); + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + duneD_[seg][seg][comp_idx][pv_idx] += accumulation_term.derivative(pv_idx + numEq); + } + } + } + + // considering the contributions from the inlet segments + { + for (const int inlet : segment_inlets_[seg]) { + for (int comp_idx = 0; comp_idx < num_comp; ++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); + } + } + } + } + + // calculating the perforation rate for each perforation that belongs to this segment + const EvalWell seg_pressure = getSegmentPressure(seg); + for (const int perf : segment_perforations_[seg]) { + const int cell_idx = well_cells_[perf]; + const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + std::vector mob(num_comp, 0.0); + getMobility(ebosSimulator, perf, mob); + std::vector cq_s(num_comp, 0.0); + computePerfRate(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s); + + for (int comp_idx = 0; comp_idx < num_comp; ++comp_idx) { + // the cq_s entering mass balance equations need to consider the efficiency factors. + const EvalWell cq_s_effective = cq_s[comp_idx] * well_efficiency_factor_; + + if (!only_wells) { + // subtract sum of component fluxes in the reservoir equation. + // need to consider the efficiency factor + // TODO: the name of the function flowPhaseToEbosCompIdx is prolematic, since the input + // is a component index :D + ebosResid[cell_idx][flowPhaseToEbosCompIdx(comp_idx)] -= cq_s_effective.value(); + } + + // subtract sum of phase fluxes in the well equations. + resWell_[seg][comp_idx] -= cq_s_effective.value(); + + // assemble the jacobians + for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + duneC_[seg][cell_idx][pv_idx][flowPhaseToEbosCompIdx(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); + } + + for (int pv_idx = 0; pv_idx < numEq; ++pv_idx) { + if (!only_wells) { + // also need to consider the efficiency factor when manipulating the jacobians. + ebosJac[cell_idx][cell_idx][flowPhaseToEbosCompIdx(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 + if (seg == 0) { // top segment, pressure equation is the control equation + assembleControlEq(); + } else { + assemblePressureEq(seg); + } + } + } + } diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index ffa9adb42..f1744decc 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -111,6 +111,7 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const; virtual void assembleWellEq(Simulator& ebosSimulator, + const ModelParameters& param, const double dt, WellState& well_state, bool only_wells); @@ -121,7 +122,7 @@ namespace Opm WellState& xw) const; /// check whether the well equations get converged for this well - virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, + virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator, const std::vector& B_avg, const ModelParameters& param) const; diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index f42159cc4..7873d0895 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -514,6 +514,7 @@ namespace Opm void StandardWell:: assembleWellEq(Simulator& ebosSimulator, + const ModelParameters& /* param */, const double dt, WellState& well_state, bool only_wells) @@ -1364,7 +1365,7 @@ namespace Opm template typename StandardWell::ConvergenceReport StandardWell:: - getWellConvergence(Simulator& ebosSimulator, + getWellConvergence(const Simulator& /* ebosSimulator */, const std::vector& B_avg, const ModelParameters& param) const { diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp index e0fa8924f..901c7b3c9 100644 --- a/opm/autodiff/WellInterface.hpp +++ b/opm/autodiff/WellInterface.hpp @@ -144,7 +144,7 @@ namespace Opm } }; - virtual ConvergenceReport getWellConvergence(Simulator& ebosSimulator, + virtual ConvergenceReport getWellConvergence(const Simulator& ebosSimulator, const std::vector& B_avg, const ModelParameters& param) const = 0; @@ -152,6 +152,7 @@ namespace Opm WellState& well_state) = 0; virtual void assembleWellEq(Simulator& ebosSimulator, + const ModelParameters& param, const double dt, WellState& well_state, bool only_wells) = 0;