From 70ebd21cfd94e73fe3576fc7ce5799647ec9df9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 23 May 2013 15:08:21 +0200 Subject: [PATCH] Moved some calculations to computeExplicitData() method. --- opm/autodiff/ImpesTPFAAD.cpp | 116 +++++++++++++++++++++++++++++++---- opm/autodiff/ImpesTPFAAD.hpp | 5 ++ 2 files changed, 109 insertions(+), 12 deletions(-) diff --git a/opm/autodiff/ImpesTPFAAD.cpp b/opm/autodiff/ImpesTPFAAD.cpp index 916f98720..49bd1d165 100644 --- a/opm/autodiff/ImpesTPFAAD.cpp +++ b/opm/autodiff/ImpesTPFAAD.cpp @@ -28,6 +28,13 @@ #include + +// Repeated from inside ImpesTPFAAD for convenience. +typedef AutoDiff::ForwardBlock ADB; +typedef ADB::V V; +typedef ADB::M M; + + namespace { std::vector buildAllCells(const int nc) @@ -86,14 +93,37 @@ namespace { return G; } + + V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav) + { + const int nw = wells.number_of_wells; + const int nperf = wells.well_connpos[nw]; + const int dim = grid.dimensions; + V wdp = V::Zero(nperf,1); + ASSERT(wdp.size() == rho.size()); + + // Main loop, iterate over all perforations, + // using the following formula: + // wdp(perf) = g*(perf_z - well_ref_z)*rho(perf) + // where the total density rho(perf) is taken to be + // sum_p (rho_p*saturation_p) in the perforation cell. + // [although this is computed on the outside of this function]. + for (int w = 0; w < nw; ++w) { + const double ref_depth = wells.depth_ref[w]; + for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) { + const int cell = wells.well_cells[j]; + const double cell_depth = grid.cell_centroids[dim * cell + dim - 1]; + wdp[j] = rho[j]*grav*(cell_depth - ref_depth); + } + } + return wdp; + } + } // anonymous namespace namespace Opm { - // Repeated from inside ImpesTPFAAD for convenience. - typedef AutoDiff::ForwardBlock ADB; - typedef ADB::V V; - typedef ADB::M M; + ImpesTPFAAD::ImpesTPFAAD(const UnstructuredGrid& grid, @@ -113,9 +143,14 @@ namespace Opm { , well_flow_residual_ () , well_residual_ (ADB::null()) , total_residual_ (ADB::null()) + , qs_ (ADB::null()) { } + + + + void ImpesTPFAAD::solve(const double dt, BlackoilState& state, @@ -126,6 +161,8 @@ namespace Opm { well_flow_residual_.resize(np, ADB::null()); + // Compute dynamic data that are treated explicitly. + computeExplicitData(dt, state, well_state); // Compute relperms once and for all (since saturations are explicit). DataBlock s = Eigen::Map(state.saturation().data(), nc, np); ASSERT(np == 2); @@ -163,12 +200,12 @@ namespace Opm { const double r = residualNorm(); - std::cout << std::setw(9) << it - << std::setw(18) << r << std::endl; - resTooLarge = (r > atol) && (r > rtol*r0); it += 1; + + std::cout << std::setw(9) << it + << std::setw(18) << r << std::endl; } if (resTooLarge) { @@ -179,6 +216,64 @@ namespace Opm { } } + + + + + void + ImpesTPFAAD::computeExplicitData(const double dt, + const BlackoilState& state, + const WellState& well_state) + { + const int nc = grid_.number_of_cells; + const int np = state.numPhases(); + const int nw = wells_.number_of_wells; + const int nperf = wells_.well_connpos[nw]; + const int dim = grid_.dimensions; + + const std::vector cells = buildAllCells(nc); + + // Compute relperms. + 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. + 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); + + // Compute well pressure differentials. + // Construct pressure difference vector for wells. + 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); + } + } + V cell_rho_total = V::Zero(nc,1); + const Eigen::Map p(state.pressure().data(), nc, 1); + for (int phase = 0; phase < np; ++phase) { + const V cell_rho = fluidRho(phase, p, cells); + const V cell_s = s.col(phase); + cell_rho_total += cell_s * cell_rho; + } + V rho_perf = subset(cell_rho_total, well_cells); + well_perf_dp_ = computePerfPress(grid_, wells_, rho_perf, g ? g[dim-1] : 0.0); + } + + + + + void ImpesTPFAAD::assemble(const double dt, const BlackoilState& state, @@ -231,10 +326,8 @@ namespace Opm { } well_to_perf.setFromTriplets(w2p.begin(), w2p.end()); const M perf_to_well = well_to_perf.transpose(); - // Construct pressure difference vector for wells. - const V well_perf_dp = V::Zero(well_cells.size()); // No gravity yet! // Finally construct well perforation pressures and well flows. - const ADB p_perfwell = well_to_perf*bhp + well_perf_dp; + 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()); @@ -406,8 +499,7 @@ namespace Opm { ops_.internal_faces); const V nkgradp = transi * (ops_.ngrad * p.matrix()).array(); - const V well_perf_dp = V::Zero(well_cells.size()); // No gravity yet! - const V p_perfwell = (well_to_perf*bhp.matrix()).array() + well_perf_dp; + const V p_perfwell = (well_to_perf*bhp.matrix()).array() + well_perf_dp_; const V nkgradp_well = transw * (p_perfcell - p_perfwell); const Selector cell_to_well_selector(nkgradp_well); diff --git a/opm/autodiff/ImpesTPFAAD.hpp b/opm/autodiff/ImpesTPFAAD.hpp index 65917c86a..298f0543a 100644 --- a/opm/autodiff/ImpesTPFAAD.hpp +++ b/opm/autodiff/ImpesTPFAAD.hpp @@ -88,8 +88,13 @@ namespace Opm { ADB total_residual_; std::vector kr_; std::vector well_kr_; + ADB qs_; + V well_perf_dp_; // Methods for assembling and solving. + void computeExplicitData(const double dt, + const BlackoilState& state, + const WellState& well_state); void assemble(const double dt, const BlackoilState& state, const WellState& well_state);