diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 9ec5d7cad..6318a3772 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -160,9 +160,6 @@ namespace Opm { void computeWellConnectionPressures(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw); - void - addOldWellEq(const SolutionState& state); - void addWellControlEq(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw, diff --git a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp index ded3d6672..f64debb52 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp @@ -1127,128 +1127,6 @@ namespace { - template - void FullyImplicitBlackoilSolver::addOldWellEq(const SolutionState& state) - { - // -------- Well equation, and well contributions to the mass balance equations -------- - - // Contribution to mass balance will have to wait. - - const int nc = numCells(grid_); - const int np = wells_.number_of_phases; - const int nw = wells_.number_of_wells; - const int nperf = wells_.well_connpos[nw]; - - const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); - const V transw = Eigen::Map(wells_.WI, nperf); - - const ADB& bhp = state.bhp; - - const DataBlock well_s = wops_.w2p * Eigen::Map(wells_.comp_frac, nw, np).matrix(); - - // Extract variables for perforation cell pressures - // and corresponding perforation well pressures. - const ADB p_perfcell = subset(state.pressure, well_cells); - // Finally construct well perforation pressures and well flows. - - // Compute well pressure differentials. - // Construct pressure difference vector for wells. - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - const int dim = dimensions(grid_); - const double* g = geo_.gravity(); - if (g) { - // Guard against gravity in anything but last dimension. - for (int dd = 0; dd < dim - 1; ++dd) { - assert(g[dd] == 0.0); - } - } - - // make a copy of the phaseConditions - std::vector cond = phaseCondition_; - - ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern()); - for (int phase = 0; phase < 3; ++phase) { - if (active_[phase]) { - const int pos = pu.phase_pos[phase]; - const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, state.rv,cond, cells_); - cell_rho_total += state.saturation[pos] * cell_rho; - } - } - ADB inj_rho_total = ADB::constant(V::Zero(nperf), state.pressure.blockPattern()); - assert(np == wells_.number_of_phases); - const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); - for (int phase = 0; phase < 3; ++phase) { - if (active_[phase]) { - const int pos = pu.phase_pos[phase]; - const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, state.rv,cond, cells_); - const V fraction = compi.col(pos); - inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells); - } - } - const V rho_perf_cell = subset(cell_rho_total, well_cells).value(); - const V rho_perf_well = inj_rho_total.value(); - V prodperfs = V::Constant(nperf, -1.0); - for (int w = 0; w < nw; ++w) { - if (wells_.type[w] == PRODUCER) { - std::fill(prodperfs.data() + wells_.well_connpos[w], - prodperfs.data() + wells_.well_connpos[w+1], 1.0); - } - } - const Selector producer(prodperfs); - const V rho_perf = producer.select(rho_perf_cell, rho_perf_well); - const V well_perf_dp = computePerfPress(grid_, wells_, rho_perf, g ? g[dim-1] : 0.0); - - const ADB p_perfwell = wops_.w2p * bhp + well_perf_dp; - const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); - // DUMP(nkgradp_well); - const Selector cell_to_well_selector(nkgradp_well.value()); - ADB well_rates_all = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern()); - ADB perf_total_mob = subset(rq_[0].mob, well_cells); - for (int phase = 1; phase < np; ++phase) { - perf_total_mob += subset(rq_[phase].mob, well_cells); - } - std::vector well_contribs(np, ADB::null()); - std::vector well_perf_rates(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const ADB& cell_b = rq_[phase].b; - const ADB perf_b = subset(cell_b, well_cells); - const ADB& cell_mob = rq_[phase].mob; - const V well_fraction = compi.col(phase); - // Using total mobilities for all phases for injection. - const ADB perf_mob_injector = (wops_.w2p * well_fraction.matrix()).array() * perf_total_mob; - const ADB perf_mob = producer.select(subset(cell_mob, well_cells), - perf_mob_injector); - const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations. - well_perf_rates[phase] = (perf_flux*perf_b); - const ADB well_rates = wops_.p2w * well_perf_rates[phase]; - well_rates_all += superset(well_rates, Span(nw, 1, phase*nw), nw*np); - - // const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc); - well_contribs[phase] = superset(perf_flux*perf_b, well_cells, nc); - // DUMP(well_contribs[phase]); - residual_.material_balance_eq[phase] += well_contribs[phase]; - } - if (active_[Gas] && active_[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const ADB rs_perf = subset(state.rs, well_cells); - const ADB rv_perf = subset(state.rv, well_cells); - well_rates_all += superset(wops_.p2w * (well_perf_rates[oilpos]*rs_perf), Span(nw, 1, gaspos*nw), nw*np); - well_rates_all += superset(wops_.p2w * (well_perf_rates[gaspos]*rv_perf), Span(nw, 1, oilpos*nw), nw*np); - // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs); - residual_.material_balance_eq[gaspos] += well_contribs[oilpos]*state.rs; - residual_.material_balance_eq[oilpos] += well_contribs[gaspos]*state.rv; - } - - // Set the well flux equation - residual_.well_flux_eq = state.qs + well_rates_all; - // DUMP(residual_.well_flux_eq); - } - - - - - template V FullyImplicitBlackoilSolver::solveJacobianSystem() const {