From adf291a30c42c7610cf3bdd478c905114d24e09d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 08:17:13 +0200 Subject: [PATCH 1/5] Work in progress on well flux equations. --- .../.#FullyImplicitBlackoilSolver.hpp | 1 + opm/autodiff/FullyImplicitBlackoilSolver.cpp | 24 +++++++++++++++---- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 1 + 3 files changed, 22 insertions(+), 4 deletions(-) create mode 120000 opm/autodiff/.#FullyImplicitBlackoilSolver.hpp diff --git a/opm/autodiff/.#FullyImplicitBlackoilSolver.hpp b/opm/autodiff/.#FullyImplicitBlackoilSolver.hpp new file mode 120000 index 000000000..a680f0a09 --- /dev/null +++ b/opm/autodiff/.#FullyImplicitBlackoilSolver.hpp @@ -0,0 +1 @@ +atgeirr@arwen.lan.6100 \ No newline at end of file diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 17f8d0c63..ee65b28bf 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -198,6 +198,7 @@ namespace Opm { , rq_ (fluid.numPhases()) , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), ADB::null(), + std::vector(fluid.numPhases(), ADB::null()), ADB::null() } ) { } @@ -643,7 +644,7 @@ namespace Opm { const ADB p_perfwell = wops_.w2p * bhp + well_perf_dp; const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); - DUMP(nkgradp_well); + // DUMP(nkgradp_well); const Selector cell_to_well_selector(nkgradp_well.value()); ADB qs = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern()); // We can safely use a dummy rs here (for well calculations) @@ -655,6 +656,7 @@ namespace Opm { 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); @@ -665,8 +667,9 @@ namespace Opm { 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. - const ADB well_rates = wops_.p2w * (perf_flux*perf_b); - qs = qs + superset(well_rates, Span(nw, 1, phase*nw), nw*np); + well_perf_rates[phase] = (perf_flux*perf_b); + const ADB well_rates = wops_.p2w * well_perf_rates[phase]; + qs += 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); @@ -676,6 +679,8 @@ namespace Opm { 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); + qs += superset(well_perf_rates[oilpos]*rs_perf, Span(nw, 1, gaspos*nw), nw*np); // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.Rs); residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.Rs; } @@ -703,6 +708,7 @@ namespace Opm { // Choose bhp residual for positive bhp targets. Selector bhp_selector(bhp_targets); residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual); + DUMP(residual_.well_eq); } @@ -779,7 +785,7 @@ namespace Opm { const V p_old = Eigen::Map(&state.pressure()[0], nc, 1); const V absdpmax = dpmaxrel*p_old.abs(); const V dpsign = dp/dp.abs(); - const V dp_limited = dpsign * dp.abs().max(absdpmax); + const V dp_limited = dpsign * dp.abs().min(absdpmax); const V p = (p_old - dp_limited).max(zero); std::copy(&p[0], &p[0] + nc, state.pressure().begin()); @@ -978,6 +984,16 @@ namespace Opm { { r = std::max(r, (*b).value().matrix().norm()); } + if (active_[Oil] && active_[Gas]) { + r = std::max(r, residual_.rs_or_sg_eq.value().matrix().norm()); + } + for (std::vector::const_iterator + b = residual_.well_flux_eq.begin(), + e = residual_.well_flux_eq.end(); + b != e; ++b) + { + r = std::max(r, (*b).value().matrix().norm()); + } r = std::max(r, residual_.well_eq.value().matrix().norm()); return r; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 040173c11..3a8e571c9 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -116,6 +116,7 @@ namespace Opm { struct { std::vector mass_balance; ADB rs_or_sg_eq; // Only used if both gas and oil present + std::vector well_flux_eq; ADB well_eq; } residual_; From 74a5e10f7b876e1173a1a640d394f1667523a2e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 08:19:21 +0200 Subject: [PATCH 2/5] Renames variables dtpv->pvdt and Rs->rs. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 44 ++++++++++---------- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 2 +- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index ee65b28bf..747d0f562 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -213,7 +213,7 @@ namespace Opm { BlackoilState& x , WellState& xw) { - const V dtpv = geo_.poreVolume() / dt; + const V pvdt = geo_.poreVolume() / dt; { const SolutionState state = constantState(x, xw); @@ -224,7 +224,7 @@ namespace Opm { const double rtol = 5.0e-8; const int maxit = 15; - assemble(dtpv, x, xw); + assemble(pvdt, x, xw); const double r0 = residualNorm(); int it = 0; @@ -237,7 +237,7 @@ namespace Opm { updateState(dx, x, xw); - assemble(dtpv, x, xw); + assemble(pvdt, x, xw); const double r = residualNorm(); @@ -274,7 +274,7 @@ namespace Opm { FullyImplicitBlackoilSolver::SolutionState::SolutionState(const int np) : pressure ( ADB::null()) , saturation(np, ADB::null()) - , Rs ( ADB::null()) + , rs ( ADB::null()) , bhp ( ADB::null()) { } @@ -370,13 +370,13 @@ namespace Opm { } } - // Gas-oil ratio (Rs). + // Gas-oil ratio (rs). if (active_[ Oil ] && active_[ Gas ]) { const V rs = Eigen::Map(& x.gasoilratio()[0], x.gasoilratio().size()); - state.Rs = ADB::constant(rs, bpat); + state.rs = ADB::constant(rs, bpat); } else { const V Rs = V::Zero(nc, 1); - state.Rs = ADB::constant(Rs, bpat); + state.rs = ADB::constant(Rs, bpat); } // Well bottom-hole pressure. @@ -465,9 +465,9 @@ namespace Opm { // Rs. if (active_[ Oil ] && active_[ Gas ]) { - state.Rs = vars[ nextvar++ ]; + state.rs = vars[ nextvar++ ]; } else { - state.Rs = ADB::constant(V::Zero(nc), bpat); + state.rs = ADB::constant(V::Zero(nc), bpat); } // Bhp. @@ -490,7 +490,7 @@ namespace Opm { const ADB& press = state.pressure; const std::vector& sat = state.saturation; - const ADB& rs = state.Rs; + const ADB& rs = state.rs; const int maxnp = Opm::BlackoilPhases::MaxNumPhases; for (int phase = 0; phase < maxnp; ++phase) { @@ -508,7 +508,7 @@ namespace Opm { const int po = pu.phase_pos[ Oil ]; const int pg = pu.phase_pos[ Gas ]; - rq_[pg].accum[aix] += state.Rs * rq_[po].accum[aix]; + rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix]; // DUMP(rq_[pg].accum[aix]); } } @@ -519,7 +519,7 @@ namespace Opm { void FullyImplicitBlackoilSolver:: - assemble(const V& dtpv, + assemble(const V& pvdt, const BlackoilState& x , const WellState& xw ) { @@ -548,7 +548,7 @@ namespace Opm { // std::cout << rq_[phase].mflux; residual_.mass_balance[ phase ] = - dtpv*(rq_[phase].accum[1] - rq_[phase].accum[0]) + pvdt*(rq_[phase].accum[1] - rq_[phase].accum[0]) + ops_.div*rq_[phase].mflux; // DUMP(residual_.mass_balance[phase]); @@ -563,7 +563,7 @@ namespace Opm { const int po = fluid_.phaseUsage().phase_pos[ Oil ]; const UpwindSelector upwind(grid_, ops_, rq_[po].head.value()); - const ADB rs_face = upwind.select(state.Rs); + const ADB rs_face = upwind.select(state.rs); residual_.mass_balance[ Gas ] += ops_.div * (rs_face * rq_[po].mflux); // DUMP(residual_.mass_balance[ Gas ]); @@ -572,7 +572,7 @@ namespace Opm { const int pg = fluid_.phaseUsage().phase_pos[ Gas ]; const ADB sg_eq = state.saturation[pg]; const ADB rs_max = fluidRsMax(state.pressure, cells_); - const ADB rs_eq = state.Rs - rs_max; + const ADB rs_eq = state.rs - rs_max; Selector use_rs_eq(rs_eq.value()); residual_.rs_or_sg_eq = use_rs_eq.select(rs_eq, sg_eq); // DUMP(residual_.rs_or_sg_eq); @@ -614,7 +614,7 @@ namespace Opm { 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, cells_); + const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, cells_); cell_rho_total += state.saturation[pos] * cell_rho; } } @@ -624,7 +624,7 @@ namespace Opm { 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, cells_); + const ADB cell_rho = fluidDensity(phase, state.pressure, state.rs, cells_); const V fraction = compi.col(pos); inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells); } @@ -679,10 +679,10 @@ namespace Opm { 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 rs_perf = subset(state.rs, well_cells); qs += superset(well_perf_rates[oilpos]*rs_perf, Span(nw, 1, gaspos*nw), nw*np); - // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.Rs); - residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.Rs; + // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs); + residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs; } // Handling BHP and SURFACE_RATE wells. V bhp_targets(nw); @@ -950,11 +950,11 @@ namespace Opm { const SolutionState& state ) { const int phase = canph_[ actph ]; - const ADB mu = fluidViscosity(phase, state.pressure, state.Rs, cells_); + const ADB mu = fluidViscosity(phase, state.pressure, state.rs, cells_); rq_[ actph ].mob = kr[ phase ] / mu; - const ADB rho = fluidDensity(phase, state.pressure, state.Rs, cells_); + const ADB rho = fluidDensity(phase, state.pressure, state.rs, cells_); const ADB gflux = grav_ * rho; ADB& head = rq_[ actph ].head; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index 3a8e571c9..adab07bbc 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -79,7 +79,7 @@ namespace Opm { SolutionState(const int np); ADB pressure; std::vector saturation; - ADB Rs; + ADB rs; ADB bhp; }; From 765ce23c3e8d9a75ce6f172985614ede55c96a5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 08:58:30 +0200 Subject: [PATCH 3/5] Work in progress on well flux equations. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 44 ++++++++++++-------- opm/autodiff/FullyImplicitBlackoilSolver.hpp | 3 +- 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 747d0f562..8d805851f 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -198,7 +198,7 @@ namespace Opm { , rq_ (fluid.numPhases()) , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), ADB::null(), - std::vector(fluid.numPhases(), ADB::null()), + ADB::null(), ADB::null() } ) { } @@ -275,6 +275,7 @@ namespace Opm { : pressure ( ADB::null()) , saturation(np, ADB::null()) , rs ( ADB::null()) + , qs ( ADB::null()) , bhp ( ADB::null()) { } @@ -326,6 +327,7 @@ namespace Opm { // water saturation (if water present) // gas saturation (if gas present) // gas solution factor (if both gas and oil present) + // well rates per active phase and well // well bottom-hole pressure // Note that oil is assumed to always be present, but is never // a primary variable. @@ -335,6 +337,7 @@ namespace Opm { if (gasandoil) { bpat.push_back(nc); } + bpat.push_back(xw.bhp().size() * np); bpat.push_back(xw.bhp().size()); SolutionState state(np); @@ -379,6 +382,11 @@ namespace Opm { state.rs = ADB::constant(Rs, bpat); } + // Well rates. + assert (not xw.wellRates().empty()); + const V qs = Eigen::Map(& xw.wellRates()[0], xw.wellRates().size()); + state.qs = ADB::constant(qs, bpat); + // Well bottom-hole pressure. assert (not xw.bhp().empty()); const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); @@ -427,6 +435,11 @@ namespace Opm { vars0.push_back(rs); } + // Initial well rates. + assert (not xw.wellRates().empty()); + const V qs = Eigen::Map(& xw.wellRates()[0], xw.wellRates().size()); + vars0.push_back(qs); + // Initial well bottom-hole pressure. assert (not xw.bhp().empty()); const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); @@ -470,8 +483,11 @@ namespace Opm { state.rs = ADB::constant(V::Zero(nc), bpat); } + // Qs. + state.qs = vars[ nextvar++ ]; + // Bhp. - state.bhp = vars[ nextvar ++]; + state.bhp = vars[ nextvar++ ]; ASSERT(nextvar == int(vars.size())); @@ -646,11 +662,7 @@ namespace Opm { const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); // DUMP(nkgradp_well); const Selector cell_to_well_selector(nkgradp_well.value()); - ADB qs = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern()); - // We can safely use a dummy rs here (for well calculations) - // as long as we do not inject oil. - const ADB rs_perfwell = ADB::constant(V::Zero(nperf), state.bhp.blockPattern()); - const std::vector well_kr = computeRelPermWells(state, well_s, well_cells); + 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); @@ -669,7 +681,7 @@ namespace Opm { 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]; - qs += superset(well_rates, Span(nw, 1, phase*nw), nw*np); + 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); @@ -680,10 +692,14 @@ namespace Opm { const int oilpos = pu.phase_pos[Oil]; const int gaspos = pu.phase_pos[Gas]; const ADB rs_perf = subset(state.rs, well_cells); - qs += superset(well_perf_rates[oilpos]*rs_perf, Span(nw, 1, gaspos*nw), nw*np); + well_rates_all += superset(well_perf_rates[oilpos]*rs_perf, Span(nw, 1, gaspos*nw), nw*np); // DUMP(well_contribs[gaspos] + well_contribs[oilpos]*state.rs); residual_.mass_balance[gaspos] += well_contribs[oilpos]*state.rs; } + + // Set the well flux equation + residual_.well_flux_eq = state.qs - well_rates_all; + // Handling BHP and SURFACE_RATE wells. V bhp_targets(nw); V rate_targets(nw); @@ -704,7 +720,7 @@ namespace Opm { } } const ADB bhp_residual = bhp - bhp_targets; - const ADB rate_residual = rate_distr * qs - rate_targets; + const ADB rate_residual = rate_distr * state.qs - rate_targets; // Choose bhp residual for positive bhp targets. Selector bhp_selector(bhp_targets); residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual); @@ -987,13 +1003,7 @@ namespace Opm { if (active_[Oil] && active_[Gas]) { r = std::max(r, residual_.rs_or_sg_eq.value().matrix().norm()); } - for (std::vector::const_iterator - b = residual_.well_flux_eq.begin(), - e = residual_.well_flux_eq.end(); - b != e; ++b) - { - r = std::max(r, (*b).value().matrix().norm()); - } + r = std::max(r, residual_.well_flux_eq.value().matrix().norm()); r = std::max(r, residual_.well_eq.value().matrix().norm()); return r; diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.hpp b/opm/autodiff/FullyImplicitBlackoilSolver.hpp index adab07bbc..1c6e4b874 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.hpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.hpp @@ -80,6 +80,7 @@ namespace Opm { ADB pressure; std::vector saturation; ADB rs; + ADB qs; ADB bhp; }; @@ -116,7 +117,7 @@ namespace Opm { struct { std::vector mass_balance; ADB rs_or_sg_eq; // Only used if both gas and oil present - std::vector well_flux_eq; + ADB well_flux_eq; ADB well_eq; } residual_; From b23e622a611201316eccce6b80588d3900f1d684 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 08:59:00 +0200 Subject: [PATCH 4/5] Remove accidentally added file. --- opm/autodiff/.#FullyImplicitBlackoilSolver.hpp | 1 - 1 file changed, 1 deletion(-) delete mode 120000 opm/autodiff/.#FullyImplicitBlackoilSolver.hpp diff --git a/opm/autodiff/.#FullyImplicitBlackoilSolver.hpp b/opm/autodiff/.#FullyImplicitBlackoilSolver.hpp deleted file mode 120000 index a680f0a09..000000000 --- a/opm/autodiff/.#FullyImplicitBlackoilSolver.hpp +++ /dev/null @@ -1 +0,0 @@ -atgeirr@arwen.lan.6100 \ No newline at end of file From 0376cb0ffff8b54c1925ed9c04008e3021d870df Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 2 Jun 2013 23:50:21 +0200 Subject: [PATCH 5/5] Initialize and update qs primary variable. --- opm/autodiff/FullyImplicitBlackoilSolver.cpp | 38 ++++++++++++++++---- 1 file changed, 31 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/FullyImplicitBlackoilSolver.cpp b/opm/autodiff/FullyImplicitBlackoilSolver.cpp index 8d805851f..f28525d2b 100644 --- a/opm/autodiff/FullyImplicitBlackoilSolver.cpp +++ b/opm/autodiff/FullyImplicitBlackoilSolver.cpp @@ -384,7 +384,12 @@ namespace Opm { // Well rates. assert (not xw.wellRates().empty()); - const V qs = Eigen::Map(& xw.wellRates()[0], xw.wellRates().size()); + // Need to reshuffle well rates, from ordered by wells, then phase, + // to ordered by phase, then wells. + const int nw = wells_.number_of_wells; + // The transpose() below switches the ordering. + const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); + const V qs = Eigen::Map(wrates.data(), nw*np); state.qs = ADB::constant(qs, bpat); // Well bottom-hole pressure. @@ -437,7 +442,12 @@ namespace Opm { // Initial well rates. assert (not xw.wellRates().empty()); - const V qs = Eigen::Map(& xw.wellRates()[0], xw.wellRates().size()); + // Need to reshuffle well rates, from ordered by wells, then phase, + // to ordered by phase, then wells. + const int nw = wells_.number_of_wells; + // The transpose() below switches the ordering. + const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); + const V qs = Eigen::Map(wrates.data(), nw*np); vars0.push_back(qs); // Initial well bottom-hole pressure. @@ -698,7 +708,8 @@ namespace Opm { } // Set the well flux equation - residual_.well_flux_eq = state.qs - well_rates_all; + residual_.well_flux_eq = state.qs + well_rates_all; + // DUMP(residual_.well_flux_eq); // Handling BHP and SURFACE_RATE wells. V bhp_targets(nw); @@ -724,7 +735,7 @@ namespace Opm { // Choose bhp residual for positive bhp targets. Selector bhp_selector(bhp_targets); residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual); - DUMP(residual_.well_eq); + // DUMP(residual_.well_eq); } @@ -741,8 +752,9 @@ namespace Opm { if (active_[Oil] && active_[Gas]) { mass_res = vertcat(mass_res, residual_.rs_or_sg_eq); } - const ADB total_residual = collapseJacs(vertcat(mass_res, residual_.well_eq)); - DUMP(total_residual); + const ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq); + const ADB total_residual = collapseJacs(vertcat(mass_res, well_res)); + // DUMP(total_residual); const Eigen::SparseMatrix matr = total_residual.derivative()[0]; @@ -792,6 +804,8 @@ namespace Opm { varstart += dsg.size(); const V drs = (active_[Water] && active_[Gas]) ? subset(dx, Span(nc, 1, varstart)) : null; varstart += drs.size(); + const V dqs = subset(dx, Span(np*nw, 1, varstart)); + varstart += dqs.size(); const V dbhp = subset(dx, Span(nw, 1, varstart)); varstart += dbhp.size(); ASSERT(varstart == dx.size()); @@ -889,10 +903,20 @@ namespace Opm { } } + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const V dwr = Eigen::Map(wwr.data(), nw*np); + const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const V wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + // Bhp update. const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); const V bhp = bhp_old - dbhp; - std::copy(&bhp[0], &bhp[0] + nw, well_state.bhp().begin()); + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); }