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/master'

Browse Source
This commit is contained in:
Bård Skaflestad
2013-05-08 13:02:59 +02:00
parent ea1cdcd478 a797ce5f32
commit dfd33d8973
2 changed files with 49 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
AutoDiffHelpers.hpp
View File

@@ -194,16 +194,57 @@ namespace {
return xf;
}
/// Apply selector to single per-cell quantity.
/// Apply selector to single per-cell ADB quantity.
ADB select(const ADB& xc) const
{
return select_*xc;
}
/// Apply selector to single per-cell constant quantity.
typename ADB::V select(const typename ADB::V& xc) const
{
return (select_*xc.matrix()).array();
}
private:
typename ADB::M select_;
};
}
namespace {
template <typename Scalar, class IntVec>
Eigen::SparseMatrix<Scalar>
constructSubsetSparseMatrix(const int full_size, const IntVec& indices)
{
typedef Eigen::Triplet<Scalar> Tri;
const int subset_size = indices.size();
std::vector<Tri> triplets(subset_size);
for (int i = 0; i < subset_size; ++i) {
triplets[i] = Tri(i, indices[i], 1);
}
Eigen::SparseMatrix<Scalar> sub(subset_size, full_size);
sub.setFromTriplets(triplets.begin(), triplets.end());
return sub;
}
} // anon namespace
template <typename Scalar, class IntVec>
AutoDiff::ForwardBlock<Scalar>
subset(const AutoDiff::ForwardBlock<Scalar>& x,
const IntVec& indices)
{
return ::constructSubsetSparseMatrix<Scalar>(x.value().size(), indices) * x;
}
template <typename Scalar, class IntVec>
Eigen::Array<Scalar, Eigen::Dynamic, 1>
subset(const Eigen::Array<Scalar, Eigen::Dynamic, 1>& x,
const IntVec& indices)
{
return (::constructSubsetSparseMatrix<Scalar>(x.size(), indices) * x.matrix()).array();
}
#endif // OPM_AUTODIFFHELPERS_HEADER_INCLUDED

21
TransportSolverTwophaseAd.cpp
View File

@@ -65,11 +65,7 @@ namespace Opm
tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid), props.permeability(), htrans.data());
V trans(grid_.number_of_faces);
tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid), htrans.data(), trans.data());
const int num_internal = ops_.internal_faces.size();
transi_.resize(num_internal);
for (int fi = 0; fi < num_internal; ++fi) {
transi_[fi] = trans[ops_.internal_faces[fi]];
}
transi_ = subset(trans, ops_.internal_faces);
}
}
@@ -174,16 +170,13 @@ namespace Opm
V sw1 = 0.5*V::Ones(nc,1);
const V dflux_all = Vec(state.faceflux().data(), grid_.number_of_faces, 1);
const int num_internal = ops_.internal_faces.size();
V dflux(num_internal);
for (int fi = 0; fi < num_internal; ++fi) {
dflux[fi] = dflux_all[ops_.internal_faces[fi]];
}
V dflux = subset(dflux_all, ops_.internal_faces);
// Upwind selection of mobilities by phase.
// We have that for a phase P
// v_P = lambda_P K (-grad p + rho_P g grad z)
// and we assume that this has the same direction as
// v_P' = -grad p + rho_P g grad z.
// dh_P = -grad p + rho_P g grad z.
// This may not be true for arbitrary anisotropic situations,
// but for scalar lambda and using TPFA it holds.
const V p1 = Vec(state.pressure().data(), nc, 1);
@@ -193,10 +186,10 @@ namespace Opm
const V ndz = (ops_.ngrad * z.matrix()).array();
ASSERT(num_internal == ndp.size());
const double* density = props_.density();
const V vw = ndp - ndz*(gravity_*density[0]);
const V vo = ndp - ndz*(gravity_*density[1]);
const UpwindSelector<double> upwind_w(grid_, ops_, vw);
const UpwindSelector<double> upwind_o(grid_, ops_, vo);
const V dhw = ndp - ndz*(gravity_*density[0]);
const V dho = ndp - ndz*(gravity_*density[1]);
const UpwindSelector<double> upwind_w(grid_, ops_, dhw);
const UpwindSelector<double> upwind_o(grid_, ops_, dho);
// Compute more explicit and constant terms used in the equations.
const V pv = Vec(porevolume, nc, 1);