mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
Merge remote-tracking branch 'atgeirr/fully-implicit' into fully-implicit
Conflicts: opm/autodiff/ImplicitBOStep.hpp
This commit is contained in:
Bård Skaflestad
commit
cf1885fb04
@ -28,6 +28,7 @@
|
||||
list (APPEND MAIN_SOURCE_FILES
|
||||
opm/autodiff/BlackoilPropsAd.cpp
|
||||
opm/autodiff/BlackoilPropsAdInterface.cpp
|
||||
opm/autodiff/FullyImplicitBlackoilSolver.cpp
|
||||
opm/autodiff/ImpesTPFAAD.cpp
|
||||
opm/autodiff/SimulatorCompressibleAd.cpp
|
||||
opm/autodiff/SimulatorIncompTwophaseAdfi.cpp
|
||||
@ -73,7 +74,7 @@ list (APPEND PUBLIC_HEADER_FILES
|
||||
opm/autodiff/BlackoilPropsAdInterface.hpp
|
||||
opm/autodiff/GeoProps.hpp
|
||||
opm/autodiff/ImpesTPFAAD.hpp
|
||||
opm/autodiff/ImplicitBOStep.hpp
|
||||
opm/autodiff/FullyImplicitBlackoilSolver.hpp
|
||||
opm/autodiff/SimulatorCompressibleAd.hpp
|
||||
opm/autodiff/SimulatorIncompTwophaseAdfi.hpp
|
||||
opm/autodiff/TransportSolverTwophaseAd.hpp
|
||||
|
||||
@ -20,7 +20,7 @@
|
||||
|
||||
#include <config.h>
|
||||
|
||||
#include <opm/autodiff/ImplicitBOStep.hpp>
|
||||
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
|
||||
#include <opm/autodiff/GeoProps.hpp>
|
||||
#include <opm/autodiff/BlackoilPropsAd.hpp>
|
||||
|
||||
@ -50,18 +50,18 @@ main(int argc, char* argv[])
|
||||
const Opm::BlackoilPropertiesBasic props0(param, 2, nc);
|
||||
const Opm::BlackoilPropsAd props(props0);
|
||||
|
||||
typedef Opm::ImplicitBOStep<Opm::DerivedGeology> IStep;
|
||||
typedef Opm::FullyImplicitBlackoilSolver BOSolver;
|
||||
|
||||
double grav[] = { 1.0, 0.0 };
|
||||
Opm::DerivedGeology geo(*g, props, grav);
|
||||
|
||||
IStep is(*g, props, geo);
|
||||
BOSolver solver(*g, props, geo);
|
||||
|
||||
Opm::BlackoilState state;
|
||||
initStateBasic(*g, props0, param, 0.0, state);
|
||||
initBlackoilSurfvol(*g, props0, state);
|
||||
|
||||
is.step(1.0, state);
|
||||
solver.step(1.0, state);
|
||||
|
||||
#if 0
|
||||
Opm::WellState well_state;
|
||||
|
||||
@ -171,6 +171,9 @@ namespace Opm
|
||||
ADB BlackoilPropsAd::muWat(const ADB& pw,
|
||||
const Cells& cells) const
|
||||
{
|
||||
#if 1
|
||||
return ADB::constant(muWat(pw.value(), cells), pw.blockPattern());
|
||||
#else
|
||||
if (!pu_.phase_used[Water]) {
|
||||
THROW("Cannot call muWat(): water phase not present.");
|
||||
}
|
||||
@ -188,6 +191,7 @@ namespace Opm
|
||||
jacs[block] = dmu_diag * pw.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu.col(pu_.phase_pos[Water]), jacs);
|
||||
#endif
|
||||
}
|
||||
|
||||
/// Oil viscosity.
|
||||
@ -199,6 +203,9 @@ namespace Opm
|
||||
const ADB& rs,
|
||||
const Cells& cells) const
|
||||
{
|
||||
#if 1
|
||||
return ADB::constant(muOil(po.value(), rs.value(), cells), po.blockPattern());
|
||||
#else
|
||||
if (!pu_.phase_used[Oil]) {
|
||||
THROW("Cannot call muOil(): oil phase not present.");
|
||||
}
|
||||
@ -225,6 +232,7 @@ namespace Opm
|
||||
jacs[block] = dmu_diag * po.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu.col(pu_.phase_pos[Oil]), jacs);
|
||||
#endif
|
||||
}
|
||||
|
||||
/// Gas viscosity.
|
||||
@ -234,6 +242,9 @@ namespace Opm
|
||||
ADB BlackoilPropsAd::muGas(const ADB& pg,
|
||||
const Cells& cells) const
|
||||
{
|
||||
#if 1
|
||||
return ADB::constant(muGas(pg.value(), cells), pg.blockPattern());
|
||||
#else
|
||||
if (!pu_.phase_used[Gas]) {
|
||||
THROW("Cannot call muGas(): gas phase not present.");
|
||||
}
|
||||
@ -251,6 +262,7 @@ namespace Opm
|
||||
jacs[block] = dmu_diag * pg.derivative()[block];
|
||||
}
|
||||
return ADB::function(mu.col(pu_.phase_pos[Gas]), jacs);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
|
||||
483
opm/autodiff/FullyImplicitBlackoilSolver.cpp
Normal file
483
opm/autodiff/FullyImplicitBlackoilSolver.cpp
Normal file
@ -0,0 +1,483 @@
|
||||
/*
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
|
||||
|
||||
|
||||
#include <opm/autodiff/AutoDiffBlock.hpp>
|
||||
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
||||
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
||||
#include <opm/autodiff/GeoProps.hpp>
|
||||
|
||||
#include <opm/core/simulator/BlackoilState.hpp>
|
||||
#include <opm/core/grid.h>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
|
||||
typedef AutoDiff::ForwardBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef ADB::M M;
|
||||
typedef Eigen::Array<double,
|
||||
Eigen::Dynamic,
|
||||
Eigen::Dynamic,
|
||||
Eigen::RowMajor> DataBlock;
|
||||
|
||||
|
||||
namespace {
|
||||
std::vector<int>
|
||||
buildAllCells(const int nc)
|
||||
{
|
||||
std::vector<int> all_cells(nc);
|
||||
|
||||
for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
|
||||
|
||||
return all_cells;
|
||||
}
|
||||
|
||||
template <class GeoProps>
|
||||
AutoDiff::ForwardBlock<double>::M
|
||||
gravityOperator(const UnstructuredGrid& grid,
|
||||
const HelperOps& ops ,
|
||||
const GeoProps& geo )
|
||||
{
|
||||
const int nc = grid.number_of_cells;
|
||||
|
||||
std::vector<int> 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<double>::V V;
|
||||
typedef AutoDiff::ForwardBlock<double>::M M;
|
||||
|
||||
const V& gpot = geo.gravityPotential();
|
||||
const V& trans = geo.transmissibility();
|
||||
|
||||
const HelperOps::IFaces::Index ni = ops.internal_faces.size();
|
||||
|
||||
typedef Eigen::Triplet<double> Tri;
|
||||
std::vector<Tri> 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;
|
||||
}
|
||||
|
||||
template <class PU>
|
||||
std::vector<bool>
|
||||
activePhases(const PU& pu)
|
||||
{
|
||||
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
|
||||
std::vector<bool> active(maxnp, false);
|
||||
|
||||
for (int p = 0; p < pu.MaxNumPhases; ++p) {
|
||||
active[ p ] = pu.phase_used[ p ] != 0;
|
||||
}
|
||||
|
||||
return active;
|
||||
}
|
||||
|
||||
template <class PU>
|
||||
std::vector<int>
|
||||
active2Canonical(const PU& pu)
|
||||
{
|
||||
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
|
||||
std::vector<int> act2can(maxnp, -1);
|
||||
|
||||
for (int phase = 0; phase < maxnp; ++phase) {
|
||||
if (pu.phase_used[ phase ]) {
|
||||
act2can[ pu.phase_pos[ phase ] ] = phase;
|
||||
}
|
||||
}
|
||||
|
||||
return act2can;
|
||||
}
|
||||
} // Anonymous namespace
|
||||
|
||||
|
||||
namespace Opm {
|
||||
|
||||
|
||||
FullyImplicitBlackoilSolver::FullyImplicitBlackoilSolver(const UnstructuredGrid& grid ,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
const DerivedGeology& geo )
|
||||
: grid_ (grid)
|
||||
, fluid_ (fluid)
|
||||
, geo_ (geo)
|
||||
, active_(activePhases(fluid.phaseUsage()))
|
||||
, canph_ (active2Canonical(fluid.phaseUsage()))
|
||||
, cells_ (buildAllCells(grid.number_of_cells))
|
||||
, ops_ (grid)
|
||||
, grav_ (gravityOperator(grid_, ops_, geo_))
|
||||
, rq_ (fluid.numPhases())
|
||||
{
|
||||
allocateResidual();
|
||||
}
|
||||
|
||||
void
|
||||
FullyImplicitBlackoilSolver::step(const double dt,
|
||||
BlackoilState& x)
|
||||
{
|
||||
const V dtpv = geo_.poreVolume() / dt;
|
||||
|
||||
{
|
||||
const SolutionState state = constantState(x);
|
||||
computeAccum(state, 0);
|
||||
}
|
||||
|
||||
#if 0
|
||||
const double atol = 1.0e-15;
|
||||
const double rtol = 5.0e-10;
|
||||
const int maxit = 15;
|
||||
#endif
|
||||
|
||||
assemble(dtpv, x);
|
||||
|
||||
#if 0
|
||||
const double r0 = residualNorm();
|
||||
int it = 0;
|
||||
bool resTooLarge = r0 > atol;
|
||||
while (resTooLarge && (it < maxit)) {
|
||||
solveJacobianSystem(x);
|
||||
|
||||
assemble(dtpv, x);
|
||||
|
||||
const double r = residualNorm();
|
||||
|
||||
resTooLarge = (r > atol) && (r > rtol*r0);
|
||||
|
||||
it += 1;
|
||||
}
|
||||
|
||||
if (resTooLarge) {
|
||||
THROW("Failed to compute converge solution");
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
FullyImplicitBlackoilSolver::ReservoirResidualQuant::ReservoirResidualQuant()
|
||||
: accum(2, ADB::null())
|
||||
, mflux( ADB::null())
|
||||
, b ( ADB::null())
|
||||
, head ( ADB::null())
|
||||
, mob ( ADB::null())
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
FullyImplicitBlackoilSolver::SolutionState::SolutionState(const int np)
|
||||
: pressure ( ADB::null())
|
||||
, saturation(np, ADB::null())
|
||||
, Rs ( ADB::null())
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
void
|
||||
FullyImplicitBlackoilSolver::allocateResidual()
|
||||
{
|
||||
residual_.reservoir.resize(fluid_.numPhases(), ADB::null());
|
||||
}
|
||||
|
||||
FullyImplicitBlackoilSolver::SolutionState
|
||||
FullyImplicitBlackoilSolver::constantState(const BlackoilState& x)
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const int np = x.numPhases();
|
||||
|
||||
const std::vector<int> bpat(np, nc);
|
||||
|
||||
SolutionState state(np);
|
||||
|
||||
assert (not x.pressure().empty());
|
||||
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
|
||||
state.pressure = ADB::constant(p, bpat);
|
||||
|
||||
assert (not x.saturation().empty());
|
||||
const DataBlock s =
|
||||
Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
|
||||
const Opm::PhaseUsage pu = fluid_.phaseUsage();
|
||||
|
||||
{
|
||||
V so = V::Ones(nc, 1);
|
||||
if (active_[ Water ]) {
|
||||
const int pos = pu.phase_pos[ Water ];
|
||||
const V sw = s.col(pos);
|
||||
so -= sw;
|
||||
|
||||
state.saturation[pos] = ADB::constant(sw, bpat);
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
const int pos = pu.phase_pos[ Gas ];
|
||||
const V sg = s.col(pos);
|
||||
so -= sg;
|
||||
|
||||
state.saturation[pos] = ADB::constant(sg, bpat);
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
const int pos = pu.phase_pos[ Oil ];
|
||||
state.saturation[pos] = ADB::constant(so, bpat);
|
||||
}
|
||||
}
|
||||
|
||||
// Ignore miscibility effects (no dissolved gas) for now!
|
||||
const V Rs = V::Zero(nc, 1);
|
||||
state.Rs = ADB::constant(Rs, bpat);
|
||||
|
||||
return state;
|
||||
}
|
||||
|
||||
FullyImplicitBlackoilSolver::SolutionState
|
||||
FullyImplicitBlackoilSolver::variableState(const BlackoilState& x)
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const int np = x.numPhases();
|
||||
|
||||
std::vector<V> vars0;
|
||||
vars0.reserve(np);
|
||||
|
||||
assert (not x.pressure().empty());
|
||||
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
|
||||
vars0.push_back(p);
|
||||
|
||||
assert (not x.saturation().empty());
|
||||
const DataBlock s =
|
||||
Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
|
||||
const Opm::PhaseUsage pu = fluid_.phaseUsage();
|
||||
|
||||
if (active_[ Water ]) {
|
||||
const V sw = s.col(pu.phase_pos[ Water ]);
|
||||
vars0.push_back(sw);
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
const V sg = s.col(pu.phase_pos[ Gas ]);
|
||||
vars0.push_back(sg);
|
||||
}
|
||||
|
||||
std::vector<ADB> vars = ADB::variables(vars0);
|
||||
|
||||
SolutionState state(np);
|
||||
state.pressure = vars[0];
|
||||
|
||||
const std::vector<int>& bpat = vars[0].blockPattern();
|
||||
{
|
||||
ADB so = ADB::constant(V::Ones(nc, 1), bpat);
|
||||
int off = 1; // First saturation variable at offset 1.
|
||||
|
||||
if (active_[ Water ]) {
|
||||
ADB& sw = vars[ off++ ];
|
||||
state.saturation[ pu.phase_pos[ Water ] ] = sw;
|
||||
|
||||
so = so - sw;
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
ADB& sg = vars[ off++ ];
|
||||
state.saturation[ pu.phase_pos[ Gas ] ] = sg;
|
||||
|
||||
so = so - sg;
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
state.saturation[ pu.phase_pos[ Oil ] ] = so;
|
||||
}
|
||||
}
|
||||
|
||||
// Ignore miscibility effects (no dissolved gas) for now!
|
||||
V Rs = V::Zero(nc, 1);
|
||||
state.Rs = ADB::constant(Rs, bpat);
|
||||
|
||||
return state;
|
||||
}
|
||||
|
||||
void
|
||||
FullyImplicitBlackoilSolver::computeAccum(const SolutionState& state,
|
||||
const int aix )
|
||||
{
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
|
||||
const ADB& press = state.pressure;
|
||||
const std::vector<ADB>& sat = state.saturation;
|
||||
|
||||
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
|
||||
for (int phase = 0; phase < maxnp; ++phase) {
|
||||
if (active_[ phase ]) {
|
||||
const int pos = pu.phase_pos[ phase ];
|
||||
rq_[pos].b = fluidReciprocFVF(phase, press, cells_);
|
||||
rq_[pos].accum[aix] = rq_[pos].b * sat[pos];
|
||||
}
|
||||
}
|
||||
|
||||
if (active_[ Oil ] && active_[ Gas ]) {
|
||||
// Account for gas dissolved in oil.
|
||||
const int po = pu.phase_pos[ Oil ];
|
||||
const int pg = pu.phase_pos[ Gas ];
|
||||
|
||||
rq_[pg].accum[aix] += state.Rs * rq_[po].accum[aix];
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
FullyImplicitBlackoilSolver::assemble(const V& dtpv, const BlackoilState& x)
|
||||
{
|
||||
const V transi = subset(geo_.transmissibility(),
|
||||
ops_.internal_faces);
|
||||
|
||||
const SolutionState state = variableState(x);
|
||||
const std::vector<ADB> kr = computeRelPerm(state);
|
||||
|
||||
computeAccum(state, 1);
|
||||
|
||||
for (int phase = 0; phase < fluid_.numPhases(); ++phase) {
|
||||
computeMassFlux(phase, transi, kr, state);
|
||||
|
||||
residual_.reservoir[ phase ] =
|
||||
dtpv*(rq_[phase].accum[1] - rq_[phase].accum[0])
|
||||
+ ops_.div*rq_[phase].mflux;
|
||||
}
|
||||
|
||||
if (active_[ Oil ] && active_[ Gas ]) {
|
||||
const int po = fluid_.phaseUsage().phase_pos[ Oil ];
|
||||
const UpwindSelector<double> upwind(grid_, ops_,
|
||||
rq_[po].head.value());
|
||||
const ADB Rs = upwind.select(state.Rs);
|
||||
|
||||
residual_.reservoir[ Gas ] += ops_.div * (Rs * rq_[po].mflux);
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<ADB>
|
||||
FullyImplicitBlackoilSolver::computeRelPerm(const SolutionState& state)
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const std::vector<int>& bpat = state.pressure.blockPattern();
|
||||
|
||||
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
||||
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
const ADB sw = (active_[ Water ]
|
||||
? state.saturation[ pu.phase_pos[ Water ] ]
|
||||
: null);
|
||||
|
||||
const ADB so = (active_[ Oil ]
|
||||
? state.saturation[ pu.phase_pos[ Oil ] ]
|
||||
: null);
|
||||
|
||||
const ADB sg = (active_[ Gas ]
|
||||
? state.saturation[ pu.phase_pos[ Gas ] ]
|
||||
: null);
|
||||
|
||||
return fluid_.relperm(sw, so, sg, cells_);
|
||||
}
|
||||
|
||||
void
|
||||
FullyImplicitBlackoilSolver::computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const std::vector<ADB>& kr ,
|
||||
const SolutionState& state )
|
||||
{
|
||||
const int phase = canph_[ actph ];
|
||||
const ADB mu = fluidViscosity(phase, state.pressure, cells_);
|
||||
|
||||
rq_[ actph ].mob = kr[ phase ] / mu;
|
||||
|
||||
const ADB rho = fluidDensity(phase, state.pressure, cells_);
|
||||
const ADB gflux = grav_ * rho;
|
||||
|
||||
ADB& head = rq_[ actph ].head;
|
||||
head = transi*(ops_.ngrad * state.pressure) + gflux;
|
||||
|
||||
UpwindSelector<double> upwind(grid_, ops_, head.value());
|
||||
|
||||
const ADB& b = rq_[ actph ].b;
|
||||
const ADB& mob = rq_[ actph ].mob;
|
||||
rq_[ actph ].mflux = upwind.select(b * mob) * head;
|
||||
}
|
||||
|
||||
ADB
|
||||
FullyImplicitBlackoilSolver::fluidViscosity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& 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
|
||||
FullyImplicitBlackoilSolver::fluidReciprocFVF(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& 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
|
||||
FullyImplicitBlackoilSolver::fluidDensity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const
|
||||
{
|
||||
const double* rhos = fluid_.surfaceDensity();
|
||||
ADB b = fluidReciprocFVF(phase, p, cells);
|
||||
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
|
||||
return rho;
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
147
opm/autodiff/FullyImplicitBlackoilSolver.hpp
Normal file
147
opm/autodiff/FullyImplicitBlackoilSolver.hpp
Normal file
@ -0,0 +1,147 @@
|
||||
/*
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
|
||||
#define OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
|
||||
|
||||
#include <opm/autodiff/AutoDiffBlock.hpp>
|
||||
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
||||
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
||||
|
||||
struct UnstructuredGrid;
|
||||
struct Wells;
|
||||
|
||||
namespace Opm {
|
||||
|
||||
class DerivedGeology;
|
||||
class LinearSolverInterface;
|
||||
class BlackoilState;
|
||||
class WellState;
|
||||
|
||||
|
||||
/// A fully implicit TPFA-based solver for the black-oil problem.
|
||||
class FullyImplicitBlackoilSolver
|
||||
{
|
||||
public:
|
||||
FullyImplicitBlackoilSolver(const UnstructuredGrid& grid ,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
const DerivedGeology& geo );
|
||||
|
||||
/// Take a single forward step, modifiying
|
||||
/// state.pressure()
|
||||
/// state.faceflux()
|
||||
/// state.saturation()
|
||||
/// state.surfacevol()
|
||||
void
|
||||
step(const double dt,
|
||||
BlackoilState& state);
|
||||
|
||||
private:
|
||||
// Types and enums
|
||||
typedef AutoDiff::ForwardBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef ADB::M M;
|
||||
typedef Eigen::Array<double,
|
||||
Eigen::Dynamic,
|
||||
Eigen::Dynamic,
|
||||
Eigen::RowMajor> DataBlock;
|
||||
|
||||
struct ReservoirResidualQuant {
|
||||
ReservoirResidualQuant();
|
||||
std::vector<ADB> accum; // Accumulations
|
||||
ADB mflux; // Mass flux (surface conditions)
|
||||
ADB b; // Reciprocal FVF
|
||||
ADB head; // Pressure drop across int. interfaces
|
||||
ADB mob; // Phase mobility (per cell)
|
||||
};
|
||||
|
||||
struct SolutionState {
|
||||
SolutionState(const int np);
|
||||
ADB pressure;
|
||||
std::vector<ADB> saturation;
|
||||
ADB Rs;
|
||||
};
|
||||
|
||||
enum { Water = BlackoilPropsAdInterface::Water,
|
||||
Oil = BlackoilPropsAdInterface::Oil ,
|
||||
Gas = BlackoilPropsAdInterface::Gas };
|
||||
|
||||
// Member data
|
||||
const UnstructuredGrid& grid_;
|
||||
const BlackoilPropsAdInterface& fluid_;
|
||||
const DerivedGeology& geo_;
|
||||
// For each canonical phase -> true if active
|
||||
const std::vector<bool> active_;
|
||||
// Size = # active faces. Maps active -> canonical phase indices.
|
||||
const std::vector<int> canph_;
|
||||
const std::vector<int> cells_; // All grid cells
|
||||
HelperOps ops_;
|
||||
const M grav_;
|
||||
|
||||
std::vector<ReservoirResidualQuant> rq_;
|
||||
|
||||
struct {
|
||||
std::vector<ADB> reservoir;
|
||||
} residual_;
|
||||
|
||||
// Private methods.
|
||||
void
|
||||
allocateResidual();
|
||||
|
||||
SolutionState
|
||||
constantState(const BlackoilState& x);
|
||||
|
||||
SolutionState
|
||||
variableState(const BlackoilState& x);
|
||||
|
||||
void
|
||||
computeAccum(const SolutionState& state,
|
||||
const int aix );
|
||||
|
||||
void
|
||||
assemble(const V& dtpv, const BlackoilState& x);
|
||||
|
||||
std::vector<ADB>
|
||||
computeRelPerm(const SolutionState& state);
|
||||
|
||||
void
|
||||
computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const std::vector<ADB>& kr ,
|
||||
const SolutionState& state );
|
||||
|
||||
ADB
|
||||
fluidViscosity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const;
|
||||
|
||||
ADB
|
||||
fluidReciprocFVF(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const;
|
||||
|
||||
ADB
|
||||
fluidDensity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const;
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
#endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
|
||||
@ -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.
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
@ -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]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -1,535 +0,0 @@
|
||||
/*
|
||||
Copyright 2013 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2013 Statoil ASA.
|
||||
|
||||
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 <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_IMPLICITBO_HPP_HEADER_INCLUDED
|
||||
#define OPM_IMPLICITBO_HPP_HEADER_INCLUDED
|
||||
|
||||
#include <opm/autodiff/AutoDiffBlock.hpp>
|
||||
#include <opm/autodiff/AutoDiffHelpers.hpp>
|
||||
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
|
||||
|
||||
#include <opm/core/grid.h>
|
||||
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <cmath>
|
||||
|
||||
namespace {
|
||||
std::vector<int>
|
||||
buildAllCells(const int nc)
|
||||
{
|
||||
std::vector<int> all_cells(nc);
|
||||
|
||||
for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
|
||||
|
||||
return all_cells;
|
||||
}
|
||||
|
||||
template <class GeoProps>
|
||||
AutoDiff::ForwardBlock<double>::M
|
||||
gravityOperator(const UnstructuredGrid& grid,
|
||||
const HelperOps& ops ,
|
||||
const GeoProps& geo )
|
||||
{
|
||||
const int nc = grid.number_of_cells;
|
||||
|
||||
std::vector<int> 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<double>::V V;
|
||||
typedef AutoDiff::ForwardBlock<double>::M M;
|
||||
|
||||
const V& gpot = geo.gravityPotential();
|
||||
const V& trans = geo.transmissibility();
|
||||
|
||||
const HelperOps::IFaces::Index ni = ops.internal_faces.size();
|
||||
|
||||
typedef Eigen::Triplet<double> Tri;
|
||||
std::vector<Tri> 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;
|
||||
}
|
||||
|
||||
template <class PU>
|
||||
std::vector<bool>
|
||||
activePhases(const PU& pu)
|
||||
{
|
||||
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
|
||||
std::vector<bool> active(maxnp, false);
|
||||
|
||||
for (int p = 0; p < pu.MaxNumPhases; ++p) {
|
||||
active[ p ] = pu.phase_used[ p ] != 0;
|
||||
}
|
||||
|
||||
return active;
|
||||
}
|
||||
|
||||
template <class PU>
|
||||
std::vector<int>
|
||||
active2Canonical(const PU& pu)
|
||||
{
|
||||
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
|
||||
std::vector<int> act2can(maxnp, -1);
|
||||
|
||||
for (int phase = 0; phase < maxnp; ++phase) {
|
||||
if (pu.phase_used[ phase ]) {
|
||||
act2can[ pu.phase_pos[ phase ] ] = phase;
|
||||
}
|
||||
}
|
||||
|
||||
return act2can;
|
||||
}
|
||||
} // Anonymous namespace
|
||||
|
||||
|
||||
namespace Opm {
|
||||
|
||||
template <class GeoProps>
|
||||
class ImplicitBOStep {
|
||||
public:
|
||||
ImplicitBOStep(const UnstructuredGrid& grid ,
|
||||
const BlackoilPropsAdInterface& fluid,
|
||||
const GeoProps& geo )
|
||||
: grid_ (grid)
|
||||
, fluid_ (fluid)
|
||||
, geo_ (geo)
|
||||
, active_(activePhases(fluid.phaseUsage()))
|
||||
, canph_ (active2Canonical(fluid.phaseUsage()))
|
||||
, cells_ (buildAllCells(grid.number_of_cells))
|
||||
, ops_ (grid)
|
||||
, grav_ (gravityOperator(grid_, ops_, geo_))
|
||||
, rq_ (fluid.numPhases())
|
||||
{
|
||||
allocateResidual();
|
||||
}
|
||||
|
||||
template <class ResSol>
|
||||
void
|
||||
step(const double dt,
|
||||
ResSol& x)
|
||||
{
|
||||
const V dtpv = geo_.poreVolume() / dt;
|
||||
|
||||
{
|
||||
const SolutionState state = constantState(x);
|
||||
computeAccum(state, 0);
|
||||
}
|
||||
|
||||
const double atol = 1.0e-15;
|
||||
#if 0
|
||||
const double rtol = 5.0e-10;
|
||||
const int maxit = 15;
|
||||
#endif
|
||||
|
||||
assemble(dtpv, x);
|
||||
|
||||
const double r0 = residualNorm();
|
||||
bool resTooLarge = r0 > atol;
|
||||
#if 0
|
||||
int it = 0;
|
||||
while (resTooLarge && (it < maxit)) {
|
||||
solveJacobianSystem(x);
|
||||
|
||||
assemble(dtpv, x);
|
||||
|
||||
const double r = residualNorm();
|
||||
|
||||
resTooLarge = (r > atol) && (r > rtol*r0);
|
||||
|
||||
it += 1;
|
||||
}
|
||||
#endif
|
||||
if (resTooLarge) {
|
||||
THROW("Failed to compute converge solution");
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
typedef AutoDiff::ForwardBlock<double> ADB;
|
||||
typedef ADB::V V;
|
||||
typedef ADB::M M;
|
||||
typedef Eigen::Array<double,
|
||||
Eigen::Dynamic,
|
||||
Eigen::Dynamic,
|
||||
Eigen::RowMajor> DataBlock;
|
||||
|
||||
struct ReservoirResidualQuant {
|
||||
ReservoirResidualQuant()
|
||||
: accum(2, ADB::null())
|
||||
, mflux( ADB::null())
|
||||
, b ( ADB::null())
|
||||
, head ( ADB::null())
|
||||
, mob ( ADB::null())
|
||||
{
|
||||
}
|
||||
|
||||
std::vector<ADB> accum; // Accumulations
|
||||
ADB mflux; // Mass flux (surface conditions)
|
||||
ADB b; // Reciprocal FVF
|
||||
ADB head; // Pressure drop across int. interfaces
|
||||
ADB mob; // Phase mobility (per cell)
|
||||
};
|
||||
|
||||
struct SolutionState {
|
||||
SolutionState(const int np)
|
||||
: pressure ( ADB::null())
|
||||
, saturation(np, ADB::null())
|
||||
, Rs ( ADB::null())
|
||||
{
|
||||
}
|
||||
|
||||
ADB pressure;
|
||||
std::vector<ADB> saturation;
|
||||
ADB Rs;
|
||||
};
|
||||
|
||||
const UnstructuredGrid& grid_;
|
||||
const BlackoilPropsAdInterface& fluid_;
|
||||
const GeoProps& geo_;
|
||||
const std::vector<bool> active_; // Canonical -> active
|
||||
const std::vector<int> canph_; // Active -> canonical phases
|
||||
const std::vector<int> cells_; // All grid cells
|
||||
HelperOps ops_;
|
||||
const M grav_;
|
||||
|
||||
std::vector<ReservoirResidualQuant> rq_;
|
||||
|
||||
struct {
|
||||
std::vector<ADB> reservoir;
|
||||
} residual_;
|
||||
|
||||
enum { Water = BlackoilPropsAdInterface::Water,
|
||||
Oil = BlackoilPropsAdInterface::Oil ,
|
||||
Gas = BlackoilPropsAdInterface::Gas };
|
||||
|
||||
void
|
||||
allocateResidual()
|
||||
{
|
||||
residual_.reservoir.resize(fluid_.numPhases(), ADB::null());
|
||||
}
|
||||
|
||||
template <class ResSol>
|
||||
SolutionState
|
||||
constantState(const ResSol& x)
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const int np = x.numPhases();
|
||||
|
||||
const std::vector<int> bpat(np, nc);
|
||||
|
||||
SolutionState state(np);
|
||||
|
||||
assert (not x.pressure().empty());
|
||||
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
|
||||
state.pressure = ADB::constant(p, bpat);
|
||||
|
||||
assert (not x.saturation().empty());
|
||||
const DataBlock s =
|
||||
Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
|
||||
const Opm::PhaseUsage pu = fluid_.phaseUsage();
|
||||
|
||||
{
|
||||
V so = V::Ones(nc, 1);
|
||||
if (active_[ Water ]) {
|
||||
const int pos = pu.phase_pos[ Water ];
|
||||
const V sw = s.col(pos);
|
||||
so -= sw;
|
||||
|
||||
state.saturation[pos] = ADB::constant(sw, bpat);
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
const int pos = pu.phase_pos[ Gas ];
|
||||
const V sg = s.col(pos);
|
||||
so -= sg;
|
||||
|
||||
state.saturation[pos] = ADB::constant(sg, bpat);
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
const int pos = pu.phase_pos[ Oil ];
|
||||
state.saturation[pos] = ADB::constant(so, bpat);
|
||||
}
|
||||
}
|
||||
|
||||
// Ignore miscibility effects (no dissolved gas) for now!
|
||||
const V Rs = V::Zero(nc, 1);
|
||||
state.Rs = ADB::constant(Rs, bpat);
|
||||
|
||||
return state;
|
||||
}
|
||||
|
||||
template <class ResSol>
|
||||
SolutionState
|
||||
variableState(const ResSol& x)
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const int np = x.numPhases();
|
||||
|
||||
std::vector<V> vars0;
|
||||
vars0.reserve(np);
|
||||
|
||||
assert (not x.pressure().empty());
|
||||
const V p = Eigen::Map<const V>(& x.pressure()[0], nc, 1);
|
||||
vars0.push_back(p);
|
||||
|
||||
assert (not x.saturation().empty());
|
||||
const DataBlock s =
|
||||
Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
|
||||
|
||||
const Opm::PhaseUsage pu = fluid_.phaseUsage();
|
||||
|
||||
if (active_[ Water ]) {
|
||||
const V sw = s.col(pu.phase_pos[ Water ]);
|
||||
vars0.push_back(sw);
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
const V sg = s.col(pu.phase_pos[ Gas ]);
|
||||
vars0.push_back(sg);
|
||||
}
|
||||
|
||||
std::vector<ADB> vars = ADB::variables(vars0);
|
||||
|
||||
SolutionState state(np);
|
||||
state.pressure = vars[0];
|
||||
|
||||
const std::vector<int>& bpat = vars[0].blockPattern();
|
||||
{
|
||||
ADB so = ADB::constant(V::Ones(nc, 1), bpat);
|
||||
int off = 1; // First saturation variable at offset 1.
|
||||
|
||||
if (active_[ Water ]) {
|
||||
ADB& sw = vars[ off++ ];
|
||||
state.saturation[ pu.phase_pos[ Water ] ] = sw;
|
||||
|
||||
so = so - sw;
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
ADB& sg = vars[ off++ ];
|
||||
state.saturation[ pu.phase_pos[ Gas ] ] = sg;
|
||||
|
||||
so = so - sg;
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
state.saturation[ pu.phase_pos[ Oil ] ] = so;
|
||||
}
|
||||
}
|
||||
|
||||
// Ignore miscibility effects (no dissolved gas) for now!
|
||||
V Rs = V::Zero(nc, 1);
|
||||
state.Rs = ADB::constant(Rs, bpat);
|
||||
|
||||
return state;
|
||||
}
|
||||
|
||||
void
|
||||
computeAccum(const SolutionState& state,
|
||||
const int aix )
|
||||
{
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
|
||||
const ADB& press = state.pressure;
|
||||
const std::vector<ADB>& sat = state.saturation;
|
||||
|
||||
const int maxnp = Opm::BlackoilPhases::MaxNumPhases;
|
||||
for (int phase = 0; phase < maxnp; ++phase) {
|
||||
if (active_[ phase ]) {
|
||||
const int pos = pu.phase_pos[ phase ];
|
||||
rq_[pos].b = fluidReciprocFVF(phase, press, cells_);
|
||||
rq_[pos].accum[aix] = rq_[pos].b * sat[pos];
|
||||
}
|
||||
}
|
||||
|
||||
if (active_[ Oil ] && active_[ Gas ]) {
|
||||
// Account for gas dissolved in oil.
|
||||
const int po = pu.phase_pos[ Oil ];
|
||||
const int pg = pu.phase_pos[ Gas ];
|
||||
|
||||
rq_[pg].accum[aix] += state.Rs * rq_[po].accum[aix];
|
||||
}
|
||||
}
|
||||
|
||||
template <class ResSol>
|
||||
void
|
||||
assemble(const V& dtpv, const ResSol& x)
|
||||
{
|
||||
const V transi = subset(geo_.transmissibility(),
|
||||
ops_.internal_faces);
|
||||
|
||||
const SolutionState state = variableState(x);
|
||||
const std::vector<ADB> kr = computeRelPerm(state);
|
||||
|
||||
computeAccum(state, 1);
|
||||
|
||||
for (int phase = 0; phase < fluid_.numPhases(); ++phase) {
|
||||
computeMassFlux(phase, transi, kr, state);
|
||||
|
||||
residual_.reservoir[ phase ] =
|
||||
dtpv*(rq_[phase].accum[1] - rq_[phase].accum[0])
|
||||
+ ops_.div*rq_[phase].mflux;
|
||||
}
|
||||
|
||||
if (active_[ Oil ] && active_[ Gas ]) {
|
||||
const int po = fluid_.phaseUsage().phase_pos[ Oil ];
|
||||
const UpwindSelector<double> upwind(grid_, ops_,
|
||||
rq_[po].head.value());
|
||||
const ADB Rs = upwind.select(state.Rs);
|
||||
|
||||
residual_.reservoir[ Gas ] += ops_.div * (Rs * rq_[po].mflux);
|
||||
}
|
||||
}
|
||||
|
||||
std::vector<ADB>
|
||||
computeRelPerm(const SolutionState& state)
|
||||
{
|
||||
const int nc = grid_.number_of_cells;
|
||||
const std::vector<int>& bpat = state.pressure.blockPattern();
|
||||
|
||||
const ADB null = ADB::constant(V::Zero(nc, 1), bpat);
|
||||
|
||||
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
|
||||
const ADB sw = (active_[ Water ]
|
||||
? state.saturation[ pu.phase_pos[ Water ] ]
|
||||
: null);
|
||||
|
||||
const ADB so = (active_[ Oil ]
|
||||
? state.saturation[ pu.phase_pos[ Oil ] ]
|
||||
: null);
|
||||
|
||||
const ADB sg = (active_[ Gas ]
|
||||
? state.saturation[ pu.phase_pos[ Gas ] ]
|
||||
: null);
|
||||
|
||||
return fluid_.relperm(sw, so, sg, cells_);
|
||||
}
|
||||
|
||||
void
|
||||
computeMassFlux(const int actph ,
|
||||
const V& transi,
|
||||
const std::vector<ADB>& kr ,
|
||||
const SolutionState& state )
|
||||
{
|
||||
const int phase = canph_[ actph ];
|
||||
const ADB mu = fluidViscosity(phase, state.pressure, cells_);
|
||||
|
||||
rq_[ actph ].mob = kr[ phase ] / mu;
|
||||
|
||||
const ADB rho = fluidDensity(phase, state.pressure, cells_);
|
||||
const ADB gflux = grav_ * rho;
|
||||
|
||||
ADB& head = rq_[ actph ].head;
|
||||
head = transi*(ops_.ngrad * state.pressure) + gflux;
|
||||
|
||||
UpwindSelector<double> upwind(grid_, ops_, head.value());
|
||||
|
||||
const ADB& b = rq_[ actph ].b;
|
||||
const ADB& mob = rq_[ actph ].mob;
|
||||
rq_[ actph ].mflux = upwind.select(b * mob) * head;
|
||||
}
|
||||
|
||||
double
|
||||
residualNorm() const
|
||||
{
|
||||
double r = 0;
|
||||
for (std::vector<ADB>::const_iterator
|
||||
b = residual_.reservoir.begin(),
|
||||
e = residual_.reservoir.end();
|
||||
b != e; ++b)
|
||||
{
|
||||
r = std::max(r, (*b).value().matrix().norm());
|
||||
}
|
||||
|
||||
return r;
|
||||
}
|
||||
|
||||
ADB
|
||||
fluidViscosity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& 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
|
||||
fluidReciprocFVF(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& 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
|
||||
fluidDensity(const int phase,
|
||||
const ADB& p ,
|
||||
const std::vector<int>& cells) const
|
||||
{
|
||||
const double* rhos = fluid_.surfaceDensity();
|
||||
ADB b = fluidReciprocFVF(phase, p, cells);
|
||||
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
|
||||
return rho;
|
||||
}
|
||||
};
|
||||
} // namespace Opm
|
||||
|
||||
#endif /* OPM_IMPLICITBO_HPP_HEADER_INCLUDED */
|
||||
Loading…
Reference in New Issue
Block a user