diff --git a/opm/autodiff/ImpesTPFAAD.hpp b/opm/autodiff/ImpesTPFAAD.hpp index 6ecab0c9e..a2ff4dc2f 100644 --- a/opm/autodiff/ImpesTPFAAD.hpp +++ b/opm/autodiff/ImpesTPFAAD.hpp @@ -33,6 +33,7 @@ #include #include #include +#include #include @@ -129,9 +130,24 @@ namespace Opm { { const int nc = grid_.number_of_cells; const int np = state.numPhases(); + + // Compute relperms once and for all (since saturations are explicit). DataBlock s = Eigen::Map(state.saturation().data(), nc, np); ASSERT(np == 2); kr_ = fluid_.relperm(s.col(0), s.col(1), V::Zero(nc,1), buildAllCells(nc)); + // Compute relperms for wells. This must be revisited for crossflow. + const int nw = wells_.number_of_wells; + const int nperf = wells_.well_connpos[nw]; + DataBlock well_s(nperf, np); + for (int w = 0; w < nw; ++w) { + const double* comp_frac = &wells_.comp_frac[np*w]; + for (int j = wells_.well_connpos[w]; j < wells_.well_connpos[w+1]; ++j) { + well_s.row(j) = Eigen::Map(comp_frac, 1, np); + } + } + const std::vector well_cells(wells_.well_cells, + wells_.well_cells + nperf); + well_kr_ = fluid_.relperm(well_s.col(0), well_s.col(1), V::Zero(nperf,1), well_cells); const double atol = 1.0e-15; const double rtol = 5.0e-10; @@ -177,22 +193,26 @@ namespace Opm { Eigen::Dynamic, Eigen::RowMajor> DataBlock; - const UnstructuredGrid& grid_; - const BOFluid& fluid_; - const GeoProps& geo_ ; - const Wells& wells_; + const UnstructuredGrid& grid_; + const BOFluid& fluid_; + const GeoProps& geo_ ; + const Wells& wells_; const LinearSolverInterface& linsolver_; - // PDepFData pdepfdata_; - HelperOps ops_; - const M grav_; - ADB cell_residual_; - ADB well_residual_; - std::vector kr_; + HelperOps ops_; + const M grav_; + ADB cell_residual_; + ADB well_residual_; + std::vector kr_; + std::vector well_kr_; enum { Water = BOFluid::Water, Oil = BOFluid::Oil, Gas = BOFluid::Gas }; + + + + void assemble(const double dt, const BlackoilState& state, @@ -203,8 +223,7 @@ namespace Opm { const int nc = grid_.number_of_cells; const int np = state.numPhases(); const int nw = wells_.number_of_wells; - - // pdepfdata_.computePressQuant(state); + const int nperf = wells_.well_connpos[nw]; const std::vector cells = buildAllCells(nc); @@ -213,10 +232,9 @@ namespace Opm { const V delta_t = dt * V::Ones(nc, 1); const V transi = subset(geo_.transmissibility(), ops_.internal_faces); - const int num_perf = wells_.well_connpos[nw]; const std::vector well_cells(wells_.well_cells, - wells_.well_cells + num_perf); - const V transw = Eigen::Map(wells_.WI, num_perf, 1); + wells_.well_cells + nperf); + const V transw = Eigen::Map(wells_.WI, nperf, 1); // Initialize AD variables: p (cell pressures) and bhp (well bhp). const V p0 = Eigen::Map(&state.pressure()[0], nc, 1); @@ -246,10 +264,10 @@ namespace Opm { well_to_perf.setFromTriplets(w2p.begin(), w2p.end()); // Construct pressure difference vector for wells. const V well_perf_dp = V::Zero(well_cells.size()); // No gravity yet! - // Finally construct well perforation pressures. + // Finally construct well perforation pressures and well flows. const ADB p_perfwell = well_to_perf*bhp + well_perf_dp; - const ADB nkgradp_well = transw * (p_perfcell - p_perfwell); + const Selector cell_to_well_selector(nkgradp_well.value()); cell_residual_ = ADB::constant(pv, bpat); #define COMPENSATE_FLOAT_PRECISION 0 @@ -259,22 +277,35 @@ namespace Opm { for (int phase = 0; phase < np; ++phase) { const ADB cell_b = fluidFvf(phase, p, cells); const ADB cell_rho = fluidRho(phase, p, cells); + const ADB well_b = fluidFvf(phase, p_perfwell, well_cells); const V kr = fluidKr(phase); + // Explicitly not asking for derivatives of viscosity, + // since they are not available yet. const V mu = fluidMu(phase, p.value(), cells); - const V mf = upwind.select(kr / mu); - const ADB flux = mf * (nkgradp + (grav_ * cell_rho)); + const V cell_mob = kr / mu; + const V face_mob = upwind.select(cell_mob); + const V well_kr = fluidKrWell(phase); + const V well_mu = fluidMu(phase, p_perfwell.value(), well_cells); + const V well_mob = well_kr / well_mu; + const V perf_mob = cell_to_well_selector.select(subset(cell_mob, well_cells), well_mob); + const ADB flux = face_mob * (nkgradp + (grav_ * cell_rho)); + const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations. const ADB face_b = upwind.select(cell_b); + const ADB perf_b = cell_to_well_selector.select(subset(cell_b, well_cells), well_b); const V z0 = z0all.block(0, phase, nc, 1); const V q = qall .block(0, phase, nc, 1); + const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc); #if COMPENSATE_FLOAT_PRECISION - const ADB divcontrib = delta_t * (ops_.div * (flux * face_b)); + const ADB divcontrib = delta_t * (ops_.div * (flux * face_b) + + well_contrib); const V qcontrib = delta_t * q; const ADB pvcontrib = ADB::constant(pv*z0, bpat); const ADB component_contrib = pvcontrib + qcontrib; divcontrib_sum = divcontrib_sum - divcontrib/cell_b; cell_residual_ = cell_residual_ - (component_contrib/cell_b); #else - const ADB component_contrib = pv*z0 + delta_t*(q - (ops_.div * (flux * face_b))); + const ADB component_contrib = pv*z0 + + delta_t*(q - (ops_.div * (flux * face_b) + well_contrib)); cell_residual_ = cell_residual_ - (component_contrib / cell_b); #endif } @@ -283,6 +314,10 @@ namespace Opm { #endif } + + + + void solveJacobianSystem(BlackoilState& state) const { @@ -306,12 +341,20 @@ namespace Opm { std::copy(&p[0], &p[0] + nc, state.pressure().begin()); } + + + + double residualNorm() const { return cell_residual_.value().matrix().norm(); } + + + + void computeFluxes(BlackoilState& state) const { @@ -344,6 +387,9 @@ namespace Opm { } + + + V fluidMu(const int phase, const V& p, const std::vector& cells) const { switch (phase) { @@ -359,6 +405,11 @@ namespace Opm { THROW("Unknown phase index " << phase); } } + + + + + ADB fluidMu(const int phase, const ADB& p, const std::vector& cells) const { switch (phase) { @@ -375,6 +426,10 @@ namespace Opm { } } + + + + ADB fluidFvf(const int phase, const ADB& p, const std::vector& cells) const { switch (phase) { @@ -391,6 +446,10 @@ namespace Opm { } } + + + + ADB fluidRho(const int phase, const ADB& p, const std::vector& cells) const { const double* rhos = fluid_.surfaceDensity(); @@ -400,10 +459,23 @@ namespace Opm { } + + + V fluidKr(const int phase) const { return kr_[phase]; } + + + + + + V fluidKrWell(const int phase) const + { + return well_kr_[phase]; + } + }; } // namespace Opm