24 #ifndef OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
25 #define OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
27 #include <opm/autodiff/BlackoilModelBase.hpp>
28 #include <opm/autodiff/BlackoilDetails.hpp>
29 #include <opm/autodiff/BlackoilLegacyDetails.hpp>
31 #include <opm/autodiff/AutoDiffBlock.hpp>
32 #include <opm/autodiff/AutoDiffHelpers.hpp>
33 #include <opm/autodiff/GridHelpers.hpp>
34 #include <opm/autodiff/WellHelpers.hpp>
35 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
36 #include <opm/autodiff/GeoProps.hpp>
37 #include <opm/autodiff/WellDensitySegmented.hpp>
38 #include <opm/autodiff/VFPProperties.hpp>
39 #include <opm/autodiff/VFPProdProperties.hpp>
40 #include <opm/autodiff/VFPInjProperties.hpp>
42 #include <opm/core/grid.h>
43 #include <opm/core/linalg/LinearSolverInterface.hpp>
44 #include <opm/core/linalg/ParallelIstlInformation.hpp>
45 #include <opm/core/props/rock/RockCompressibility.hpp>
46 #include <opm/common/ErrorMacros.hpp>
47 #include <opm/common/Exceptions.hpp>
48 #include <opm/common/OpmLog/OpmLog.hpp>
49 #include <opm/parser/eclipse/Units/Units.hpp>
50 #include <opm/core/well_controls.h>
51 #include <opm/core/utility/parameters/ParameterGroup.hpp>
52 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
53 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
55 #include <dune/common/timer.hh>
67 #define OPM_AD_DUMP(foo) \
69 std::cout << "==========================================\n" \
71 << collapseJacs(foo) << std::endl; \
74 #define OPM_AD_DUMPVAL(foo) \
76 std::cout << "==========================================\n" \
78 << foo.value() << std::endl; \
81 #define OPM_AD_DISKVAL(foo) \
83 std::ofstream os(#foo); \
85 os << foo.value() << std::endl; \
91 typedef AutoDiffBlock<double> ADB;
94 typedef Eigen::Array<double,
97 Eigen::RowMajor> DataBlock;
99 template <
class Gr
id,
class WellModel,
class Implementation>
105 const RockCompressibility* rock_comp_props,
106 const WellModel& well_model,
108 std::shared_ptr< const Opm::EclipseState > eclState,
109 const bool has_disgas,
110 const bool has_vapoil,
111 const bool terminal_output)
115 , rock_comp_props_(rock_comp_props)
117 eclState->getTableManager().getVFPInjTables(),
118 eclState->getTableManager().getVFPProdTables())
119 , linsolver_ (linsolver)
120 , active_(detail::activePhases(fluid.phaseUsage()))
121 , canph_ (detail::active2Canonical(fluid.phaseUsage()))
122 , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
123 , ops_ (grid, geo.nnc())
124 , has_disgas_(has_disgas)
125 , has_vapoil_(has_vapoil)
127 , use_threshold_pressure_(false)
128 , sd_ (fluid.numPhases())
129 , phaseCondition_(AutoDiffGrid::numCells(grid))
130 , well_model_ (well_model)
131 , isRs_(V::Zero(AutoDiffGrid::numCells(grid)))
132 , isRv_(V::Zero(AutoDiffGrid::numCells(grid)))
133 , isSg_(V::Zero(AutoDiffGrid::numCells(grid)))
137 { 1.1169, 1.0031, 0.0031 },
139 , terminal_output_ (terminal_output)
141 , current_relaxation_(1.0)
146 , rate_converter_(fluid_.phaseUsage(), std::vector<int>(AutoDiffGrid::numCells(grid_),0))
148 if (active_[Water]) {
149 material_name_.push_back(
"Water");
152 material_name_.push_back(
"Oil");
155 material_name_.push_back(
"Gas");
158 assert(numMaterials() == std::accumulate(active_.begin(), active_.end(), 0));
160 const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
161 const V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_);
163 well_model_.init(&fluid_, &active_, &phaseCondition_, &vfp_properties_, gravity, depth);
168 const Wells* wells_arg = asImpl().well_model_.wellsPointer();
170 if ( linsolver_.parallelInformation().type() ==
typeid(ParallelISTLInformation) )
172 const ParallelISTLInformation& info =
173 boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation());
174 if ( terminal_output_ ) {
176 terminal_output_ = (info.communicator().rank()==0);
178 int local_number_of_wells = localWellsActive() ? wells().number_of_wells : 0;
179 int global_number_of_wells = info.communicator().sum(local_number_of_wells);
180 const bool wells_active = ( wells_arg && global_number_of_wells > 0 );
181 wellModel().setWellsActive(wells_active);
183 std::vector<int> v( Opm::AutoDiffGrid::numCells(grid_), 1);
185 info.computeReduction(v, Opm::Reduction::makeGlobalSumFunctor<int>(), global_nc_);
189 wellModel().setWellsActive( localWellsActive() );
190 global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
195 template <
class Gr
id,
class WellModel,
class Implementation>
201 if ( linsolver_.parallelInformation().type() !=
202 typeid(ParallelISTLInformation) )
208 const auto& comm =boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation()).communicator();
209 return comm.size() > 1;
217 template <
class Gr
id,
class WellModel,
class Implementation>
221 const ReservoirState& reservoir_state,
226 pvdt_ = geo_.poreVolume() / dt;
228 updatePrimalVariableFromState(reservoir_state);
236 template <
class Gr
id,
class WellModel,
class Implementation>
237 template <
class NonlinearSolverType>
242 NonlinearSolverType& nonlinear_solver,
243 ReservoirState& reservoir_state,
244 WellState& well_state)
246 SimulatorReport report;
247 Dune::Timer perfTimer;
252 if (iteration == 0) {
255 residual_norms_history_.clear();
256 current_relaxation_ = 1.0;
257 dx_old_ = V::Zero(sizeNonLinear());
260 report += asImpl().assemble(reservoir_state, well_state, iteration == 0);
261 report.assemble_time += perfTimer.stop();
264 report.assemble_time += perfTimer.stop();
268 report.total_linearizations = 1;
271 report.converged = asImpl().getConvergence(timer, iteration);
272 residual_norms_history_.push_back(asImpl().computeResidualNorms());
273 report.update_time += perfTimer.stop();
275 const bool must_solve = (iteration < nonlinear_solver.minIter()) || (!report.converged);
279 report.total_newton_iterations = 1;
282 residual_.singlePrecision = ( dt < param_.maxSinglePrecisionTimeStep_ );
287 dx = asImpl().solveJacobianSystem();
288 report.linear_solve_time += perfTimer.stop();
289 report.total_linear_iterations += linearIterationsLastSolve();
292 report.linear_solve_time += perfTimer.stop();
293 report.total_linear_iterations += linearIterationsLastSolve();
300 if (param_.use_update_stabilization_) {
302 bool isOscillate =
false;
303 bool isStagnate =
false;
304 nonlinear_solver.detectOscillations(residual_norms_history_, iteration, isOscillate, isStagnate);
306 current_relaxation_ -= nonlinear_solver.relaxIncrement();
307 current_relaxation_ = std::max(current_relaxation_, nonlinear_solver.relaxMax());
308 if (terminalOutputEnabled()) {
309 std::string msg =
" Oscillating behavior detected: Relaxation set to "
310 + std::to_string(current_relaxation_);
314 nonlinear_solver.stabilizeNonlinearUpdate(dx, dx_old_, current_relaxation_);
319 asImpl().updateState(dx, reservoir_state, well_state);
320 report.update_time += perfTimer.stop();
330 template <
class Gr
id,
class WellModel,
class Implementation>
344 template <
class Gr
id,
class WellModel,
class Implementation>
349 return residual_.sizeNonLinear();
356 template <
class Gr
id,
class WellModel,
class Implementation>
361 return linsolver_.iterations();
368 template <
class Gr
id,
class WellModel,
class Implementation>
373 return terminal_output_;
380 template <
class Gr
id,
class WellModel,
class Implementation>
385 return fluid_.numPhases();
391 template <
class Gr
id,
class WellModel,
class Implementation>
396 return material_name_.size();
403 template <
class Gr
id,
class WellModel,
class Implementation>
408 assert(material_index < numMaterials());
409 return material_name_[material_index];
416 template <
class Gr
id,
class WellModel,
class Implementation>
421 const int num_faces = AutoDiffGrid::numFaces(grid_);
422 const int num_nnc = geo_.nnc().numNNC();
423 const int num_connections = num_faces + num_nnc;
424 if (
int(threshold_pressures.size()) != num_connections) {
425 OPM_THROW(std::runtime_error,
"Illegal size of threshold_pressures input ( " << threshold_pressures.size()
426 <<
" ), must be equal to number of faces + nncs ( " << num_faces <<
" + " << num_nnc <<
" ).");
428 use_threshold_pressure_ =
true;
430 const int num_ifaces = ops_.internal_faces.size();
431 threshold_pressures_by_connection_.resize(num_ifaces + num_nnc);
432 for (
int ii = 0; ii < num_ifaces; ++ii) {
433 threshold_pressures_by_connection_[ii] = threshold_pressures[ops_.internal_faces[ii]];
437 for (
int ii = 0; ii < num_nnc; ++ii) {
438 threshold_pressures_by_connection_[ii + num_ifaces] = threshold_pressures[ii + num_faces];
446 template <
class Gr
id,
class WellModel,
class Implementation>
449 : accum(2,
ADB::null())
450 , mflux(
ADB::null())
460 template <
class Gr
id,
class WellModel,
class Implementation>
461 BlackoilModelBase<Grid, WellModel, Implementation>::
462 SimulatorData::SimulatorData(
int num_phases)
481 template <
class Gr
id,
class WellModel,
class Implementation>
483 BlackoilModelBase<Grid, WellModel, Implementation>::
484 makeConstantState(SolutionState& state)
const
491 state.temperature =
ADB::constant(state.temperature.value());
494 const int num_phases = state.saturation.
size();
495 for (
int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) {
496 state.saturation[phaseIdx] =
ADB::constant(state.saturation[phaseIdx].value());
500 assert(state.canonical_phase_pressures.size() ==
static_cast<std::size_t
>(Opm::BlackoilPhases::MaxNumPhases));
501 for (
int canphase = 0; canphase < Opm::BlackoilPhases::MaxNumPhases; ++canphase) {
502 ADB& pp = state.canonical_phase_pressures[canphase];
511 template <
class Gr
id,
class WellModel,
class Implementation>
512 typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
513 BlackoilModelBase<Grid, WellModel, Implementation>::
514 variableState(
const ReservoirState& x,
515 const WellState& xw)
const
517 std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
519 return asImpl().variableStateExtractVars(x, asImpl().variableStateIndices(), vars);
526 template <
class Gr
id,
class WellModel,
class Implementation>
528 BlackoilModelBase<Grid, WellModel, Implementation>::
529 variableStateInitials(
const ReservoirState& x,
530 const WellState& xw)
const
532 assert(active_[ Oil ]);
534 const int np = x.numPhases();
536 std::vector<V> vars0;
539 vars0.reserve(np + 1);
540 variableReservoirStateInitials(x, vars0);
541 asImpl().wellModel().variableWellStateInitials(xw, vars0);
549 template <
class Gr
id,
class WellModel,
class Implementation>
551 BlackoilModelBase<Grid, WellModel, Implementation>::
552 variableReservoirStateInitials(
const ReservoirState& x, std::vector<V>& vars0)
const
554 using namespace Opm::AutoDiffGrid;
555 const int nc = numCells(grid_);
556 const int np = x.numPhases();
558 assert (not x.pressure().empty());
559 const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
563 assert (not x.saturation().empty());
564 const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
565 const Opm::PhaseUsage pu = fluid_.phaseUsage();
567 assert (active_[ Oil]);
568 if (active_[ Water ]) {
569 const V sw = s.col(pu.phase_pos[ Water ]);
573 if (active_[ Gas ]) {
576 const V sg = s.col(pu.phase_pos[ Gas ]);
577 const V rs = Eigen::Map<const V>(& x.gasoilratio()[0], x.gasoilratio().size());
578 const V rv = Eigen::Map<const V>(& x.rv()[0], x.rv().size());
579 xvar = isRs_*rs + isRv_*rv + isSg_*sg;
580 vars0.push_back(xvar);
588 template <
class Gr
id,
class WellModel,
class Implementation>
590 BlackoilModelBase<Grid, WellModel, Implementation>::
591 variableStateIndices()
const
593 assert(active_[Oil]);
594 std::vector<int> indices(5, -1);
596 indices[Pressure] = next++;
597 if (active_[Water]) {
598 indices[Sw] = next++;
601 indices[Xvar] = next++;
603 asImpl().wellModel().variableStateWellIndices(indices, next);
604 assert(next == fluid_.numPhases() + 2);
612 template <
class Gr
id,
class WellModel,
class Implementation>
613 typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
614 BlackoilModelBase<Grid, WellModel, Implementation>::
615 variableStateExtractVars(
const ReservoirState& x,
616 const std::vector<int>& indices,
617 std::vector<ADB>& vars)
const
620 const int nc = Opm::AutoDiffGrid::numCells(grid_);
621 const Opm::PhaseUsage pu = fluid_.phaseUsage();
623 SolutionState state(fluid_.numPhases());
626 state.pressure = std::move(vars[indices[Pressure]]);
629 const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
636 if (active_[ Water ]) {
637 state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]);
638 const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
642 if (active_[ Gas ]) {
645 const ADB& xvar = vars[indices[Xvar]];
646 ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
647 sg = isSg_*xvar + isRv_*so;
652 const ADB& sw = (active_[ Water ]
653 ? state.saturation[ pu.phase_pos[ Water ] ]
655 state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
658 if (active_[ Oil ]) {
660 sd_.rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
662 state.rs = (1-isRs_)*sd_.rsSat + isRs_*xvar;
664 state.rs = sd_.rsSat;
666 sd_.rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
668 state.rv = (1-isRv_)*sd_.rvSat + isRv_*xvar;
670 state.rv = sd_.rvSat;
672 sd_.soMax = fluid_.satOilMax();
673 fluid_.getGasOilHystParams(sd_.krnswdc_go, sd_.pcswmdc_go, cells_);
674 fluid_.getOilWaterHystParams(sd_.krnswdc_ow, sd_.pcswmdc_ow, cells_);
676 sd_.Pb = fluid_.bubblePointPressure(cells_,
677 state.temperature.value(),
679 sd_.Pd = fluid_.dewPointPressure(cells_,
680 state.temperature.value(),
686 const ADB& sw = (active_[ Water ]
687 ? state.saturation[ pu.phase_pos[ Water ] ]
690 state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg);
693 if (active_[ Oil ]) {
695 state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
699 asImpl().wellModel().variableStateExtractWellsVars(indices, vars, state);
707 template <
class Gr
id,
class WellModel,
class Implementation>
709 BlackoilModelBase<Grid, WellModel, Implementation>::
710 computeAccum(
const SolutionState& state,
713 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
715 const ADB& press = state.pressure;
716 const ADB& temp = state.temperature;
717 const std::vector<ADB>& sat = state.saturation;
718 const ADB& rs = state.rs;
719 const ADB& rv = state.rv;
721 const std::vector<PhasePresence> cond = phaseCondition();
723 const ADB pv_mult = poroMult(press);
725 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
726 for (
int phase = 0; phase < maxnp; ++phase) {
727 if (active_[ phase ]) {
728 const int pos = pu.phase_pos[ phase ];
729 sd_.rq[pos].b = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
730 sd_.rq[pos].accum[aix] = pv_mult * sd_.rq[pos].b * sat[pos];
736 if (active_[ Oil ] && active_[ Gas ]) {
738 const int po = pu.phase_pos[ Oil ];
739 const int pg = pu.phase_pos[ Gas ];
743 const ADB accum_gas_copy =sd_.rq[pg].accum[aix];
745 sd_.rq[pg].accum[aix] += state.rs * sd_.rq[po].accum[aix];
746 sd_.rq[po].accum[aix] += state.rv * accum_gas_copy;
755 template <
class Gr
id,
class WellModel,
class Implementation>
759 WellState& well_state,
760 const bool initial_assembly)
762 using namespace Opm::AutoDiffGrid;
764 SimulatorReport report;
770 SolutionState state = asImpl().variableState(reservoir_state, well_state);
771 SolutionState state0 = state;
772 asImpl().makeConstantState(state0);
773 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
777 if (asImpl().wellModel().wellCollection()->groupControlActive() && initial_assembly) {
778 setupGroupControl(reservoir_state, well_state);
783 asImpl().wellModel().updateWellControls(well_state);
785 if (asImpl().wellModel().wellCollection()->groupControlActive()) {
787 applyVREPGroupControl(reservoir_state, well_state);
789 asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
793 SolutionState state = asImpl().variableState(reservoir_state, well_state);
795 if (initial_assembly) {
797 SolutionState state0 = state;
798 asImpl().makeConstantState(state0);
801 asImpl().computeAccum(state0, 0);
802 asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
815 asImpl().assembleMassBalanceEq(state);
818 if ( ! wellsActive() ) {
822 std::vector<ADB> mob_perfcells;
823 std::vector<ADB> b_perfcells;
824 asImpl().wellModel().extractWellPerfProperties(state, sd_.rq, mob_perfcells, b_perfcells);
825 if (param_.solve_welleq_initially_ && initial_assembly) {
827 report += asImpl().solveWellEq(mob_perfcells, b_perfcells, reservoir_state, state, well_state);
830 std::vector<ADB> cq_s;
831 asImpl().wellModel().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);
832 asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, state, well_state);
833 asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
834 asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
835 asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
843 template <
class Gr
id,
class WellModel,
class Implementation>
854 asImpl().computeAccum(state, 1);
858 const V transi =
subset(geo_.transmissibility(), ops_.internal_faces);
859 const V trans_nnc = ops_.nnc_trans;
860 V trans_all(transi.size() + trans_nnc.size());
861 trans_all << transi, trans_nnc;
865 const std::vector<ADB> kr = asImpl().computeRelPerm(state);
866 for (
int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
867 sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
870 #pragma omp parallel for schedule(static)
871 for (
int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
872 const std::vector<PhasePresence>& cond = phaseCondition();
873 sd_.rq[phaseIdx].mu = asImpl().fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
874 sd_.rq[phaseIdx].rho = asImpl().fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
875 asImpl().computeMassFlux(phaseIdx, trans_all, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
877 residual_.material_balance_eq[ phaseIdx ] =
878 pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
879 + ops_.div*sd_.rq[phaseIdx].mflux;
887 if (active_[ Oil ] && active_[ Gas ]) {
888 const int po = fluid_.phaseUsage().phase_pos[ Oil ];
889 const int pg = fluid_.phaseUsage().phase_pos[ Gas ];
891 const UpwindSelector<double> upwindOil(grid_, ops_,
892 sd_.rq[po].dh.value());
893 const ADB rs_face = upwindOil.select(state.rs);
895 const UpwindSelector<double> upwindGas(grid_, ops_,
896 sd_.rq[pg].dh.value());
897 const ADB rv_face = upwindGas.select(state.rv);
899 residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * sd_.rq[po].mflux);
900 residual_.material_balance_eq[ po ] += ops_.div * (rv_face * sd_.rq[pg].mflux);
907 if (param_.update_equations_scaling_) {
908 asImpl().updateEquationsScaling();
917 template <
class Gr
id,
class WellModel,
class Implementation>
922 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
923 for (
int idx=0; idx<MaxNumPhases; ++idx )
926 const int pos = pu.phase_pos[idx];
927 const ADB& temp_b = sd_.rq[pos].b;
928 B = 1. / temp_b.
value();
930 if ( linsolver_.parallelInformation().type() ==
typeid(ParallelISTLInformation) )
932 const ParallelISTLInformation& real_info =
933 boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation());
934 double B_global_sum = 0;
935 real_info.computeReduction(B, Reduction::makeGlobalSumFunctor<double>(), B_global_sum);
936 residual_.matbalscale[idx] = B_global_sum / global_nc_;
941 residual_.matbalscale[idx] = B.mean();
951 template <
class Gr
id,
class WellModel,
class Implementation>
955 const SolutionState&,
958 if ( !asImpl().localWellsActive() )
966 const int nc = Opm::AutoDiffGrid::numCells(grid_);
967 const int np = asImpl().numPhases();
968 const V& efficiency_factors = wellModel().wellPerfEfficiencyFactors();
969 for (
int phase = 0; phase < np; ++phase) {
970 residual_.material_balance_eq[phase] -=
superset(efficiency_factors * cq_s[phase],
971 wellModel().wellOps().well_cells, nc);
979 template <
class Gr
id,
class WellModel,
class Implementation>
981 BlackoilModelBase<Grid, WellModel, Implementation>::
984 if( ! localWellsActive() ) {
988 if ( vfp_properties_.getProd()->empty() && vfp_properties_.getInj()->empty() ) {
992 const int nw = wells().number_of_wells;
994 for (
int w = 0; w < nw; ++w) {
995 const WellControls* wc = wells().ctrls[w];
997 const int nwc = well_controls_get_num(wc);
1000 for (
int c=0; c < nwc; ++c) {
1001 const WellControlType ctrl_type = well_controls_iget_type(wc, c);
1003 if (ctrl_type == THP) {
1015 template <
class Gr
id,
class WellModel,
class Implementation>
1017 BlackoilModelBase<Grid, WellModel, Implementation>::
1018 solveWellEq(
const std::vector<ADB>& mob_perfcells,
1019 const std::vector<ADB>& b_perfcells,
1020 const ReservoirState& reservoir_state,
1021 SolutionState& state,
1022 WellState& well_state)
1025 const int np = wells().number_of_phases;
1027 std::vector<int> indices = asImpl().wellModel().variableWellStateIndices();
1028 SolutionState state0 = state;
1029 WellState well_state0 = well_state;
1030 asImpl().makeConstantState(state0);
1032 std::vector<ADB> mob_perfcells_const(np,
ADB::null());
1033 std::vector<ADB> b_perfcells_const(np,
ADB::null());
1035 if (asImpl().localWellsActive() ){
1038 for (
int phase = 0; phase < np; ++phase) {
1039 mob_perfcells_const[phase] =
ADB::constant(mob_perfcells[phase].value());
1040 b_perfcells_const[phase] =
ADB::constant(b_perfcells[phase].value());
1048 std::vector<V> vars0;
1050 asImpl().wellModel().variableWellStateInitials(well_state, vars0);
1053 SolutionState wellSolutionState = state0;
1054 asImpl().wellModel().variableStateExtractWellsVars(indices, vars, wellSolutionState);
1055 asImpl().wellModel().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s);
1056 asImpl().wellModel().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state);
1057 asImpl().wellModel().addWellFluxEq(cq_s, wellSolutionState, residual_);
1058 asImpl().wellModel().addWellControlEq(wellSolutionState, well_state, aliveWells, residual_);
1059 converged = getWellConvergence(it);
1066 if( localWellsActive() )
1068 std::vector<ADB> eqs;
1070 eqs.push_back(residual_.well_flux_eq);
1071 eqs.push_back(residual_.well_eq);
1073 const std::vector<M>& Jn = total_residual.derivative();
1074 typedef Eigen::SparseMatrix<double> Sp;
1076 Jn[0].toSparse(Jn0);
1077 const Eigen::SparseLU< Sp > solver(Jn0);
1078 ADB::V total_residual_v = total_residual.value();
1079 const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix());
1080 assert(dx.size() == total_residual_v.size());
1081 asImpl().wellModel().updateWellState(dx.array(), dbhpMaxRel(), well_state);
1086 asImpl().wellModel().updateWellControls(well_state);
1088 if (asImpl().wellModel().wellCollection()->groupControlActive()) {
1090 applyVREPGroupControl(reservoir_state, well_state);
1091 asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
1096 if (terminalOutputEnabled())
1098 OpmLog::note(
"well converged iter: " + std::to_string(it));
1101 const int nw = wells().number_of_wells;
1105 ADB::V new_bhp = Eigen::Map<ADB::V>(well_state.bhp().data(), nw);
1108 std::vector<ADB::M> old_derivs = state.bhp.derivative();
1109 state.bhp =
ADB::function(std::move(new_bhp), std::move(old_derivs));
1115 const DataBlock wrates = Eigen::Map<const DataBlock>(well_state.wellRates().data(), nw, np).transpose();
1116 ADB::V new_qs = Eigen::Map<const V>(wrates.data(), nw*np);
1117 std::vector<ADB::M> old_derivs = state.qs.derivative();
1118 state.qs =
ADB::function(std::move(new_qs), std::move(old_derivs));
1120 asImpl().computeWellConnectionPressures(state, well_state);
1124 well_state = well_state0;
1127 SimulatorReport report;
1128 report.total_well_iterations = it;
1129 report.converged = converged;
1137 template <
class Gr
id,
class WellModel,
class Implementation>
1142 return linsolver_.computeNewtonIncrement(residual_);
1145 template <
class Gr
id,
class WellModel,
class Implementation>
1149 ReservoirState& reservoir_state,
1150 WellState& well_state)
1152 using namespace Opm::AutoDiffGrid;
1153 const int np = fluid_.numPhases();
1154 const int nc = numCells(grid_);
1156 assert(null.size() == 0);
1157 const V zero = V::Zero(nc);
1158 const V ones = V::Constant(nc,1.0);
1163 const V dsw = active_[Water] ?
subset(dx,
Span(nc, 1, varstart)) : null;
1164 varstart += dsw.size();
1166 const V dxvar = active_[Gas] ?
subset(dx,
Span(nc, 1, varstart)): null;
1167 varstart += dxvar.size();
1170 const V dwells =
subset(dx,
Span(asImpl().wellModel().numWellVars(), 1, varstart));
1171 varstart += dwells.size();
1173 assert(varstart == dx.size());
1176 const double dpmaxrel = dpMaxRel();
1177 const V p_old = Eigen::Map<const V>(&reservoir_state.pressure()[0], nc, 1);
1178 const V absdpmax = dpmaxrel*p_old.abs();
1179 const V dp_limited =
sign(dp) * dp.abs().min(absdpmax);
1180 const V p = (p_old - dp_limited).max(zero);
1181 std::copy(&p[0], &p[0] + nc, reservoir_state.pressure().begin());
1184 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
1185 const DataBlock s_old = Eigen::Map<const DataBlock>(& reservoir_state.saturation()[0], nc, np);
1186 const double dsmax = dsMax();
1195 if (active_[Water]){
1196 maxVal = dsw.abs().max(maxVal);
1202 dsg = isSg_ * dxvar - isRv_ * dsw;
1203 maxVal = dsg.abs().max(maxVal);
1207 maxVal = dso.abs().max(maxVal);
1209 V step = dsmax/maxVal;
1210 step = step.min(1.);
1212 if (active_[Water]) {
1213 const int pos = pu.phase_pos[ Water ];
1214 const V sw_old = s_old.col(pos);
1215 sw = sw_old - step * dsw;
1219 const int pos = pu.phase_pos[ Gas ];
1220 const V sg_old = s_old.col(pos);
1221 sg = sg_old - step * dsg;
1224 assert(active_[Oil]);
1225 const int pos = pu.phase_pos[ Oil ];
1226 const V so_old = s_old.col(pos);
1227 so = so_old - step * dso;
1232 for (
int c = 0; c < nc; ++c) {
1234 if (active_[Water]) {
1235 sw[c] = sw[c] / (1-sg[c]);
1237 so[c] = so[c] / (1-sg[c]);
1245 for (
int c = 0; c < nc; ++c) {
1247 if (active_[Water]) {
1248 sw[c] = sw[c] / (1-so[c]);
1251 sg[c] = sg[c] / (1-so[c]);
1258 if (active_[Water]) {
1260 for (
int c = 0; c < nc; ++c) {
1262 so[c] = so[c] / (1-sw[c]);
1264 sg[c] = sg[c] / (1-sw[c]);
1272 const double drmaxrel = drMaxRel();
1275 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
1276 const V drs = isRs_ * dxvar;
1277 const V drs_limited =
sign(drs) * drs.abs().min( (rs_old.abs()*drmaxrel).max( ones*1.0));
1278 rs = rs_old - drs_limited;
1283 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
1284 const V drv = isRv_ * dxvar;
1285 const V drv_limited =
sign(drv) * drv.abs().min( (rv_old.abs()*drmaxrel).max( ones*1e-3));
1286 rv = rv_old - drv_limited;
1291 const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
1292 auto watOnly = sw > (1 - epsilon);
1295 std::vector<HydroCarbonState>& hydroCarbonState = reservoir_state.hydroCarbonState();
1296 std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
1299 const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_);
1300 const V rsSat = fluidRsSat(p, so, cells_);
1303 auto hasGas = (sg > 0 && isRs_ == 0);
1306 const V rs_old = Eigen::Map<const V>(&reservoir_state.gasoilratio()[0], nc);
1307 auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs_ == 1 ) && (rs_old > rsSat0 * (1-epsilon)) );
1308 auto useSg = watOnly || hasGas || gasVaporized;
1309 for (
int c = 0; c < nc; ++c) {
1318 hydroCarbonState[c] = HydroCarbonState::OilOnly;
1328 const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas));
1329 const V gaspress = computeGasPressure(p, sw, so, sg);
1330 const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_);
1331 const V rvSat = fluidRvSat(gaspress, so, cells_);
1335 auto hasOil = (so > 0 && isRv_ == 0);
1338 const V rv_old = Eigen::Map<const V>(&reservoir_state.rv()[0], nc);
1339 auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv_ == 1) && (rv_old > rvSat0 * (1-epsilon)) );
1340 auto useSg = watOnly || hasOil || oilCondensed;
1341 for (
int c = 0; c < nc; ++c) {
1350 hydroCarbonState[c] = HydroCarbonState::GasOnly;
1357 if (active_[Water]) {
1358 for (
int c = 0; c < nc; ++c) {
1359 reservoir_state.saturation()[c*np + pu.phase_pos[ Water ]] = sw[c];
1364 for (
int c = 0; c < nc; ++c) {
1365 reservoir_state.saturation()[c*np + pu.phase_pos[ Gas ]] = sg[c];
1369 if (active_[ Oil ]) {
1370 for (
int c = 0; c < nc; ++c) {
1371 reservoir_state.saturation()[c*np + pu.phase_pos[ Oil ]] = so[c];
1376 std::copy(&rs[0], &rs[0] + nc, reservoir_state.gasoilratio().begin());
1380 std::copy(&rv[0], &rv[0] + nc, reservoir_state.rv().begin());
1383 asImpl().wellModel().updateWellState(dwells, dbhpMaxRel(), well_state);
1386 updatePhaseCondFromPrimalVariable(reservoir_state);
1393 template <
class Gr
id,
class WellModel,
class Implementation>
1398 using namespace Opm::AutoDiffGrid;
1399 const int nc = numCells(grid_);
1403 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
1404 const ADB& sw = (active_[ Water ]
1405 ? state.saturation[ pu.phase_pos[ Water ] ]
1408 const ADB& so = (active_[ Oil ]
1409 ? state.saturation[ pu.phase_pos[ Oil ] ]
1412 const ADB& sg = (active_[ Gas ]
1413 ? state.saturation[ pu.phase_pos[ Gas ] ]
1416 return fluid_.relperm(sw, so, sg, cells_);
1423 template <
class Gr
id,
class WellModel,
class Implementation>
1425 BlackoilModelBase<Grid, WellModel, Implementation>::
1426 computePressures(
const ADB& po,
1429 const ADB& sg)
const
1432 std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
1433 for (
int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
1435 if (phaseIdx == BlackoilPhases::Liquid)
1437 if (active_[phaseIdx]) {
1438 pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
1447 for (
int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
1448 if (active_[phaseIdx]) {
1449 if (phaseIdx == BlackoilPhases::Aqua) {
1450 pressure[phaseIdx] = po - pressure[phaseIdx];
1452 pressure[phaseIdx] += po;
1464 template <
class Gr
id,
class WellModel,
class Implementation>
1466 BlackoilModelBase<Grid, WellModel, Implementation>::
1467 computeGasPressure(
const V& po,
1472 assert (active_[Gas]);
1477 return cp[Gas].value() + po;
1482 template <
class Gr
id,
class WellModel,
class Implementation>
1484 BlackoilModelBase<Grid, WellModel, Implementation>::
1485 computeMassFlux(
const int actph ,
1490 const ADB& phasePressure,
1491 const SolutionState& state)
1494 const ADB tr_mult = transMult(state.pressure);
1495 sd_.rq[ actph ].mob = tr_mult * kr / mu;
1498 const ADB rhoavg = ops_.caver * rho;
1499 sd_.rq[ actph ].dh = ops_.ngrad * phasePressure - geo_.gravity()[2] * (rhoavg * (ops_.ngrad * geo_.z().matrix()));
1500 if (use_threshold_pressure_) {
1501 applyThresholdPressures(sd_.rq[ actph ].dh);
1505 const ADB& b = sd_.rq[ actph ].b;
1506 const ADB& mob = sd_.rq[ actph ].mob;
1507 const ADB& dh = sd_.rq[ actph ].dh;
1508 UpwindSelector<double> upwind(grid_, ops_, dh.value());
1509 sd_.rq[ actph ].mflux = upwind.select(b * mob) * (transi * dh);
1516 template <
class Gr
id,
class WellModel,
class Implementation>
1518 BlackoilModelBase<Grid, WellModel, Implementation>::
1519 applyThresholdPressures(ADB& dp)
1532 const V high_potential = (dp.value().abs() >= threshold_pressures_by_connection_).
template cast<double>();
1535 const M keep_high_potential(high_potential.matrix().asDiagonal());
1538 const V sign_dp =
sign(dp.value());
1539 const V threshold_modification = sign_dp * threshold_pressures_by_connection_;
1542 dp = keep_high_potential * (dp - threshold_modification);
1549 template <
class Gr
id,
class WellModel,
class Implementation>
1554 std::vector<double> residualNorms;
1556 std::vector<ADB>::const_iterator massBalanceIt = residual_.material_balance_eq.begin();
1557 const std::vector<ADB>::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
1559 for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
1560 const double massBalanceResid = detail::infinityNorm( (*massBalanceIt),
1561 linsolver_.parallelInformation() );
1562 if (!std::isfinite(massBalanceResid)) {
1563 OPM_THROW(Opm::NumericalProblem,
1564 "Encountered a non-finite residual");
1566 residualNorms.push_back(massBalanceResid);
1570 const double wellFluxResid = detail::infinityNormWell( residual_.well_flux_eq,
1571 linsolver_.parallelInformation() );
1572 if (!std::isfinite(wellFluxResid)) {
1573 OPM_THROW(Opm::NumericalProblem,
1574 "Encountered a non-finite residual");
1576 residualNorms.push_back(wellFluxResid);
1578 const double wellResid = detail::infinityNormWell( residual_.well_eq,
1579 linsolver_.parallelInformation() );
1580 if (!std::isfinite(wellResid)) {
1581 OPM_THROW(Opm::NumericalProblem,
1582 "Encountered a non-finite residual");
1584 residualNorms.push_back(wellResid);
1586 return residualNorms;
1590 template <
class Gr
id,
class WellModel,
class Implementation>
1594 const SimulationDataContainer& current )
const
1596 std::vector< double > p0 ( previous.pressure() );
1597 std::vector< double > sat0( previous.saturation() );
1599 const std::size_t pSize = p0.size();
1600 const std::size_t satSize = sat0.size();
1603 for( std::size_t i=0; i<pSize; ++i ) {
1604 p0[ i ] -= current.pressure()[ i ];
1607 for( std::size_t i=0; i<satSize; ++i ) {
1608 sat0[ i ] -= current.saturation()[ i ];
1612 const double stateOld = detail::euclidianNormSquared( p0.begin(), p0.end(), 1, linsolver_.parallelInformation() ) +
1613 detail::euclidianNormSquared( sat0.begin(), sat0.end(),
1614 current.numPhases(),
1615 linsolver_.parallelInformation() );
1618 const double stateNew = detail::euclidianNormSquared( current.pressure().begin(), current.pressure().end(), 1, linsolver_.parallelInformation() ) +
1619 detail::euclidianNormSquared( current.saturation().begin(), current.saturation().end(),
1620 current.numPhases(),
1621 linsolver_.parallelInformation() );
1623 if( stateNew > 0.0 ) {
1624 return stateOld / stateNew ;
1631 template <
class Gr
id,
class WellModel,
class Implementation>
1635 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
1636 const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
1637 std::vector<double>& R_sum,
1638 std::vector<double>& maxCoeff,
1639 std::vector<double>& B_avg,
1640 std::vector<double>& maxNormWell,
1643 const int np = asImpl().numPhases();
1644 const int nm = asImpl().numMaterials();
1645 const int nw = residual_.well_flux_eq.size() / np;
1646 assert(nw * np ==
int(residual_.well_flux_eq.size()));
1650 if ( linsolver_.parallelInformation().type() ==
typeid(ParallelISTLInformation) )
1652 const ParallelISTLInformation& info =
1653 boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation());
1656 std::vector<int> v(nc, 1);
1657 auto nc_and_pv = std::tuple<int, double>(0, 0.0);
1658 auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<int>(),
1659 Opm::Reduction::makeGlobalSumFunctor<double>());
1660 auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume());
1661 info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv);
1663 for (
int idx = 0; idx < nm; ++idx )
1665 auto values = std::tuple<double,double,double>(0.0 ,0.0 ,0.0);
1666 auto containers = std::make_tuple(B.col(idx),
1669 auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
1670 Opm::Reduction::makeGlobalMaxFunctor<double>(),
1671 Opm::Reduction::makeGlobalSumFunctor<double>());
1672 info.computeReduction(containers, operators, values);
1673 B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
1674 maxCoeff[idx] = std::get<1>(values);
1675 R_sum[idx] = std::get<2>(values);
1678 maxNormWell[idx] = 0.0;
1679 for (
int w = 0; w < nw; ++w ) {
1680 maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
1684 info.communicator().max(maxNormWell.data(), np);
1686 return std::get<1>(nc_and_pv);
1692 maxCoeff.resize(nm);
1694 maxNormWell.resize(np);
1695 for (
int idx = 0; idx < nm; ++idx )
1697 B_avg[idx] = B.col(idx).sum()/nc;
1698 maxCoeff[idx] = tempV.col(idx).maxCoeff();
1699 R_sum[idx] = R.col(idx).sum();
1703 maxNormWell[idx] = 0.0;
1704 for (
int w = 0; w < nw; ++w ) {
1705 maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
1710 return geo_.poreVolume().sum();
1718 template <
class Gr
id,
class WellModel,
class Implementation>
1724 const double tol_mb = param_.tolerance_mb_;
1725 const double tol_cnv = param_.tolerance_cnv_;
1726 const double tol_wells = param_.tolerance_wells_;
1727 const double tol_well_control = param_.tolerance_well_control_;
1729 const int nc = Opm::AutoDiffGrid::numCells(grid_);
1730 const int np = asImpl().numPhases();
1731 const int nm = asImpl().numMaterials();
1732 assert(
int(sd_.rq.size()) == nm);
1734 const V& pv = geo_.poreVolume();
1736 std::vector<double> R_sum(nm);
1737 std::vector<double> B_avg(nm);
1738 std::vector<double> maxCoeff(nm);
1739 std::vector<double> maxNormWell(np);
1740 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
1741 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
1742 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
1744 for (
int idx = 0; idx < nm; ++idx )
1746 const ADB& tempB = sd_.rq[idx].b;
1747 B.col(idx) = 1./tempB.
value();
1748 R.col(idx) = residual_.material_balance_eq[idx].value();
1749 tempV.col(idx) = R.col(idx).abs()/pv;
1752 const double pvSum = convergenceReduction(B, tempV, R,
1753 R_sum, maxCoeff, B_avg, maxNormWell,
1756 std::vector<double> CNV(nm);
1757 std::vector<double> mass_balance_residual(nm);
1758 std::vector<double> well_flux_residual(np);
1760 bool converged_MB =
true;
1761 bool converged_CNV =
true;
1762 bool converged_Well =
true;
1764 for (
int idx = 0; idx < nm; ++idx )
1766 CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
1767 mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
1768 converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
1769 converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
1774 well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
1775 converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
1779 const double residualWell = detail::infinityNormWell(residual_.well_eq,
1780 linsolver_.parallelInformation());
1781 converged_Well = converged_Well && (residualWell < tol_well_control);
1783 const bool converged = converged_MB && converged_CNV && converged_Well;
1786 const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed();
1788 if ( terminal_output_ )
1791 if (iteration == 0) {
1792 std::string msg =
"Iter";
1793 for (
int idx = 0; idx < nm; ++idx) {
1794 msg +=
" MB(" + materialName(idx).substr(0, 3) +
") ";
1796 for (
int idx = 0; idx < nm; ++idx) {
1797 msg +=
" CNV(" + materialName(idx).substr(0, 1) +
") ";
1799 for (
int idx = 0; idx < np; ++idx) {
1800 msg +=
" W-FLUX(" + materialName(idx).substr(0, 1) +
")";
1802 msg +=
" WELL-CONT";
1806 std::ostringstream ss;
1807 const std::streamsize oprec = ss.precision(3);
1808 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1809 ss << std::setw(4) << iteration;
1810 for (
int idx = 0; idx < nm; ++idx) {
1811 ss << std::setw(11) << mass_balance_residual[idx];
1813 for (
int idx = 0; idx < nm; ++idx) {
1814 ss << std::setw(11) << CNV[idx];
1816 for (
int idx = 0; idx < np; ++idx) {
1817 ss << std::setw(11) << well_flux_residual[idx];
1819 ss << std::setw(11) << residualWell;
1821 ss.precision(oprec);
1823 OpmLog::note(ss.str());
1826 for (
int idx = 0; idx < nm; ++idx) {
1827 if (std::isnan(mass_balance_residual[idx])
1828 || std::isnan(CNV[idx])
1829 || (idx < np && std::isnan(well_flux_residual[idx]))) {
1830 const auto msg = std::string(
"NaN residual for phase ") + materialName(idx);
1831 if (terminal_output_) {
1834 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1836 if (mass_balance_residual[idx] > maxResidualAllowed()
1837 || CNV[idx] > maxResidualAllowed()
1838 || (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
1839 const auto msg = std::string(
"Too large residual for phase ") + materialName(idx);
1840 if (terminal_output_) {
1841 OpmLog::problem(msg);
1843 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1846 if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) {
1847 const auto msg = std::string(
"NaN or too large residual for well control equation");
1848 if (terminal_output_) {
1849 OpmLog::problem(msg);
1851 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1861 template <
class Gr
id,
class WellModel,
class Implementation>
1866 const double tol_wells = param_.tolerance_wells_;
1867 const double tol_well_control = param_.tolerance_well_control_;
1869 const int nc = Opm::AutoDiffGrid::numCells(grid_);
1870 const int np = asImpl().numPhases();
1871 const int nm = asImpl().numMaterials();
1873 const V& pv = geo_.poreVolume();
1874 std::vector<double> R_sum(nm);
1875 std::vector<double> B_avg(nm);
1876 std::vector<double> maxCoeff(nm);
1877 std::vector<double> maxNormWell(np);
1878 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
1879 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
1880 Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
1881 for (
int idx = 0; idx < nm; ++idx )
1883 const ADB& tempB = sd_.rq[idx].b;
1884 B.col(idx) = 1./tempB.
value();
1885 R.col(idx) = residual_.material_balance_eq[idx].value();
1886 tempV.col(idx) = R.col(idx).abs()/pv;
1889 convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc);
1891 std::vector<double> well_flux_residual(np);
1892 bool converged_Well =
true;
1894 for (
int idx = 0; idx < np; ++idx )
1896 well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
1897 converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
1900 const double residualWell = detail::infinityNormWell(residual_.well_eq,
1901 linsolver_.parallelInformation());
1902 converged_Well = converged_Well && (residualWell < tol_well_control);
1904 const bool converged = converged_Well;
1907 for (
int idx = 0; idx < np; ++idx) {
1908 if (std::isnan(well_flux_residual[idx])) {
1909 const auto msg = std::string(
"NaN residual for phase ") + materialName(idx);
1910 if (terminal_output_) {
1913 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1915 if (well_flux_residual[idx] > maxResidualAllowed()) {
1916 const auto msg = std::string(
"Too large residual for phase ") + materialName(idx);
1917 if (terminal_output_) {
1918 OpmLog::problem(msg);
1920 OPM_THROW_NOLOG(Opm::NumericalProblem, msg);
1924 if ( terminal_output_ )
1927 if (iteration == 0) {
1930 for (
int idx = 0; idx < np; ++idx) {
1931 msg +=
" W-FLUX(" + materialName(idx).substr(0, 1) +
")";
1933 msg +=
" WELL-CONT";
1936 std::ostringstream ss;
1937 const std::streamsize oprec = ss.precision(3);
1938 const std::ios::fmtflags oflags = ss.setf(std::ios::scientific);
1939 ss << std::setw(4) << iteration;
1940 for (
int idx = 0; idx < np; ++idx) {
1941 ss << std::setw(11) << well_flux_residual[idx];
1943 ss << std::setw(11) << residualWell;
1944 ss.precision(oprec);
1946 OpmLog::note(ss.str());
1955 template <
class Gr
id,
class WellModel,
class Implementation>
1957 BlackoilModelBase<Grid, WellModel, Implementation>::
1958 fluidViscosity(
const int phase,
1963 const std::vector<PhasePresence>& cond)
const
1967 return fluid_.muWat(p, temp, cells_);
1969 return fluid_.muOil(p, temp, rs, cond, cells_);
1971 return fluid_.muGas(p, temp, rv, cond, cells_);
1973 OPM_THROW(std::runtime_error,
"Unknown phase index " << phase);
1981 template <
class Gr
id,
class WellModel,
class Implementation>
1983 BlackoilModelBase<Grid, WellModel, Implementation>::
1984 fluidReciprocFVF(
const int phase,
1989 const std::vector<PhasePresence>& cond)
const
1993 return fluid_.bWat(p, temp, cells_);
1995 return fluid_.bOil(p, temp, rs, cond, cells_);
1997 return fluid_.bGas(p, temp, rv, cond, cells_);
1999 OPM_THROW(std::runtime_error,
"Unknown phase index " << phase);
2007 template <
class Gr
id,
class WellModel,
class Implementation>
2009 BlackoilModelBase<Grid, WellModel, Implementation>::
2010 fluidDensity(
const int phase,
2013 const ADB& rv)
const
2015 const V& rhos = fluid_.surfaceDensity(phase, cells_);
2016 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
2018 if (phase == Oil && active_[Gas]) {
2019 rho += fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * rs * b;
2021 if (phase == Gas && active_[Oil]) {
2022 rho += fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * rv * b;
2031 template <
class Gr
id,
class WellModel,
class Implementation>
2033 BlackoilModelBase<Grid, WellModel, Implementation>::
2034 fluidRsSat(
const V& p,
2036 const std::vector<int>& cells)
const
2045 template <
class Gr
id,
class WellModel,
class Implementation>
2047 BlackoilModelBase<Grid, WellModel, Implementation>::
2048 fluidRsSat(
const ADB& p,
2050 const std::vector<int>& cells)
const
2052 return fluid_.rsSat(p, satOil, cells);
2059 template <
class Gr
id,
class WellModel,
class Implementation>
2061 BlackoilModelBase<Grid, WellModel, Implementation>::
2062 fluidRvSat(
const V& p,
2064 const std::vector<int>& cells)
const
2073 template <
class Gr
id,
class WellModel,
class Implementation>
2075 BlackoilModelBase<Grid, WellModel, Implementation>::
2076 fluidRvSat(
const ADB& p,
2078 const std::vector<int>& cells)
const
2080 return fluid_.rvSat(p, satOil, cells);
2087 template <
class Gr
id,
class WellModel,
class Implementation>
2089 BlackoilModelBase<Grid, WellModel, Implementation>::
2090 poroMult(
const ADB& p)
const
2092 const int n = p.
size();
2093 if (rock_comp_props_ && rock_comp_props_->isActive()) {
2096 #pragma omp parallel for schedule(static)
2097 for (
int i = 0; i < n; ++i) {
2098 pm[i] = rock_comp_props_->poroMult(p.value()[i]);
2099 dpm[i] = rock_comp_props_->poroMultDeriv(p.value()[i]);
2101 ADB::M dpm_diag(dpm.matrix().asDiagonal());
2102 const int num_blocks = p.numBlocks();
2103 std::vector<ADB::M> jacs(num_blocks);
2104 #pragma omp parallel for schedule(dynamic)
2105 for (
int block = 0; block < num_blocks; ++block) {
2118 template <
class Gr
id,
class WellModel,
class Implementation>
2120 BlackoilModelBase<Grid, WellModel, Implementation>::
2121 transMult(
const ADB& p)
const
2123 const int n = p.
size();
2124 if (rock_comp_props_ && rock_comp_props_->isActive()) {
2127 #pragma omp parallel for schedule(static)
2128 for (
int i = 0; i < n; ++i) {
2129 tm[i] = rock_comp_props_->transMult(p.value()[i]);
2130 dtm[i] = rock_comp_props_->transMultDeriv(p.value()[i]);
2132 ADB::M dtm_diag(dtm.matrix().asDiagonal());
2133 const int num_blocks = p.numBlocks();
2134 std::vector<ADB::M> jacs(num_blocks);
2135 #pragma omp parallel for schedule(dynamic)
2136 for (
int block = 0; block < num_blocks; ++block) {
2149 template <
class Gr
id,
class WellModel,
class Implementation>
2151 BlackoilModelBase<Grid, WellModel, Implementation>::
2152 classifyCondition(
const ReservoirState& state)
2154 using namespace Opm::AutoDiffGrid;
2155 const int nc = numCells(grid_);
2156 const int np = state.numPhases();
2158 const PhaseUsage& pu = fluid_.phaseUsage();
2159 const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
2160 if (active_[ Gas ]) {
2162 const V so = s.col(pu.phase_pos[ Oil ]);
2163 const V sg = s.col(pu.phase_pos[ Gas ]);
2165 for (V::Index c = 0, e = sg.size(); c != e; ++c) {
2166 if (so[c] > 0) { phaseCondition_[c].setFreeOil (); }
2167 if (sg[c] > 0) { phaseCondition_[c].setFreeGas (); }
2168 if (active_[ Water ]) { phaseCondition_[c].setFreeWater(); }
2173 assert (active_[ Water ]);
2175 const V so = s.col(pu.phase_pos[ Oil ]);
2178 for (V::Index c = 0, e = so.size(); c != e; ++c) {
2179 phaseCondition_[c].setFreeWater();
2181 if (so[c] > 0) { phaseCondition_[c].setFreeOil(); }
2190 template <
class Gr
id,
class WellModel,
class Implementation>
2195 updatePhaseCondFromPrimalVariable(state);
2203 template <
class Gr
id,
class WellModel,
class Implementation>
2208 const int nc = Opm::AutoDiffGrid::numCells(grid_);
2209 isRs_ = V::Zero(nc);
2210 isRv_ = V::Zero(nc);
2211 isSg_ = V::Zero(nc);
2213 if (! (active_[Gas] && active_[Oil])) {
2215 phaseCondition_.assign(nc, PhasePresence());
2218 for (
int c = 0; c < nc; ++c) {
2219 phaseCondition_[c] = PhasePresence();
2220 phaseCondition_[c].setFreeWater();
2221 switch (state.hydroCarbonState()[c]) {
2222 case HydroCarbonState::GasAndOil:
2223 phaseCondition_[c].setFreeOil();
2224 phaseCondition_[c].setFreeGas();
2227 case HydroCarbonState::OilOnly:
2228 phaseCondition_[c].setFreeOil();
2231 case HydroCarbonState::GasOnly:
2232 phaseCondition_[c].setFreeGas();
2236 OPM_THROW(std::logic_error,
"Unknown primary variable enum value in cell " << c <<
": " << state.hydroCarbonState()[c]);
2245 template <
class Gr
id,
class WellModel,
class Implementation>
2249 const WellState& well_state)
2251 asImpl().
wellModel().computeWellConnectionPressures(state, well_state);
2258 template <
class Gr
id,
class WellModel,
class Implementation>
2259 std::vector<std::vector<double> >
2262 const std::vector<int>& fipnum)
2264 using namespace Opm::AutoDiffGrid;
2265 const int nc = numCells(grid_);
2266 std::vector<ADB> saturation(3,
ADB::null());
2267 const DataBlock s = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, x.numPhases());
2268 const ADB pressure =
ADB::constant(Eigen::Map<const V>(& x.pressure()[0], nc, 1));
2269 const ADB temperature =
ADB::constant(Eigen::Map<const V>(& x.temperature()[0], nc, 1));
2271 saturation[Oil] = active_[Oil] ?
ADB::constant(s.col(Oil)) : ADB::constant(V::Zero(nc));
2272 saturation[Gas] = active_[Gas] ?
ADB::constant(s.col(Gas)) : ADB::constant(V::Zero(nc));
2273 const ADB rs =
ADB::constant(Eigen::Map<const V>(& x.gasoilratio()[0], nc, 1));
2275 const auto canonical_phase_pressures = computePressures(pressure, saturation[Water], saturation[Oil], saturation[Gas]);
2276 const Opm::PhaseUsage& pu = fluid_.phaseUsage();
2277 const std::vector<PhasePresence> cond = phaseCondition();
2279 const ADB pv_mult = poroMult(pressure);
2280 const V& pv = geo_.poreVolume();
2281 const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
2282 for (
int phase = 0; phase < maxnp; ++phase) {
2283 if (active_[ phase ]) {
2284 const int pos = pu.phase_pos[ phase ];
2285 const auto& b = asImpl().fluidReciprocFVF(phase, canonical_phase_pressures[phase], temperature, rs, rv, cond);
2286 sd_.fip[phase] = ((pv_mult * b * saturation[pos] * pv).value());
2290 if (active_[ Oil ] && active_[ Gas ]) {
2292 sd_.fip[SimulatorData::FIP_DISSOLVED_GAS] = rs.value() * sd_.fip[SimulatorData::FIP_LIQUID];
2293 sd_.fip[SimulatorData::FIP_VAPORIZED_OIL] = rv.value() * sd_.fip[SimulatorData::FIP_VAPOUR];
2297 int dims = *std::max_element(fipnum.begin(), fipnum.end());
2298 std::vector<std::vector<double> > values(dims);
2299 for (
int i=0; i < dims; ++i) {
2300 values[i].resize(7, 0.0);
2303 const V hydrocarbon = saturation[Oil].value() + saturation[Gas].value();
2307 if ( !isParallel() )
2310 for (
int phase = 0; phase < maxnp; ++phase) {
2311 if (active_[ phase ]) {
2312 for (
int c = 0; c < nc; ++c) {
2313 const int region = fipnum[c] - 1;
2315 values[region][phase] += sd_.fip[phase][c];
2322 if (active_[ Oil ] && active_[ Gas ]) {
2323 for (
int c = 0; c < nc; ++c) {
2324 const int region = fipnum[c] - 1;
2326 values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
2327 values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
2333 hcpv = V::Zero(dims);
2334 pres = V::Zero(dims);
2336 for (
int c = 0; c < nc; ++c) {
2337 const int region = fipnum[c] - 1;
2339 hcpv[region] += pv[c] * hydrocarbon[c];
2340 pres[region] += pv[c] * pressure.value()[c];
2344 sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
2345 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
2347 for (
int c = 0; c < nc; ++c) {
2348 const int region = fipnum[c] - 1;
2350 sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
2354 if (hcpv[region] != 0) {
2355 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
2357 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
2360 values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
2361 values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
2369 const auto & pinfo =
2370 boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation());
2371 const auto& mask = pinfo.getOwnerMask();
2372 auto comm = pinfo.communicator();
2374 dims = comm.max(dims);
2375 values.resize(dims);
2376 for (
int i=0; i < dims; ++i) {
2377 values[i].resize(7);
2378 std::fill(values[i].begin(), values[i].end(), 0.0);
2382 for (
int phase = 0; phase < maxnp; ++phase) {
2383 for (
int c = 0; c < nc; ++c) {
2384 const int region = fipnum[c] - 1;
2385 if (region != -1 && mask[c]) {
2386 values[region][phase] += sd_.fip[phase][c];
2392 if (active_[ Oil ] && active_[ Gas ]) {
2393 for (
int c = 0; c < nc; ++c) {
2394 const int region = fipnum[c] - 1;
2395 if (region != -1 && mask[c]) {
2396 values[region][SimulatorData::FIP_DISSOLVED_GAS] += sd_.fip[SimulatorData::FIP_DISSOLVED_GAS][c];
2397 values[region][SimulatorData::FIP_VAPORIZED_OIL] += sd_.fip[SimulatorData::FIP_VAPORIZED_OIL][c];
2402 hcpv = V::Zero(dims);
2403 pres = V::Zero(dims);
2405 for (
int c = 0; c < nc; ++c) {
2406 const int region = fipnum[c] - 1;
2407 if (region != -1 && mask[c]) {
2408 hcpv[region] += pv[c] * hydrocarbon[c];
2409 pres[region] += pv[c] * pressure.value()[c];
2413 comm.sum(hcpv.data(), hcpv.size());
2414 comm.sum(pres.data(), pres.size());
2416 sd_.fip[SimulatorData::FIP_PV] = V::Zero(nc);
2417 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE] = V::Zero(nc);
2419 for (
int c = 0; c < nc; ++c) {
2420 const int region = fipnum[c] - 1;
2421 if (region != -1 && mask[c]) {
2422 sd_.fip[SimulatorData::FIP_PV][c] = pv[c];
2424 if (hcpv[region] != 0) {
2425 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pv[c] * pressure.value()[c] * hydrocarbon[c] / hcpv[region];
2427 sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c] = pres[region] / pv[c];
2430 values[region][SimulatorData::FIP_PV] += sd_.fip[SimulatorData::FIP_PV][c];
2431 values[region][SimulatorData::FIP_WEIGHTED_PRESSURE] += sd_.fip[SimulatorData::FIP_WEIGHTED_PRESSURE][c];
2438 for(
int reg=0; reg < dims; ++reg)
2440 comm.sum(values[reg].data(), values[reg].size());
2444 OPM_THROW(std::logic_error,
"HAVE_MPI should be defined if we are running in parallel");
2455 template <
class Gr
id,
class WellModel,
class Implementation>
2459 const WellState& well_state,
2460 std::vector<double>& well_voidage_rates,
2461 std::vector<double>& voidage_conversion_coeffs)
2468 const int nw = well_state.numWells();
2469 const int np = numPhases();
2471 const Wells* wells = asImpl().wellModel().wellsPointer();
2474 well_voidage_rates.resize(nw, 0);
2476 voidage_conversion_coeffs.resize(nw * np, 1.0);
2478 int global_number_wells = nw;
2481 if ( linsolver_.parallelInformation().type() ==
typeid(ParallelISTLInformation) )
2484 boost::any_cast<
const ParallelISTLInformation&>(linsolver_.parallelInformation());
2485 global_number_wells = info.communicator().sum(global_number_wells);
2486 if ( global_number_wells )
2488 rate_converter_.defineState(reservoir_state, boost::any_cast<const ParallelISTLInformation&>(linsolver_.parallelInformation()));
2494 if ( global_number_wells )
2496 rate_converter_.defineState(reservoir_state);
2500 std::vector<double> well_rates(np, 0.0);
2501 std::vector<double> convert_coeff(np, 1.0);
2504 if ( !well_voidage_rates.empty() ) {
2505 for (
int w = 0; w < nw; ++w) {
2506 const bool is_producer = wells->type[w] == PRODUCER;
2510 std::transform(well_state.wellRates().begin() + np * w,
2511 well_state.wellRates().begin() + np * (w + 1),
2512 well_rates.begin(), std::negate<double>());
2515 const int fipreg = 0;
2516 const int well_cell_top = wells->well_cells[wells->well_connpos[w]];
2517 const int pvtreg = fluid_.cellPvtRegionIndex()[well_cell_top];
2519 rate_converter_.calcCoeff(fipreg, pvtreg, convert_coeff);
2520 well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(),
2521 convert_coeff.begin(), 0.0);
2525 std::copy(well_state.wellRates().begin() + np * w,
2526 well_state.wellRates().begin() + np * (w + 1),
2527 well_rates.begin());
2529 const int fipreg = 0;
2530 const int well_cell_top = wells->well_cells[wells->well_connpos[w]];
2531 const int pvtreg = fluid_.cellPvtRegionIndex()[well_cell_top];
2532 rate_converter_.calcCoeff(fipreg, pvtreg, convert_coeff);
2533 std::copy(convert_coeff.begin(), convert_coeff.end(),
2534 voidage_conversion_coeffs.begin() + np * w);
2544 template <
class Gr
id,
class WellModel,
class Implementation>
2548 WellState& well_state)
2550 if (asImpl().wellModel().wellCollection()->havingVREPGroups()) {
2551 std::vector<double> well_voidage_rates;
2552 std::vector<double> voidage_conversion_coeffs;
2553 computeWellVoidageRates(reservoir_state, well_state, well_voidage_rates, voidage_conversion_coeffs);
2554 asImpl().wellModel().wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
2557 for (
const WellNode* well_node : asImpl().wellModel().wellCollection()->getLeafNodes()) {
2558 if (well_node->isInjector() && !well_node->individualControl()) {
2559 const int well_index = well_node->selfIndex();
2560 well_state.currentControls()[well_index] = well_node->groupControlIndex();
2570 template <
class Gr
id,
class WellModel,
class Implementation>
2574 WellState& well_state)
2576 if (asImpl().wellModel().wellCollection()->requireWellPotentials()) {
2577 SolutionState state = asImpl().variableState(reservoir_state, well_state);
2578 asImpl().makeConstantState(state);
2579 asImpl().wellModel().computeWellConnectionPressures(state, well_state);
2581 const int np = numPhases();
2585 const ADB& press = state.pressure;
2586 const ADB& temp = state.temperature;
2587 const ADB& rs = state.rs;
2588 const ADB& rv = state.rv;
2590 const std::vector<PhasePresence> cond = phaseCondition();
2592 const ADB pv_mult = poroMult(press);
2593 const ADB tr_mult = transMult(press);
2594 const std::vector<ADB> kr = asImpl().computeRelPerm(state);
2596 std::vector<ADB> mob_perfcells(np,
ADB::null());
2597 std::vector<ADB> b_perfcells(np,
ADB::null());
2599 for (
int phase = 0; phase < np; ++phase) {
2600 if (active_[phase]) {
2601 const std::vector<int>& well_cells = asImpl().wellModel().wellOps().well_cells;
2602 const ADB mu = asImpl().fluidViscosity(canph_[phase], state.canonical_phase_pressures[canph_[phase]],
2603 temp, rs, rv, cond);
2604 mob[phase] = tr_mult * kr[canph_[phase]] / mu;
2605 mob_perfcells[phase] =
subset(mob[phase], well_cells);
2607 b[phase] = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
2608 b_perfcells[phase] =
subset(b[phase], well_cells);
2613 std::vector<double> well_potentials;
2614 asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, well_state, state, well_potentials);
2615 asImpl().wellModel().wellCollection()->setGuideRatesWithPotentials(asImpl().wellModel().wellsPointer(),
2616 fluid_.phaseUsage(), well_potentials);
2619 applyVREPGroupControl(reservoir_state, well_state);
2621 if (asImpl().wellModel().wellCollection()->groupControlApplied()) {
2622 asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
2624 asImpl().wellModel().wellCollection()->applyGroupControls();
2628 const int nw = wells().number_of_wells;
2629 for (
int w = 0; w < nw; ++w) {
2630 const WellNode& well_node = asImpl().wellModel().wellCollection()->findWellNode(wells().name[w]);
2631 if (!well_node.individualControl()) {
2632 well_state.currentControls()[w] = well_node.groupControlIndex();
2640 #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED
static AutoDiffBlock null()
Construct an empty AutoDiffBlock.
Definition: AutoDiffBlock.hpp:104
int sizeNonLinear() const
The size (number of unknowns) of the nonlinear system of equations.
Definition: BlackoilModelBase_impl.hpp:347
double convergenceReduction(const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > &B, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > &tempV, const Eigen::Array< double, Eigen::Dynamic, Eigen::Dynamic > &R, std::vector< double > &R_sum, std::vector< double > &maxCoeff, std::vector< double > &B_avg, std::vector< double > &maxNormWell, int nc) const
Compute the reduction within the convergence check.
Definition: BlackoilModelBase_impl.hpp:1634
void computeWellVoidageRates(const ReservoirState &reservoir_state, const WellState &well_state, std::vector< double > &well_voidage_rates, std::vector< double > &voidage_conversion_coeffs)
Function to compute the resevoir voidage for the production wells.
Definition: BlackoilModelBase_impl.hpp:2458
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:438
bool getConvergence(const SimulatorTimerInterface &timer, const int iteration)
Compute convergence based on total mass balance (tol_mb) and maximum residual mass balance (tol_cnv)...
Definition: BlackoilModelBase_impl.hpp:1721
SimulatorReport assemble(const ReservoirState &reservoir_state, WellState &well_state, const bool initial_assembly)
Assemble the residual and Jacobian of the nonlinear system.
Definition: BlackoilModelBase_impl.hpp:758
AutoDiffMatrix is a wrapper class that optimizes matrix operations.
Definition: AutoDiffMatrix.hpp:43
AutoDiffMatrix M
Underlying type for jacobians.
Definition: AutoDiffBlock.hpp:101
bool terminalOutputEnabled() const
Return true if output to cout is wanted.
Definition: BlackoilModelBase_impl.hpp:371
double relativeChange(const SimulationDataContainer &previous, const SimulationDataContainer ¤t) const
compute the relative change between to simulation states
Definition: BlackoilModelBase_impl.hpp:1593
int numPhases() const
Definition: BlackoilPropsAdFromDeck.cpp:225
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:415
A model implementation for three-phase black oil.
Definition: BlackoilModelBase.hpp:76
void updateState(const V &dx, ReservoirState &reservoir_state, WellState &well_state)
Apply an update to the primary variables, chopped if appropriate.
Definition: BlackoilModelBase_impl.hpp:1148
void prepareStep(const SimulatorTimerInterface &timer, const ReservoirState &reservoir_state, const WellState &well_state)
Called once before each time step.
Definition: BlackoilModelBase_impl.hpp:220
void afterStep(const SimulatorTimerInterface &timer, ReservoirState &reservoir_state, WellState &well_state)
Called once after each time step.
Definition: BlackoilModelBase_impl.hpp:333
virtual double currentStepLength() const =0
Current step length.
std::vector< std::vector< double > > computeFluidInPlace(const ReservoirState &x, const std::vector< int > &fipnum)
Compute fluid in place.
Definition: BlackoilModelBase_impl.hpp:2261
void updatePrimalVariableFromState(const ReservoirState &state)
update the primal variable for Sg, Rv or Rs.
Definition: BlackoilModelBase_impl.hpp:2193
int numPhases() const
The number of active fluid phases in the model.
Definition: BlackoilModelBase_impl.hpp:383
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:319
WellModel & wellModel()
return the WellModel object
Definition: BlackoilModelBase.hpp:268
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
BlackoilModelBase(const ModelParameters ¶m, const Grid &grid, const BlackoilPropsAdFromDeck &fluid, const DerivedGeology &geo, const RockCompressibility *rock_comp_props, const WellModel &well_model, const NewtonIterationBlackoilInterface &linsolver, std::shared_ptr< const EclipseState > eclState, const bool has_disgas, const bool has_vapoil, const bool terminal_output)
Construct the model.
Definition: BlackoilModelBase_impl.hpp:101
void setThresholdPressures(const std::vector< double > &threshold_pressures_by_face)
Set threshold pressures that prevent or reduce flow.
Definition: BlackoilModelBase_impl.hpp:419
std::vector< double > computeResidualNorms() const
Compute the residual norms of the mass balance for each phase, the well flux, and the well equation...
Definition: BlackoilModelBase_impl.hpp:1552
void updateEquationsScaling()
Update the scaling factors for mass balance equations.
Definition: BlackoilModelBase_impl.hpp:920
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
int numMaterials() const
The number of active materials in the model.
Definition: BlackoilModelBase_impl.hpp:394
int linearIterationsLastSolve() const
Number of linear iterations used in last call to solveJacobianSystem().
Definition: BlackoilModelBase_impl.hpp:359
Definition: AutoDiffHelpers.hpp:623
bool isParallel() const
Return true if this is a parallel run.
Definition: BlackoilModelBase_impl.hpp:198
AutoDiffBlock< double > vertcatCollapseJacs(const std::vector< AutoDiffBlock< double > > &x)
Returns the vertical concatenation [ x[0]; x[1]; ...; x[n-1] ] of the inputs.
Definition: AutoDiffHelpers.hpp:530
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
const std::string & materialName(int material_index) const
The name of an active material in the model.
Definition: BlackoilModelBase_impl.hpp:406
static AutoDiffBlock function(V &&val, std::vector< M > &&jac)
Create an AutoDiffBlock by directly specifying values and jacobians.
Definition: AutoDiffBlock.hpp:184
void setupGroupControl(const ReservoirState &reservoir_state, WellState &well_state)
Set up the group control related at the beginning of each time step.
Definition: BlackoilModelBase_impl.hpp:2573
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
V solveJacobianSystem() const
Solve the Jacobian system Jx = r where J is the Jacobian and r is the residual.
Definition: BlackoilModelBase_impl.hpp:1140
void updatePhaseCondFromPrimalVariable(const ReservoirState &state)
Update the phaseCondition_ member based on the primalVariable_ member.
Definition: BlackoilModelBase_impl.hpp:2206
Eigen::Array< double, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99
void fastSparseProduct(const AutoDiffMatrix &lhs, const AutoDiffMatrix &rhs, AutoDiffMatrix &res)
Utility function to lessen code changes required elsewhere.
Definition: AutoDiffMatrix.hpp:708
static std::vector< AutoDiffBlock > variables(const std::vector< V > &initial_values)
Construct a set of primary variables, each initialized to a given vector.
Definition: AutoDiffBlock.hpp:203
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:719
SimulatorReport nonlinearIteration(const int iteration, const SimulatorTimerInterface &timer, NonlinearSolverType &nonlinear_solver, ReservoirState &reservoir_state, WellState &well_state)
Called once per nonlinear iteration.
Definition: BlackoilModelBase_impl.hpp:240