OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-28 20:13:49 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into fully-implicit

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2013-05-24 10:40:35 +02:00
commit 6f55996441
3 changed files with 138 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
opm/autodiff/GeoProps.hpp
View File

@ -40,7 +40,7 @@ namespace Opm
const double* grav = 0)
: pvol_ (grid.number_of_cells)
, trans_(grid.number_of_faces)
, gpot_ (grid.cell_facepos[ grid.number_of_cells ])
, gpot_ (Vector::Zero(grid.cell_facepos[ grid.number_of_cells ], 1))
{
// Pore volume
const typename Vector::Index nc = grid.number_of_cells;
@ -54,6 +54,8 @@ namespace Opm
tpfa_htrans_compute(ug, props.permeability(), htrans.data());
tpfa_trans_compute (ug, htrans.data() , trans_.data());
// Gravity potential
std::fill(gravity_, gravity_ + 3, 0.0);
if (grav != 0) {
const typename Vector::Index nd = grid.dimensions;
@ -71,17 +73,20 @@ namespace Opm
}
}
}
std::copy(grav, grav + nd, gravity_);
}
}
const Vector& poreVolume() const { return pvol_ ; }
const Vector& transmissibility() const { return trans_; }
const Vector& gravityPotential() const { return gpot_ ; }
const double* gravity() const { return gravity_; }
private:
Vector pvol_ ;
Vector trans_;
Vector gpot_ ;
double gravity_[3]; // Size 3 even if grid is 2-dim.
};
}

170
opm/autodiff/ImpesTPFAAD.cpp
View File

@ -28,6 +28,13 @@
#include <iomanip>
// Repeated from inside ImpesTPFAAD for convenience.
typedef AutoDiff::ForwardBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
namespace {
std::vector<int>
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<double> 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<const DataBlock>(state.saturation().data(), nc, np);
ASSERT(np == 2);
@ -153,7 +190,7 @@ namespace Opm {
const double r0 = residualNorm();
int it = 0;
std::cout << "\nIteration Residual\n"
<< std::setw(9) << it
<< std::setw(9) << it << std::setprecision(9)
<< std::setw(18) << r0 << std::endl;
bool resTooLarge = r0 > atol;
while (resTooLarge && (it < maxit)) {
@ -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::setprecision(9)
<< 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<int> cells = buildAllCells(nc);
// Compute relperms.
DataBlock s = Eigen::Map<const DataBlock>(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<const DataBlock>(comp_frac, 1, np);
}
}
const std::vector<int> 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<const V> 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,
@ -202,20 +297,17 @@ namespace Opm {
wells_.well_cells + nperf);
const V transw = Eigen::Map<const V>(wells_.WI, nperf, 1);
// Initialize AD variables: p (cell pressures), qs (well rates) and bhp (well bhp).
// Initialize AD variables: p (cell pressures) and bhp (well bhp).
const V p0 = Eigen::Map<const V>(&state.pressure()[0], nc, 1);
const V qs0 = Eigen::Map<const V>(&well_state.wellRates()[0], nw*np, 1);
const V bhp0 = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
std::vector<V> vars0 = { p0, qs0, bhp0 };
std::vector<V> vars0 = { p0, bhp0 };
std::vector<ADB> vars = ADB::variables(vars0);
const ADB& p = vars[0];
const ADB& qs = vars[1];
const ADB& bhp = vars[2];
const ADB& bhp = vars[1];
std::vector<int> bpat = p.blockPattern();
// Compute T_ij * (p_i - p_j) and use for upwinding.
// Compute T_ij * (p_i - p_j).
const ADB nkgradp = transi * (ops_.ngrad * p);
const UpwindSelector<double> upwind(grid_, ops_, nkgradp.value());
// Extract variables for perforation cell pressures
// and corresponding perforation well pressures.
@ -231,16 +323,15 @@ 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<double> 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);
qs_ = ADB::constant(V::Zero(nw*np, 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);
@ -250,20 +341,22 @@ namespace Opm {
// since they are not available yet.
const V mu = fluidMu(phase, p.value(), cells);
const V cell_mob = kr / mu;
const ADB head_diff_grav = (grav_ * cell_rho);
const ADB head = nkgradp + (grav_ * cell_rho);
const UpwindSelector<double> upwind(grid_, ops_, head.value());
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 flux = face_mob * head;
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 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;
@ -274,7 +367,7 @@ namespace Opm {
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);
qs_ = qs_ + superset(well_rates, well_flow_res_phase_idx, nw*np);
}
cell_residual_ = cell_residual_ + divcontrib_sum;
// Handling BHP and SURFACE_RATE wells.
@ -297,15 +390,17 @@ namespace Opm {
}
}
const ADB bhp_residual = bhp - bhp_targets;
const ADB rate_residual = rate_distr * qs - rate_targets;
const ADB rate_residual = rate_distr * qs_ - rate_targets;
// Choose bhp residual for positive bhp targets.
Selector<double> 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));
// Build full residual by concatenation of residual arrays and
// jacobian matrices.
total_residual_ = collapseJacs(vertcat(cell_residual_, well_residual_));
// std::cout.precision(16);
// std::cout << total_residual_;
}
@ -318,11 +413,11 @@ namespace Opm {
{
const int nc = grid_.number_of_cells;
const int nw = wells_.number_of_wells;
const int np = state.numPhases();
// const int np = state.numPhases();
Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_residual_.derivative()[0];
V dx(total_residual_.size());
V dx(V::Zero(total_residual_.size()));
Opm::LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(matr.rows(), matr.nonZeros(),
matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
@ -334,20 +429,10 @@ namespace Opm {
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<const V>(&well_state.wellRates()[0], nw*np, 1);
std::vector<int> 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<const V>(&well_state.bhp()[0], nw, 1);
std::vector<int> bhp_dofs(nw);
for (int w = 0; w < nw; ++w) {
bhp_dofs[w] = nc + np*nw + w;
bhp_dofs[w] = nc + w;
}
ASSERT(bhp_dofs.back() + 1 == total_residual_.size());
const V dbhp = subset(dx, bhp_dofs);
@ -406,8 +491,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<double> cell_to_well_selector(nkgradp_well);
@ -438,7 +522,7 @@ namespace Opm {
std::copy(perf_flux.data(), perf_flux.data() + nperf, well_state.perfRates().begin());
std::copy(p_perfwell.data(), p_perfwell.data() + nperf, well_state.perfPress().begin());
std::copy(qs_.value().data(), qs_.value().data() + np*nw, &well_state.wellRates()[0]);
}

5
opm/autodiff/ImpesTPFAAD.hpp
View File

@ -88,8 +88,13 @@ namespace Opm {
ADB total_residual_;
std::vector<V> kr_;
std::vector<V> 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);