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

Merge pull request #394 from totto82/addNNC

Browse Source 
Add support for NNC in the simulator
This commit is contained in:
Atgeirr Flø Rasmussen 2015-07-10 13:07:01 +02:00
commit 5be3ae9ad6
5 changed files with 98 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

97
opm/autodiff/AutoDiffHelpers.hpp
View File

@ -1,5 +1,6 @@
/*
Copyright 2013 SINTEF ICT, Applied Mathematics.
Copyright 2015 IRIS
This file is part of the Open Porous Media project (OPM).
@ -24,7 +25,8 @@
#include <opm/autodiff/GridHelpers.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <iostream>
#include <vector>
@ -58,41 +60,91 @@ struct HelperOps
/// Extract for each cell the sum of all its adjacent faces' (signed) values.
M fulldiv;
/// Non-neighboring connections
typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
TwoColInt nnc_cells;
/// The NNC transmissibilities
V nnc_trans;
/// Constructs all helper vectors and matrices.
template<class Grid>
HelperOps(const Grid& grid)
HelperOps(const Grid& grid, Opm::EclipseStateConstPtr eclState = EclipseStateConstPtr (nullptr) )
{
using namespace AutoDiffGrid;
const int nc = numCells(grid);
const int nf = numFaces(grid);
// Define some neighbourhood-derived helper arrays.
typedef Eigen::Array<int, Eigen::Dynamic, 2, Eigen::RowMajor> TwoColInt;
TwoColInt nbi;
extractInternalFaces(grid, internal_faces, nbi);
int num_internal=internal_faces.size();
int num_internal = internal_faces.size();
// num_connections may also include non-neighboring connections
int num_connections = num_internal;
int numNNC = 0;
// handle non-neighboring connections
std::shared_ptr<const NNC> nnc = eclState ? eclState->getNNC() : std::make_shared<const NNC> (nullptr, nullptr);
const bool has_nnc = nnc && nnc->hasNNC();
if (has_nnc) {
numNNC = nnc->numNNC();
num_connections += numNNC;
//std::cout << "Added " << numNNC << " NNC" <<std::endl;
nbi.resize(num_internal, 2);
// the nnc's acts on global indicies and must be mapped to cell indicies
size_t cartesianSize = eclState->getEclipseGrid()->getCartesianSize();
std::vector<int> global2localIdx(cartesianSize,0);
for (int i = 0; i< nc; ++i) {
global2localIdx[ globalCell( grid )[i] ] = i;
}
const std::vector<size_t>& NNC1 = nnc->nnc1();
const std::vector<size_t>& NNC2 = nnc->nnc2();
nnc_cells.resize(numNNC,2);
for (int i = 0; i < numNNC; ++i) {
nnc_cells(i,0) = global2localIdx[NNC1[i]];
nnc_cells(i,1) = global2localIdx[NNC2[i]];
}
// store the nnc transmissibilities for later usage.
nnc_trans = Eigen::Map<const V>(nnc->trans().data(), numNNC);
} else {
nnc_trans.resize(0);
nnc_cells.resize(0,0);
}
// std::cout << "nbi = \n" << nbi << std::endl;
// Create matrices.
ngrad.resize(num_internal, nc);
caver.resize(num_internal, nc);
ngrad.resize(num_connections, nc);
caver.resize(num_connections, nc);
typedef Eigen::Triplet<double> Tri;
std::vector<Tri> ngrad_tri;
std::vector<Tri> caver_tri;
ngrad_tri.reserve(2*num_internal);
caver_tri.reserve(2*num_internal);
ngrad_tri.reserve(2*num_connections);
caver_tri.reserve(2*num_connections);
for (int i = 0; i < num_internal; ++i) {
ngrad_tri.emplace_back(i, nbi(i,0), 1.0);
ngrad_tri.emplace_back(i, nbi(i,1), -1.0);
caver_tri.emplace_back(i, nbi(i,0), 0.5);
caver_tri.emplace_back(i, nbi(i,1), 0.5);
}
// add contribution from NNC
if (has_nnc) {
for (int i = 0; i < numNNC; ++i) {
ngrad_tri.emplace_back(i+num_internal, nnc_cells(i,0), 1.0);
ngrad_tri.emplace_back(i+num_internal, nnc_cells(i,1), -1.0);
caver_tri.emplace_back(i+num_internal, nnc_cells(i,0), 0.5);
caver_tri.emplace_back(i+num_internal, nnc_cells(i,1), 0.5);
}
}
ngrad.setFromTriplets(ngrad_tri.begin(), ngrad_tri.end());
caver.setFromTriplets(caver_tri.begin(), caver_tri.end());
grad = -ngrad;
div = ngrad.transpose();
std::vector<Tri> fullngrad_tri;
fullngrad_tri.reserve(2*nf);
typename ADFaceCellTraits<Grid>::Type nb=faceCellsToEigen(grid);
fullngrad_tri.reserve(2*(nf+numNNC));
typename ADFaceCellTraits<Grid>::Type nb = faceCellsToEigen(grid);
for (int i = 0; i < nf; ++i) {
if (nb(i,0) >= 0) {
fullngrad_tri.emplace_back(i, nb(i,0), 1.0);
@ -101,7 +153,14 @@ struct HelperOps
fullngrad_tri.emplace_back(i, nb(i,1), -1.0);
}
}
fullngrad.resize(nf, nc);
// add contribution from NNC
if (nnc->hasNNC()) {
for (int i = 0; i < numNNC; ++i) {
fullngrad_tri.emplace_back(i+nf, nnc_cells(i,0), 1.0);
fullngrad_tri.emplace_back(i+nf, nnc_cells(i,1), -1.0);
}
}
fullngrad.resize(nf+numNNC, nc);
fullngrad.setFromTriplets(fullngrad_tri.begin(), fullngrad_tri.end());
fulldiv = fullngrad.transpose();
}
@ -127,11 +186,15 @@ struct HelperOps
const IFIndex nif = h.internal_faces.size();
typename ADFaceCellTraits<Grid>::Type
face_cells = faceCellsToEigen(g);
assert(nif == ifaceflux.size());
// num connections may possibly include NNCs
int num_nnc = h.nnc_trans.size();
int num_connections = nif + num_nnc;
assert(num_connections == ifaceflux.size());
// Define selector structure.
typedef typename Eigen::Triplet<Scalar> Triplet;
std::vector<Triplet> s; s.reserve(nif);
std::vector<Triplet> s; s.reserve(num_connections);
for (IFIndex iface = 0; iface < nif; ++iface) {
const int f = h.internal_faces[iface];
const int c1 = face_cells(f,0);
@ -144,9 +207,15 @@ struct HelperOps
s.push_back(Triplet(iface, c, Scalar(1)));
}
if (num_nnc > 0) {
for (int i = 0; i < num_nnc; ++i) {
const int c = (ifaceflux[i+nif] >= 0) ? h.nnc_cells(i,0) : h.nnc_cells(i,1);
s.push_back(Triplet(i+nif,c,Scalar(1)));
}
}
// Assemble explicit selector operator.
select_.resize(nif, numCells(g));
select_.resize(num_connections, numCells(g));
select_.setFromTriplets(s.begin(), s.end());
}

4
opm/autodiff/BlackoilModel.hpp
View File

@ -55,6 +55,7 @@ namespace Opm {
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
@ -65,11 +66,12 @@ namespace Opm {
const RockCompressibility* rock_comp_props,
const Wells* wells,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output)
: Base(param, grid, fluid, geo, rock_comp_props, wells, linsolver,
has_disgas, has_vapoil, terminal_output)
eclState, has_disgas, has_vapoil, terminal_output)
{
}
};

3
opm/autodiff/BlackoilModelBase.hpp
View File

@ -31,6 +31,7 @@
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <opm/autodiff/BlackoilModelEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
#include <array>
@ -117,6 +118,7 @@ namespace Opm {
/// \param[in] rock_comp_props if non-null, rock compressibility properties
/// \param[in] wells well structure
/// \param[in] linsolver linear solver
/// \param[in] eclState eclipse state
/// \param[in] has_disgas turn on dissolved gas
/// \param[in] has_vapoil turn on vaporized oil feature
/// \param[in] terminal_output request output to cout/cerr
@ -127,6 +129,7 @@ namespace Opm {
const RockCompressibility* rock_comp_props,
const Wells* wells,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output);

10
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -42,6 +42,7 @@
#include <opm/core/utility/Units.hpp>
#include <opm/core/well_controls.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <cassert>
#include <cmath>
@ -144,6 +145,7 @@ namespace detail {
const RockCompressibility* rock_comp_props,
const Wells* wells,
const NewtonIterationBlackoilInterface& linsolver,
Opm::EclipseStateConstPtr eclState,
const bool has_disgas,
const bool has_vapoil,
const bool terminal_output)
@ -156,7 +158,7 @@ namespace detail {
, active_(detail::activePhases(fluid.phaseUsage()))
, canph_ (detail::active2Canonical(fluid.phaseUsage()))
, cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
, ops_ (grid)
, ops_ (grid, eclState)
, wops_ (wells_)
, has_disgas_(has_disgas)
, has_vapoil_(has_vapoil)
@ -837,9 +839,13 @@ namespace detail {
// Set up the common parts of the mass balance equations
// for each active phase.
const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
const V trans_nnc = ops_.nnc_trans;
V trans_all(transi.size() + trans_nnc.size());
trans_all << transi, trans_nnc;
const std::vector<ADB> kr = computeRelPerm(state);
for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
asImpl().computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
asImpl().computeMassFlux(phaseIdx, trans_all, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state);
residual_.material_balance_eq[ phaseIdx ] =
pvdt_ * (rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0])

1
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -341,6 +341,7 @@ namespace Opm
rock_comp_props_,
wells,
solver_,
eclipse_state_,
has_disgas_,
has_vapoil_,
terminal_output_));