diff --git a/opm/autodiff/BlackoilReorderingTransportModel.hpp b/opm/autodiff/BlackoilReorderingTransportModel.hpp index c2d633b6c..1f1780995 100644 --- a/opm/autodiff/BlackoilReorderingTransportModel.hpp +++ b/opm/autodiff/BlackoilReorderingTransportModel.hpp @@ -276,7 +276,7 @@ namespace Opm { for (int cell = 0; cell < num_cells; ++cell) { computeCellState(cell, state0_, cstate0_[cell]); } - cstate_.resize(num_cells); + cstate_ = cstate0_; } @@ -377,6 +377,24 @@ namespace Opm { { return s[phaseIdx]; } + + template + CellState flatten() const + { + return CellState{ + { s[0].value, s[1].value, s[2].value }, + rs.value, + rv.value, + { p[0].value, p[1].value, p[2].value }, + { kr[0].value, kr[1].value, kr[2].value }, + { pc[0].value, pc[1].value, pc[2].value }, + temperature.value, + { mu[0].value, mu[1].value, mu[2].value }, + { b[0].value, b[1].value, b[2].value }, + { lambda[0].value, lambda[1].value, lambda[2].value }, + { rho[0].value, rho[1].value, rho[2].value } + }; + } }; @@ -395,7 +413,7 @@ namespace Opm { State state_; std::vector> cstate0_; - std::vector> cstate_; + std::vector> cstate_; V total_flux_; V total_wellperf_flux_; @@ -566,7 +584,7 @@ namespace Opm { while (!getConvergence(res)) { Vec2 dx; jac.solve(dx, res); - updateState(dx); + updateState(cell, dx); assembleSingleCell(cell, res, jac); } } @@ -621,13 +639,15 @@ namespace Opm { { assert(numPhases() == 3); // I apologize for this to my future self, that will have to fix it. - computeCellState(cell, state_, cstate_[cell]); + CellState st; + computeCellState(cell, state_, st); + cstate_[cell] = st.template flatten(); // Accumulation terms. - const Eval ao0 = oilAccumulation(cstate0_[cell]); - const Eval ao = oilAccumulation(cstate_[cell]); - const Eval ag0 = gasAccumulation(cstate0_[cell]); - const Eval ag = gasAccumulation(cstate_[cell]); + const double ao0 = oilAccumulation(cstate0_[cell]); + const Eval ao = oilAccumulation(st); + const double ag0 = gasAccumulation(cstate0_[cell]); + const Eval ag = gasAccumulation(st); // Flux terms. Eval div_oilflux = Eval::createConstant(0.0); @@ -650,8 +670,8 @@ namespace Opm { Eval dh[3]; Eval dh_sat[3]; for (int phase : { Water, Oil, Gas }) { - const Eval gradp = cstate_[other].p[phase].value - cstate_[cell].p[phase]; - const Eval rhoavg = 0.5 * (cstate_[cell].rho[phase] + cstate_[other].rho[phase].value); + const Eval gradp = cstate_[other].p[phase] - st.p[phase]; + const Eval rhoavg = 0.5 * (st.rho[phase] + cstate_[other].rho[phase]); dh[phase] = gradp - rhoavg * gdz_[conn.index]; dh_sat[phase] = rhoavg * gdz_[conn.index]; if (Base::use_threshold_pressure_) { @@ -659,22 +679,22 @@ namespace Opm { } } const double tran = trans_all_[conn.index]; - const auto& m1 = cstate_[cell].lambda; + const auto& m1 = st.lambda; const auto& m2 = cstate_[other].lambda; const auto upw = connectionMultiPhaseUpwind({{ dh_sat[Water].value, dh_sat[Oil].value, dh_sat[Gas].value }}, {{ m1[Water].value, m1[Oil].value, m1[Gas].value }}, - {{ m2[Water].value, m2[Oil].value, m2[Gas].value }}, + {{ m2[Water], m2[Oil], m2[Gas] }}, tran, vt); Eval b[3]; Eval mob[3]; Eval tot_mob = Eval::createConstant(0.0); for (int phase : { Water, Oil, Gas }) { - b[phase] = upw[phase] > 0.0 ? cstate_[cell].b[phase] : cstate_[other].b[phase].value; - mob[phase] = upw[phase] > 0.0 ? m1[phase] : m2[phase].value; + b[phase] = upw[phase] > 0.0 ? st.b[phase] : cstate_[other].b[phase]; + mob[phase] = upw[phase] > 0.0 ? m1[phase] : m2[phase]; tot_mob += mob[phase]; } - Eval rs = upw[Oil] > 0.0 ? cstate_[cell].rs : cstate_[other].rs.value; - Eval rv = upw[Gas] > 0.0 ? cstate_[cell].rv : cstate_[other].rv.value; + Eval rs = upw[Oil] > 0.0 ? st.rs : cstate_[other].rs; + Eval rv = upw[Gas] > 0.0 ? st.rv : cstate_[other].rv; Eval flux[3]; for (int phase : { Oil, Gas }) { @@ -714,9 +734,67 @@ namespace Opm { - void updateState(const Vec2& dx) + void updateState(const int cell, + const Vec2& dx) { - // TODO: update... + // Get saturation updates. + const double dsw = dx[0]; + double dso = -dsw; + double dsg = 0.0; + auto& hcstate = state_.reservoir_state.hydroCarbonState()[cell]; + if (hcstate == HydroCarbonState::GasAndOil) { + dsg = dx[1]; + dso -= dsg; + } + + // Handle too large saturation changes. + const double maxval = std::max(std::fabs(dsw), std::max(std::fabs(dso), std::fabs(dsg))); + const double sfactor = std::min(1.0, Base::dsMax() / maxval); + double* s = state_.reservoir_state.saturation().data() + 3*cell; + s[Water] -= sfactor*dsw; + s[Gas] -= sfactor*dsg; + s[Oil] = 1.0 - s[Water] - s[Oil]; + + // Handle < 0 saturations. + for (int phase : { Gas, Oil, Water }) { // TODO: check if ordering here is significant + if (s[phase] < 0.0) { + for (int other_phase : { Water, Oil, Gas }) { + if (phase != other_phase) { + s[other_phase] /= (1.0 - s[phase]); + } + } + s[phase] = 0.0; + } + } + + // Update rs. + double& rs = state_.reservoir_state.gasoilratio()[cell]; + if (hcstate == HydroCarbonState::OilOnly) { + const double rs_old = rs; + const double max_allowed_change = std::fabs(rs_old) * Base::drMaxRel(); + const double drs = dx[1]; + const double factor = std::min(1.0, max_allowed_change / std::fabs(drs)); + rs -= factor*drs; + } + + // Update rv. + double& rv = state_.reservoir_state.rv()[cell]; + if (hcstate == HydroCarbonState::GasOnly) { + const double rv_old = rv; + const double max_allowed_change = std::fabs(rv_old) * Base::drMaxRel(); + const double drv = dx[1]; + const double factor = std::min(1.0, max_allowed_change / std::fabs(drv)); + rv -= factor*drv; + } + + const double epsilon = std::sqrt(std::numeric_limits::epsilon()); + const bool water_only = s[Water] > (1 - epsilon); + + // hcstate = HydroCarbonState::GasAndOil; + + // rssat0 = ; + // rssat = + } };