OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'atgeirr/spe9-fixes'

Browse Source
This commit is contained in:
Kai Bao
2014-04-28 14:09:06 +02:00
parent 3bdfbfd78c 17400cebaf
commit 820e65e772
17 changed files with 283 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
opm/autodiff/AutoDiffBlock.hpp
View File

@@ -20,8 +20,13 @@
#ifndef OPM_AUTODIFFBLOCK_HEADER_INCLUDED
#define OPM_AUTODIFFBLOCK_HEADER_INCLUDED
#include "disable_warning_pragmas.h"
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include "reenable_warning_pragmas.h"
#include <vector>
#include <cassert>

2
opm/autodiff/BlackoilPropsAdFromDeck.hpp
View File

@@ -34,7 +34,9 @@
#include <boost/scoped_ptr.hpp>
#ifdef HAVE_DUNE_CORNERPOINT
#include "disable_warning_pragmas.h"
#include <dune/grid/CpGrid.hpp>
#include "reenable_warning_pragmas.h"
#endif
namespace Opm

9
opm/autodiff/FullyImplicitBlackoilSolver.hpp
View File

@@ -165,14 +165,15 @@ namespace Opm {
void
addWellControlEq(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw);
const WellStateFullyImplicitBlackoil& xw,
const V& aliveWells);
void
addWellEq(const SolutionState& state,
WellStateFullyImplicitBlackoil& xw);
V& aliveWells);
void updateWellControls(const ADB& bhp,
const ADB& well_phase_flow_rate,
void updateWellControls(ADB& bhp,
ADB& well_phase_flow_rate,
WellStateFullyImplicitBlackoil& xw) const;
void

118
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@@ -49,6 +49,19 @@
<< collapseJacs(foo) << std::endl; \
} while (0)
#define DUMPVAL(foo) \
do { \
std::cout << "==========================================\n" \
<< #foo ":\n" \
<< foo.value() << std::endl; \
} while (0)
#define DISKVAL(foo) \
do { \
std::ofstream os(#foo); \
os.precision(16); \
os << foo.value() << std::endl; \
} while (0)
namespace Opm {
@@ -638,7 +651,7 @@ namespace {
{
using namespace Opm::AutoDiffGrid;
// Create the primary variables.
const SolutionState state = variableState(x, xw);
SolutionState state = variableState(x, xw);
// -------- Mass balance equations --------
@@ -695,8 +708,12 @@ namespace {
}
addWellEq(state, xw); // Note: this can change xw (switching well controls).
addWellControlEq(state, xw);
// Note: updateWellControls() can change all its arguments if
// a well control is switched.
updateWellControls(state.bhp, state.qs, xw);
V aliveWells;
addWellEq(state, aliveWells);
addWellControlEq(state, xw, aliveWells);
}
@@ -705,7 +722,7 @@ namespace {
template <class T>
void FullyImplicitBlackoilSolver<T>::addWellEq(const SolutionState& state,
WellStateFullyImplicitBlackoil& xw)
V& aliveWells)
{
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int np = wells_.number_of_phases;
@@ -722,6 +739,10 @@ namespace {
// and corresponding perforation well pressures.
const ADB p_perfcell = subset(state.pressure, well_cells);
// DUMPVAL(p_perfcell);
// DUMPVAL(state.bhp);
// DUMPVAL(ADB::constant(cdp));
// Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcell - (wops_.w2p * state.bhp + cdp);
@@ -804,12 +825,13 @@ namespace {
wbq[phase] = (isInj * compi.col(pos)) * qt_s - q_ps[phase];
wbqt += wbq[phase];
}
// DUMPVAL(wbqt);
// check for dead wells
V isNotDeadWells = V::Constant(nw,1.0);
aliveWells = V::Constant(nw, 1.0);
for (int w = 0; w < nw; ++w) {
if (wbqt.value()[w] == 0) {
isNotDeadWells[w] = 0;
aliveWells[w] = 0.0;
}
}
// compute wellbore mixture at std conds
@@ -827,9 +849,13 @@ namespace {
ADB mt = subset(rq_[0].mob,well_cells);
for (int phase = 1; phase < np; ++phase) {
mt += subset(rq_[phase].mob,well_cells);
}
// DUMPVAL(ADB::constant(isInjInx));
// DUMPVAL(ADB::constant(Tw));
// DUMPVAL(mt);
// DUMPVAL(drawdown);
// injection connections total volumerates
ADB cqt_i = -(isInjInx * Tw) * (mt * drawdown);
@@ -856,9 +882,11 @@ namespace {
tmp = tmp - subset(state.rs,well_cells) * cmix_s[oilpos] / d;
}
volRat += tmp / subset(rq_[phase].b,well_cells);
}
// DUMPVAL(cqt_i);
// DUMPVAL(volRat);
// injecting connections total volumerates at std cond
ADB cqt_is = cqt_i/volRat;
@@ -868,6 +896,9 @@ namespace {
cq_s[phase] = cq_ps[phase] + (wops_.w2p * mix_s[phase])*cqt_is;
}
// DUMPVAL(mix_s[2]);
// DUMPVAL(cq_ps[2]);
// Add well contributions to mass balance equations
for (int phase = 0; phase < np; ++phase) {
residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc);
@@ -879,11 +910,8 @@ namespace {
qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
}
residual_.well_flux_eq = qs;
// Updating the well controls is done from here, since we have
// access to the necessary bhp and rate data in this method.
updateWellControls(state.bhp, state.qs, xw);
residual_.well_flux_eq = qs;
}
@@ -956,15 +984,17 @@ namespace {
template<class T>
void FullyImplicitBlackoilSolver<T>::updateWellControls(const ADB& bhp,
const ADB& well_phase_flow_rate,
WellStateFullyImplicitBlackoil& xw) const
void FullyImplicitBlackoilSolver<T>::updateWellControls(ADB& bhp,
ADB& well_phase_flow_rate,
WellStateFullyImplicitBlackoil& xw) const
{
std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
bool bhp_changed = false;
bool rates_changed = false;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
// The current control in the well state overrides
@@ -1002,8 +1032,38 @@ namespace {
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
xw.currentControls()[w] = ctrl_index;
// Also updating well state and primary variables.
// We can only be switching to BHP and SURFACE_RATE
// controls since we do not support RESERVOIR_RATE.
const double target = well_controls_iget_target(wc, ctrl_index);
const double* distr = well_controls_iget_distr(wc, ctrl_index);
switch (well_controls_iget_type(wc, ctrl_index)) {
case BHP:
xw.bhp()[w] = target;
bhp_changed = true;
break;
case SURFACE_RATE:
for (int phase = 0; phase < np; ++phase) {
xw.wellRates()[np*w + phase] = target * distr[phase];
}
rates_changed = true;
break;
default:
OPM_THROW(std::logic_error, "Programmer error: should not have switched "
"to anything but BHP or SURFACE_RATE.");
}
}
}
// Update primary variables, if necessary.
if (bhp_changed) {
ADB::V new_bhp = Eigen::Map<ADB::V>(xw.bhp().data(), nw);
bhp = ADB::function(new_bhp, bhp.derivative());
}
if (rates_changed) {
ADB::V new_qs = Eigen::Map<ADB::V>(xw.wellRates().data(), np*nw);
well_phase_flow_rate = ADB::function(new_qs, well_phase_flow_rate.derivative());
}
}
@@ -1011,7 +1071,8 @@ namespace {
template<class T>
void FullyImplicitBlackoilSolver<T>::addWellControlEq(const SolutionState& state,
const WellStateFullyImplicitBlackoil& xw)
const WellStateFullyImplicitBlackoil& xw,
const V& aliveWells)
{
// Handling BHP and SURFACE_RATE wells.
@@ -1048,6 +1109,17 @@ namespace {
// Choose bhp residual for positive bhp targets.
Selector<double> bhp_selector(bhp_targets);
residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
// For wells that are dead (not flowing), and therefore not communicating
// with the reservoir, we set the equation to be equal to the well's total
// flow. This will be a solution only if the target rate is also zero.
M rate_summer(nw, np*nw);
for (int w = 0; w < nw; ++w) {
for (int phase = 0; phase < np; ++phase) {
rate_summer.insert(w, phase*nw + w) = 1.0;
}
}
Selector<double> alive_selector(aliveWells, Selector<double>::NotEqualZero);
residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs);
// DUMP(residual_.well_eq);
}
@@ -1501,15 +1573,23 @@ namespace {
// convert the pressure offsets to the capillary pressures
std::vector<ADB> pressure = fluid_.capPress(sw, so, sg, cells_);
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
#warning "what's the reference phase??"
// The reference pressure is always the liquid phase (oil) pressure.
if (phaseIdx == BlackoilPhases::Liquid)
continue;
pressure[phaseIdx] = pressure[phaseIdx] - pressure[BlackoilPhases::Liquid];
}
// add the total pressure to the capillary pressures
// Since pcow = po - pw, but pcog = pg - po,
// we have
// pw = po - pcow
// pg = po + pcgo
// This is an unfortunate inconsistency, but a convention we must handle.
for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) {
pressure[phaseIdx] += state.pressure;
if (phaseIdx == BlackoilPhases::Aqua) {
pressure[phaseIdx] = state.pressure - pressure[phaseIdx];
} else {
pressure[phaseIdx] += state.pressure;
}
}
return pressure;

6
opm/autodiff/GeoProps.hpp
View File

@@ -24,8 +24,14 @@
#include <opm/autodiff/GridHelpers.hpp>
//#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/tpfa/TransTpfa.hpp>
#include "disable_warning_pragmas.h"
#include <Eigen/Eigen>
#include "reenable_warning_pragmas.h"
namespace Opm
{

2
opm/autodiff/GridHelpers.cpp
View File

@@ -236,8 +236,6 @@ void extractInternalFaces(const Dune::CpGrid& grid,
typedef Eigen::Array<bool, Eigen::Dynamic, 1> OneColBool;
typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
typedef Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColBool;
ADFaceCellTraits<Dune::CpGrid>::Type nb = faceCells(grid);
// std::cout << "nb = \n" << nb << std::endl;
// Extracts the internal faces of the grid.
// These are stored in internal_faces.
int nf=numFaces(grid);

6
opm/autodiff/GridHelpers.hpp
View File

@@ -25,6 +25,9 @@
#include <boost/range/iterator_range.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/GridHelpers.hpp>
#include "disable_warning_pragmas.h"
#include <Eigen/Eigen>
#include <Eigen/Sparse>
@@ -32,6 +35,9 @@
#include <dune/grid/CpGrid.hpp>
#endif
#include "reenable_warning_pragmas.h"
namespace Opm
{

6
opm/autodiff/ImpesTPFAAD.cpp
View File

@@ -61,7 +61,7 @@ namespace {
std::vector<int> f2hf(2 * numFaces(grid), -1);
Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
face_cells;
face_cells = faceCells(grid);
typedef typename Opm::UgGridHelpers::Cell2FacesTraits<UnstructuredGrid>::Type
Cell2Faces;
@@ -71,7 +71,7 @@ namespace {
cell_faces = c2f[c];
typedef typename Cell2Faces::row_type::iterator Iter;
for (Iter f=cell_faces.begin(), end=cell_faces.end();
f!=end; ++end) {
f!=end; ++f) {
const int p = 0 + (face_cells(*f,0) != c);
f2hf[2*(*f) + p] = f-c2f[0].begin();
}
@@ -89,8 +89,6 @@ namespace {
std::vector<Tri> grav; grav.reserve(2 * ni);
for (HelperOps::IFaces::Index i = 0; i < ni; ++i) {
const int f = ops.internal_faces[ i ];
Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor>
face_cells=faceCells(grid);
const int c1 = face_cells(f,0);
const int c2 = face_cells(f,1);

7
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
View File

@@ -16,7 +16,10 @@
#include <boost/filesystem.hpp>
#ifdef HAVE_DUNE_CORNERPOINT
#include "disable_warning_pragmas.h"
#include <dune/common/version.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include "reenable_warning_pragmas.h"
#endif
namespace Opm
{
@@ -175,7 +178,11 @@ namespace Opm
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step;
#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafGridView(), Dune::VTK::nonconforming);
#else
Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafView(), Dune::VTK::nonconforming);
#endif
writer.addCellData(state.saturation(), "saturation", state.numPhases());
writer.addCellData(state.pressure(), "pressure", 1);