From 3e8bb53730f60414a60c4b13dc3beb1674377f34 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 22 May 2013 15:49:55 +0200 Subject: [PATCH] Moved implementation of class ImpesTPFAAD to separate file. Also cleaned up header usage in ImpesTPFAAD.hpp, making some new inclusions necessary in test program. --- CMakeLists_files.cmake | 1 + examples/test_impestpfa_ad.cpp | 5 +- opm/autodiff/ImpesTPFAAD.cpp | 486 +++++++++++++++++++++++++++++++ opm/autodiff/ImpesTPFAAD.hpp | 510 +++------------------------------ 4 files changed, 533 insertions(+), 469 deletions(-) create mode 100644 opm/autodiff/ImpesTPFAAD.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 7a415c6a6..29f89feb3 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -28,6 +28,7 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/BlackoilPropsAd.cpp opm/autodiff/BlackoilPropsAdInterface.cpp + opm/autodiff/ImpesTPFAAD.cpp opm/autodiff/SimulatorIncompTwophaseAdfi.cpp opm/autodiff/TransportSolverTwophaseAd.cpp ) diff --git a/examples/test_impestpfa_ad.cpp b/examples/test_impestpfa_ad.cpp index 2acd6b8b4..c53f7aa06 100644 --- a/examples/test_impestpfa_ad.cpp +++ b/examples/test_impestpfa_ad.cpp @@ -20,8 +20,6 @@ #include -#define HACK_INCOMPRESSIBLE_GRAVITY 0 - #include #include #include @@ -36,10 +34,11 @@ #include #include +#include +#include #include #include -// #include #include diff --git a/opm/autodiff/ImpesTPFAAD.cpp b/opm/autodiff/ImpesTPFAAD.cpp new file mode 100644 index 000000000..50283f242 --- /dev/null +++ b/opm/autodiff/ImpesTPFAAD.cpp @@ -0,0 +1,486 @@ +/* + Copyright 2013 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include + +#include +#include +#include +#include +#include + +namespace { + std::vector + buildAllCells(const int nc) + { + std::vector all_cells(nc); + + for (int c = 0; c < nc; ++c) { all_cells[c] = c; } + + return all_cells; + } + + template + AutoDiff::ForwardBlock::M + gravityOperator(const UnstructuredGrid& grid, + const HelperOps& ops , + const GeoProps& geo ) + { + const int nc = grid.number_of_cells; + + std::vector f2hf(2 * grid.number_of_faces, -1); + for (int c = 0, i = 0; c < nc; ++c) { + for (; i < grid.cell_facepos[c + 1]; ++i) { + const int f = grid.cell_faces[ i ]; + const int p = 0 + (grid.face_cells[2*f + 0] != c); + + f2hf[2*f + p] = i; + } + } + + typedef AutoDiff::ForwardBlock::V V; + typedef AutoDiff::ForwardBlock::M M; + + const V& gpot = geo.gravityPotential(); + const V& trans = geo.transmissibility(); + + const HelperOps::IFaces::Index ni = ops.internal_faces.size(); + + typedef Eigen::Triplet Tri; + std::vector grav; grav.reserve(2 * ni); + for (HelperOps::IFaces::Index i = 0; i < ni; ++i) { + const int f = ops.internal_faces[ i ]; + const int c1 = grid.face_cells[2*f + 0]; + const int c2 = grid.face_cells[2*f + 1]; + + assert ((c1 >= 0) && (c2 >= 0)); + + const double dG1 = gpot[ f2hf[2*f + 0] ]; + const double dG2 = gpot[ f2hf[2*f + 1] ]; + const double t = trans[ f ]; + + grav.push_back(Tri(i, c1, t * dG1)); + grav.push_back(Tri(i, c2, - t * dG2)); + } + + M G(ni, nc); G.setFromTriplets(grav.begin(), grav.end()); + + return G; + } +} // 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, + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo, + const Wells& wells, + const LinearSolverInterface& linsolver) + : grid_ (grid) + , fluid_ (fluid) + , geo_ (geo) + , wells_ (wells) + , linsolver_(linsolver) + // , pdepfdata_(grid.number_of_cells, fluid) + , ops_ (grid) + , grav_ (gravityOperator(grid_, ops_, geo_)) + , cell_residual_ (ADB::null()) + , well_flow_residual_ () + , well_residual_ (ADB::null()) + , total_residual_ (ADB::null()) + { + } + + void + ImpesTPFAAD::solve(const double dt, + BlackoilState& state, + WellState& well_state) + { + const int nc = grid_.number_of_cells; + const int np = state.numPhases(); + + well_flow_residual_.resize(np, ADB::null()); + + // 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; + const int maxit = 15; + + assemble(dt, state, well_state); + + const double r0 = residualNorm(); + int it = 0; + bool resTooLarge = r0 > atol; + while (resTooLarge && (it < maxit)) { + solveJacobianSystem(state, well_state); + + assemble(dt, state, well_state); + + const double r = residualNorm(); + + resTooLarge = (r > atol) && (r > rtol*r0); + + it += 1; + } + + if (resTooLarge) { + THROW("Failed to compute converged pressure solution"); + } + else { + computeFluxes(state); + } + } + + void + ImpesTPFAAD::assemble(const double dt, + const BlackoilState& state, + const WellState& well_state) + { + + const V& pv = geo_.poreVolume(); + 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 std::vector cells = buildAllCells(nc); + + const Eigen::Map z0all(&state.surfacevol()[0], nc, np); + const DataBlock qall = DataBlock::Zero(nc, np); + const V delta_t = dt * V::Ones(nc, 1); + const V transi = subset(geo_.transmissibility(), + ops_.internal_faces); + const std::vector well_cells(wells_.well_cells, + wells_.well_cells + nperf); + const V transw = Eigen::Map(wells_.WI, nperf, 1); + + // Initialize AD variables: p (cell pressures), qs (well rates) and bhp (well bhp). + const V p0 = Eigen::Map(&state.pressure()[0], nc, 1); + const V qs0 = Eigen::Map(&well_state.wellRates()[0], nw*np, 1); + const V bhp0 = Eigen::Map(&well_state.bhp()[0], nw, 1); + std::vector vars0 = { p0, qs0, bhp0 }; + std::vector vars = ADB::variables(vars0); + const ADB& p = vars[0]; + const ADB& qs = vars[1]; + const ADB& bhp = vars[2]; + std::vector bpat = p.blockPattern(); + + // Compute T_ij * (p_i - p_j) and use for upwinding. + const ADB nkgradp = transi * (ops_.ngrad * p); + const UpwindSelector upwind(grid_, ops_, nkgradp.value()); + + // Extract variables for perforation cell pressures + // and corresponding perforation well pressures. + const ADB p_perfcell = subset(p, well_cells); + // Construct matrix to map wells->perforations. + M well_to_perf(well_cells.size(), nw); + typedef Eigen::Triplet Tri; + std::vector w2p; + for (int w = 0; w < nw; ++w) { + for (int perf = wells_.well_connpos[w]; perf < wells_.well_connpos[w+1]; ++perf) { + w2p.emplace_back(perf, w, 1.0); + } + } + 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 nkgradp_well = transw * (p_perfcell - p_perfwell); + const Selector cell_to_well_selector(nkgradp_well.value()); + + cell_residual_ = ADB::constant(pv, bpat); + well_residual_ = ADB::constant(V::Zero(nw,1), bpat); + ADB divcontrib_sum = ADB::constant(V::Zero(nc,1), bpat); + 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 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); + 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); + const ADB well_rates = perf_to_well * (perf_flux*perf_b); + std::vector well_flow_res_phase_idx(nw); + for (int w = 0; w < nw; ++w) { + well_flow_res_phase_idx[w] = w + phase*nw; + } + well_flow_residual_[phase] = well_rates - subset(qs, well_flow_res_phase_idx); + } + cell_residual_ = cell_residual_ + divcontrib_sum; + // Handling BHP and SURFACE_RATE wells. + V bhp_targets(nw,1); + V rate_targets(nw,1); + M rate_distr(nw, np*nw); + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells_.ctrls[w]; + if (wc->type[wc->current] == BHP) { + bhp_targets[w] = wc->target[wc->current]; + rate_targets[w] = -1e100; + } else if (wc->type[wc->current] == SURFACE_RATE) { + bhp_targets[w] = -1e100; + rate_targets[w] = wc->target[wc->current]; + for (int phase = 0; phase < np; ++phase) { + rate_distr.insert(w, phase*nw + w) = wc->distr[phase]; + } + } else { + THROW("Can only handle BHP and SURFACE_RATE type controls."); + } + } + const ADB bhp_residual = bhp - bhp_targets; + const ADB rate_residual = rate_distr * qs - rate_targets; + // Choose bhp residual for positive bhp targets. + Selector bhp_selector(bhp_targets); + well_residual_ = bhp_selector.select(bhp_residual, rate_residual); + + ASSERT(np == 2); + const ADB well_flow_res = vertcat(well_flow_residual_[0], well_flow_residual_[1]); + const ADB well_res = vertcat(well_flow_res, well_residual_); + total_residual_ = collapseJacs(vertcat(cell_residual_, well_res)); + } + + + + + + void + ImpesTPFAAD::solveJacobianSystem(BlackoilState& state, + WellState& well_state) const + { + const int nc = grid_.number_of_cells; + const int nw = wells_.number_of_wells; + const int np = state.numPhases(); + + Eigen::SparseMatrix matr = total_residual_.derivative()[0]; + + V dx(total_residual_.size()); + Opm::LinearSolverInterface::LinearSolverReport rep + = linsolver_.solve(matr.rows(), matr.nonZeros(), + matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), + total_residual_.value().data(), dx.data()); + if (!rep.converged) { + THROW("ImpesTPFAAD::solve(): Linear solver convergence failure."); + } + const V p0 = Eigen::Map(&state.pressure()[0], nc, 1); + const V dp = subset(dx, buildAllCells(nc)); + const V p = p0 - dp; + std::copy(&p[0], &p[0] + nc, state.pressure().begin()); + const V qs0 = Eigen::Map(&well_state.wellRates()[0], nw*np, 1); + std::vector qs_dofs(np*nw); + for (int w = 0; w < nw; ++w) { + for (int phase = 0; phase < np; ++phase) { + qs_dofs[w*np + phase] = nc + w*np + phase; + } + } + const V dqs = subset(dx, qs_dofs); + const V qs = qs0 - dqs; + std::copy(&qs[0], &qs[0] + np*nw, well_state.wellRates().begin()); + const V bhp0 = Eigen::Map(&well_state.bhp()[0], nw, 1); + std::vector bhp_dofs(nw); + for (int w = 0; w < nw; ++w) { + bhp_dofs[w] = nc + np*nw + w; + } + ASSERT(bhp_dofs.back() + 1 == total_residual_.size()); + const V dbhp = subset(dx, bhp_dofs); + const V bhp = bhp0 - dbhp; + std::copy(&bhp[0], &bhp[0] + nw, well_state.bhp().begin()); + } + + + + + + double + ImpesTPFAAD::residualNorm() const + { + return total_residual_.value().matrix().norm(); + } + + + + + + void + ImpesTPFAAD::computeFluxes(BlackoilState& state) const + { + const int nc = grid_.number_of_cells; + const int np = state.numPhases(); + + const std::vector cells = buildAllCells(nc); + + const V p0 = Eigen::Map(&state.pressure()[0], nc, 1); + const ADB p = ADB::constant(p0, cell_residual_.blockPattern()); + + const V transi = subset(geo_.transmissibility(), + ops_.internal_faces); + const V nkgradp = transi * (ops_.ngrad * p0.matrix()).array(); + + V flux = V::Zero(ops_.internal_faces.size(), 1); + + for (int phase = 0; phase < np; ++phase) { + const ADB cell_rho = fluidRho(phase, p, cells); + const V head = nkgradp + (grav_ * cell_rho.value().matrix()).array(); + const UpwindSelector upwind(grid_, ops_, head); + const V kr = fluidKr(phase); + const V mu = fluidMu(phase, p.value(), cells); + const V mf = upwind.select(kr / mu); + + flux += mf * head; + } + V all_flux = superset(flux, ops_.internal_faces, grid_.number_of_faces); + std::copy(all_flux.data(), all_flux.data() + grid_.number_of_faces, state.faceflux().data()); + } + + + + + + V ImpesTPFAAD::fluidMu(const int phase, const V& p, const std::vector& cells) const + { + switch (phase) { + case Water: + return fluid_.muWat(p, cells); + case Oil: { + V dummy_rs = V::Zero(p.size(), 1) * p; + return fluid_.muOil(p, dummy_rs, cells); + } + case Gas: + return fluid_.muGas(p, cells); + default: + THROW("Unknown phase index " << phase); + } + } + + + + + + ADB ImpesTPFAAD::fluidMu(const int phase, const ADB& p, const std::vector& cells) const + { + switch (phase) { + case Water: + return fluid_.muWat(p, cells); + case Oil: { + ADB dummy_rs = V::Zero(p.size(), 1) * p; + return fluid_.muOil(p, dummy_rs, cells); + } + case Gas: + return fluid_.muGas(p, cells); + default: + THROW("Unknown phase index " << phase); + } + } + + + + + + ADB ImpesTPFAAD::fluidFvf(const int phase, const ADB& p, const std::vector& cells) const + { + switch (phase) { + case Water: + return fluid_.bWat(p, cells); + case Oil: { + ADB dummy_rs = V::Zero(p.size(), 1) * p; + return fluid_.bOil(p, dummy_rs, cells); + } + case Gas: + return fluid_.bGas(p, cells); + default: + THROW("Unknown phase index " << phase); + } + } + + + + + + ADB ImpesTPFAAD::fluidRho(const int phase, const ADB& p, const std::vector& cells) const + { + const double* rhos = fluid_.surfaceDensity(); + ADB b = fluidFvf(phase, p, cells); + ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b; + return rho; + } + + + + + + V ImpesTPFAAD::fluidKr(const int phase) const + { + return kr_[phase]; + } + + + + + + V ImpesTPFAAD::fluidKrWell(const int phase) const + { + return well_kr_[phase]; + } + +} // namespace Opm + diff --git a/opm/autodiff/ImpesTPFAAD.hpp b/opm/autodiff/ImpesTPFAAD.hpp index 9d2695e8d..c6b640550 100644 --- a/opm/autodiff/ImpesTPFAAD.hpp +++ b/opm/autodiff/ImpesTPFAAD.hpp @@ -25,168 +25,44 @@ #include #include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include - struct UnstructuredGrid; - -namespace { - std::vector - buildAllCells(const int nc) - { - std::vector all_cells(nc); - - for (int c = 0; c < nc; ++c) { all_cells[c] = c; } - - return all_cells; - } - - template - AutoDiff::ForwardBlock::M - gravityOperator(const UnstructuredGrid& grid, - const HelperOps& ops , - const GeoProps& geo ) - { - const int nc = grid.number_of_cells; - - std::vector f2hf(2 * grid.number_of_faces, -1); - for (int c = 0, i = 0; c < nc; ++c) { - for (; i < grid.cell_facepos[c + 1]; ++i) { - const int f = grid.cell_faces[ i ]; - const int p = 0 + (grid.face_cells[2*f + 0] != c); - - f2hf[2*f + p] = i; - } - } - - typedef AutoDiff::ForwardBlock::V V; - typedef AutoDiff::ForwardBlock::M M; - - const V& gpot = geo.gravityPotential(); - const V& trans = geo.transmissibility(); - - const HelperOps::IFaces::Index ni = ops.internal_faces.size(); - - typedef Eigen::Triplet Tri; - std::vector grav; grav.reserve(2 * ni); - for (HelperOps::IFaces::Index i = 0; i < ni; ++i) { - const int f = ops.internal_faces[ i ]; - const int c1 = grid.face_cells[2*f + 0]; - const int c2 = grid.face_cells[2*f + 1]; - - assert ((c1 >= 0) && (c2 >= 0)); - - const double dG1 = gpot[ f2hf[2*f + 0] ]; - const double dG2 = gpot[ f2hf[2*f + 1] ]; - const double t = trans[ f ]; - - grav.push_back(Tri(i, c1, t * dG1)); - grav.push_back(Tri(i, c2, - t * dG2)); - } - - M G(ni, nc); G.setFromTriplets(grav.begin(), grav.end()); - - return G; - } -} +struct Wells; namespace Opm { + class DerivedGeology; + class LinearSolverInterface; + class BlackoilState; + class WellState; - class ImpesTPFAAD { + /// Class for solving black-oil impes problems. + /// Current known limitations: + /// - pressure solve only + /// - no miscibility + /// - no gravity in wells or crossflow + class ImpesTPFAAD + { public: + /// Construct impes solver. ImpesTPFAAD(const UnstructuredGrid& grid, const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo, const Wells& wells, - const LinearSolverInterface& linsolver) - : grid_ (grid) - , fluid_ (fluid) - , geo_ (geo) - , wells_ (wells) - , linsolver_(linsolver) - // , pdepfdata_(grid.number_of_cells, fluid) - , ops_ (grid) - , grav_ (gravityOperator(grid_, ops_, geo_)) - , cell_residual_ (ADB::null()) - , well_flow_residual_ () - , well_residual_ (ADB::null()) - , total_residual_ (ADB::null()) - { - } - - void - solve(const double dt, - BlackoilState& state, - WellState& well_state) - { - const int nc = grid_.number_of_cells; - const int np = state.numPhases(); - - well_flow_residual_.resize(np, ADB::null()); - - // 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; - const int maxit = 15; - - assemble(dt, state, well_state); - - const double r0 = residualNorm(); - int it = 0; - bool resTooLarge = r0 > atol; - while (resTooLarge && (it < maxit)) { - solveJacobianSystem(state, well_state); - - assemble(dt, state, well_state); - - const double r = residualNorm(); - - resTooLarge = (r > atol) && (r > rtol*r0); - - it += 1; - } - - if (resTooLarge) { - THROW("Failed to compute converged pressure solution"); - } - else { - computeFluxes(state); - } - } + const LinearSolverInterface& linsolver); + /// Solve forward in time. + /// Currently this will only modify + /// state.pressure(), state.faceflux(), well_state.bhp() + /// and well_state.wellRates(). + void solve(const double dt, + BlackoilState& state, + WellState& well_state); private: // Disallow copying and assignment ImpesTPFAAD(const ImpesTPFAAD& rhs); ImpesTPFAAD& operator=(const ImpesTPFAAD& rhs); + // Types typedef AutoDiff::ForwardBlock ADB; typedef ADB::V V; typedef ADB::M M; @@ -194,7 +70,11 @@ namespace Opm { Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> DataBlock; + enum { Water = BlackoilPropsAdInterface::Water, + Oil = BlackoilPropsAdInterface::Oil, + Gas = BlackoilPropsAdInterface::Gas }; + // Data const UnstructuredGrid& grid_; const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_ ; @@ -209,327 +89,25 @@ namespace Opm { std::vector kr_; std::vector well_kr_; - enum { Water = BlackoilPropsAdInterface::Water, - Oil = BlackoilPropsAdInterface::Oil, - Gas = BlackoilPropsAdInterface::Gas }; - - - void - assemble(const double dt, - const BlackoilState& state, - const WellState& well_state) - { - - const V& pv = geo_.poreVolume(); - 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 std::vector cells = buildAllCells(nc); - - const Eigen::Map z0all(&state.surfacevol()[0], nc, np); - const DataBlock qall = DataBlock::Zero(nc, np); - const V delta_t = dt * V::Ones(nc, 1); - const V transi = subset(geo_.transmissibility(), - ops_.internal_faces); - const std::vector well_cells(wells_.well_cells, - wells_.well_cells + nperf); - const V transw = Eigen::Map(wells_.WI, nperf, 1); - - // Initialize AD variables: p (cell pressures), qs (well rates) and bhp (well bhp). - const V p0 = Eigen::Map(&state.pressure()[0], nc, 1); - const V qs0 = Eigen::Map(&well_state.wellRates()[0], nw*np, 1); - const V bhp0 = Eigen::Map(&well_state.bhp()[0], nw, 1); - std::vector vars0 = { p0, qs0, bhp0 }; - std::vector vars = ADB::variables(vars0); - const ADB& p = vars[0]; - const ADB& qs = vars[1]; - const ADB& bhp = vars[2]; - std::vector bpat = p.blockPattern(); - - // Compute T_ij * (p_i - p_j) and use for upwinding. - const ADB nkgradp = transi * (ops_.ngrad * p); - const UpwindSelector upwind(grid_, ops_, nkgradp.value()); - - // Extract variables for perforation cell pressures - // and corresponding perforation well pressures. - const ADB p_perfcell = subset(p, well_cells); - // Construct matrix to map wells->perforations. - M well_to_perf(well_cells.size(), nw); - typedef Eigen::Triplet Tri; - std::vector w2p; - for (int w = 0; w < nw; ++w) { - for (int perf = wells_.well_connpos[w]; perf < wells_.well_connpos[w+1]; ++perf) { - w2p.emplace_back(perf, w, 1.0); - } - } - 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 nkgradp_well = transw * (p_perfcell - p_perfwell); - const Selector cell_to_well_selector(nkgradp_well.value()); - - cell_residual_ = ADB::constant(pv, bpat); - well_residual_ = ADB::constant(V::Zero(nw,1), bpat); - ADB divcontrib_sum = ADB::constant(V::Zero(nc,1), bpat); - 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 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); - 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); - const ADB well_rates = perf_to_well * (perf_flux*perf_b); - std::vector well_flow_res_phase_idx(nw); - for (int w = 0; w < nw; ++w) { - well_flow_res_phase_idx[w] = w + phase*nw; - } - well_flow_residual_[phase] = well_rates - subset(qs, well_flow_res_phase_idx); - } - cell_residual_ = cell_residual_ + divcontrib_sum; - // Handling BHP and SURFACE_RATE wells. - V bhp_targets(nw,1); - V rate_targets(nw,1); - M rate_distr(nw, np*nw); - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells_.ctrls[w]; - if (wc->type[wc->current] == BHP) { - bhp_targets[w] = wc->target[wc->current]; - rate_targets[w] = -1e100; - } else if (wc->type[wc->current] == SURFACE_RATE) { - bhp_targets[w] = -1e100; - rate_targets[w] = wc->target[wc->current]; - for (int phase = 0; phase < np; ++phase) { - rate_distr.insert(w, phase*nw + w) = wc->distr[phase]; - } - } else { - THROW("Can only handle BHP and SURFACE_RATE type controls."); - } - } - const ADB bhp_residual = bhp - bhp_targets; - const ADB rate_residual = rate_distr * qs - rate_targets; - // Choose bhp residual for positive bhp targets. - Selector bhp_selector(bhp_targets); - well_residual_ = bhp_selector.select(bhp_residual, rate_residual); - - ASSERT(np == 2); - const ADB well_flow_res = vertcat(well_flow_residual_[0], well_flow_residual_[1]); - const ADB well_res = vertcat(well_flow_res, well_residual_); - total_residual_ = collapseJacs(vertcat(cell_residual_, well_res)); - } - - - - - - void - solveJacobianSystem(BlackoilState& state, - WellState& well_state) const - { - const int nc = grid_.number_of_cells; - const int nw = wells_.number_of_wells; - const int np = state.numPhases(); - - Eigen::SparseMatrix matr = total_residual_.derivative()[0]; - -#if HACK_INCOMPRESSIBLE_GRAVITY - matr.coeffRef(0, 0) *= 2; -#endif - - V dx(total_residual_.size()); - Opm::LinearSolverInterface::LinearSolverReport rep - = linsolver_.solve(matr.rows(), matr.nonZeros(), - matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(), - total_residual_.value().data(), dx.data()); - if (!rep.converged) { - THROW("ImpesTPFAAD::solve(): Linear solver convergence failure."); - } - const V p0 = Eigen::Map(&state.pressure()[0], nc, 1); - const V dp = subset(dx, buildAllCells(nc)); - const V p = p0 - dp; - std::copy(&p[0], &p[0] + nc, state.pressure().begin()); - const V qs0 = Eigen::Map(&well_state.wellRates()[0], nw*np, 1); - std::vector qs_dofs(np*nw); - for (int w = 0; w < nw; ++w) { - for (int phase = 0; phase < np; ++phase) { - qs_dofs[w*np + phase] = nc + w*np + phase; - } - } - const V dqs = subset(dx, qs_dofs); - const V qs = qs0 - dqs; - std::copy(&qs[0], &qs[0] + np*nw, well_state.wellRates().begin()); - const V bhp0 = Eigen::Map(&well_state.bhp()[0], nw, 1); - std::vector bhp_dofs(nw); - for (int w = 0; w < nw; ++w) { - bhp_dofs[w] = nc + np*nw + w; - } - ASSERT(bhp_dofs.back() + 1 == total_residual_.size()); - const V dbhp = subset(dx, bhp_dofs); - const V bhp = bhp0 - dbhp; - std::copy(&bhp[0], &bhp[0] + nw, well_state.bhp().begin()); - } - - - - - - double - residualNorm() const - { - return total_residual_.value().matrix().norm(); - } - - - - - - void - computeFluxes(BlackoilState& state) const - { - const int nc = grid_.number_of_cells; - const int np = state.numPhases(); - - const std::vector cells = buildAllCells(nc); - - const V p0 = Eigen::Map(&state.pressure()[0], nc, 1); - const ADB p = ADB::constant(p0, cell_residual_.blockPattern()); - - const V transi = subset(geo_.transmissibility(), - ops_.internal_faces); - const V nkgradp = transi * (ops_.ngrad * p0.matrix()).array(); - - V flux = V::Zero(ops_.internal_faces.size(), 1); - - for (int phase = 0; phase < np; ++phase) { - const ADB cell_rho = fluidRho(phase, p, cells); - const V head = nkgradp + (grav_ * cell_rho.value().matrix()).array(); - const UpwindSelector upwind(grid_, ops_, head); - const V kr = fluidKr(phase); - const V mu = fluidMu(phase, p.value(), cells); - const V mf = upwind.select(kr / mu); - - flux += mf * head; - } - V all_flux = superset(flux, ops_.internal_faces, grid_.number_of_faces); - std::copy(all_flux.data(), all_flux.data() + grid_.number_of_faces, state.faceflux().data()); - } - - - - - - V fluidMu(const int phase, const V& p, const std::vector& cells) const - { - switch (phase) { - case Water: - return fluid_.muWat(p, cells); - case Oil: { - V dummy_rs = V::Zero(p.size(), 1) * p; - return fluid_.muOil(p, dummy_rs, cells); - } - case Gas: - return fluid_.muGas(p, cells); - default: - THROW("Unknown phase index " << phase); - } - } - - - - - - ADB fluidMu(const int phase, const ADB& p, const std::vector& cells) const - { - switch (phase) { - case Water: - return fluid_.muWat(p, cells); - case Oil: { - ADB dummy_rs = V::Zero(p.size(), 1) * p; - return fluid_.muOil(p, dummy_rs, cells); - } - case Gas: - return fluid_.muGas(p, cells); - default: - THROW("Unknown phase index " << phase); - } - } - - - - - - ADB fluidFvf(const int phase, const ADB& p, const std::vector& cells) const - { - switch (phase) { - case Water: - return fluid_.bWat(p, cells); - case Oil: { - ADB dummy_rs = V::Zero(p.size(), 1) * p; - return fluid_.bOil(p, dummy_rs, cells); - } - case Gas: - return fluid_.bGas(p, cells); - default: - THROW("Unknown phase index " << phase); - } - } - - - - - - ADB fluidRho(const int phase, const ADB& p, const std::vector& cells) const - { - const double* rhos = fluid_.surfaceDensity(); - ADB b = fluidFvf(phase, p, cells); - ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b; - return rho; - } - - - - - - V fluidKr(const int phase) const - { - return kr_[phase]; - } - - - - - - V fluidKrWell(const int phase) const - { - return well_kr_[phase]; - } - + // Methods for assembling and solving. + void assemble(const double dt, + const BlackoilState& state, + const WellState& well_state); + void solveJacobianSystem(BlackoilState& state, + WellState& well_state) const; + double residualNorm() const; + void computeFluxes(BlackoilState& state) const; + + // Fluid interface forwarding calls to correct methods of fluid_. + V fluidMu(const int phase, const V& p, const std::vector& cells) const; + ADB fluidMu(const int phase, const ADB& p, const std::vector& cells) const; + ADB fluidFvf(const int phase, const ADB& p, const std::vector& cells) const; + ADB fluidRho(const int phase, const ADB& p, const std::vector& cells) const; + V fluidKr(const int phase) const; + V fluidKrWell(const int phase) const; }; + + } // namespace Opm #endif /* OPM_IMPESTPFAAD_HEADER_INCLUDED */