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-07 13:28:08 +02:00
commit b4a49c5844
4 changed files with 69 additions and 48 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
AutoDiffHelpers.hpp
View File

@ -141,17 +141,20 @@ namespace {
// -------------------- upwinding helper class --------------------
/// Upwind selection in absence of counter-current flow (i.e.,
/// without effects of gravity and/or capillary pressure).
template <typename Scalar>
class UpwindSelectorTotalFlux {
class UpwindSelector {
public:
typedef AutoDiff::ForwardBlock<Scalar> ADB;
UpwindSelectorTotalFlux(const UnstructuredGrid& g,
const HelperOps& h,
const typename ADB::V& ifaceflux)
UpwindSelector(const UnstructuredGrid& g,
const HelperOps& h,
const typename ADB::V& ifaceflux)
{
typedef HelperOps::IFaces::Index IFIndex;
const IFIndex nif = h.internal_faces.size();
assert(nif == ifaceflux.size());
// Define selector structure.
typedef typename Eigen::Triplet<Scalar> Triplet;
@ -174,14 +177,10 @@ namespace {
select_.setFromTriplets(s.begin(), s.end());
}
// Upwind selection in absence of counter-current flow (i.e.,
// without effects of gravity and/or capillary pressure).
/// Apply selector to multiple per-cell quantities.
std::vector<ADB>
select(const std::vector<ADB>& xc) const
{
// Apply selector.
//
// Absence of counter-current flow means that the same
// selector applies to all quantities, 'x', defined per
// cell.
@ -195,6 +194,12 @@ namespace {
return xf;
}
/// Apply selector to single per-cell quantity.
ADB select(const ADB& xc) const
{
return select_*xc;
}
private:
typename ADB::M select_;
};

1
SimulatorIncompTwophaseAdfi.cpp
View File

@ -368,6 +368,7 @@ namespace Opm
tsolver_.reset(new Opm::TransportSolverTwophaseAd(grid,
props,
linsolver,
gravity,
param));
} else {
THROW("Unknown transport solver type: " << transport_solver_type_);

56
TransportSolverTwophaseAd.cpp
View File

@ -20,6 +20,7 @@
#include "TransportSolverTwophaseAd.hpp"
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
@ -33,15 +34,18 @@ namespace Opm
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] linsolver Linear solver for Newton-Raphson scheme.
/// \param[in] gravity Gravity vector (null for no gravity).
/// \param[in] param Parameters for the solver.
TransportSolverTwophaseAd::TransportSolverTwophaseAd(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const LinearSolverInterface& linsolver,
const double* gravity,
const parameter::ParameterGroup& param)
: grid_(grid),
props_(props),
linsolver_(linsolver),
ops_(grid),
gravity_(0.0),
tol_(param.getDefault("nl_tolerance", 1e-9)),
maxit_(param.getDefault("nl_maxiter", 30))
{
@ -50,6 +54,23 @@ namespace Opm
for (int i = 0; i < nc; ++i) {
allcells_[i] = i;
}
if (gravity) {
gravity_ = gravity[grid_.dimensions - 1];
for (int dd = 0; dd < grid_.dimensions - 1; ++dd) {
if (gravity[dd] != 0.0) {
THROW("TransportSolverTwophaseAd: can only handle gravity aligned with last dimension");
}
}
V htrans(grid.cell_facepos[grid.number_of_cells]);
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]];
}
}
}
@ -148,22 +169,43 @@ namespace Opm
const TwoCol s0 = Eigen::Map<const TwoCol>(state.saturation().data(), nc, 2);
double res_norm = 1e100;
const V sw0 = s0.leftCols<1>();
// sw1 is the object that will be changed every Newton iteration.
// V sw1 = sw0;
V sw1 = 0.5*V::Ones(nc,1);
const V p1 = Vec(state.pressure().data(), nc, 1);
const V ndp = (ops_.ngrad * p1.matrix()).array();
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]];
}
const UpwindSelectorTotalFlux<double> upwind(grid_, ops_, dflux);
// 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.
// 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);
const V ndp = (ops_.ngrad * p1.matrix()).array();
typedef Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> DynArr;
const V z = Eigen::Map<DynArr>(grid_.cell_centroids, nc, grid_.dimensions).rightCols<1>();
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);
// Compute more explicit and constant terms used in the equations.
const V pv = Vec(porevolume, nc, 1);
const V dtpv = dt/pv;
const V q = Vec(source, nc, 1);
const V qneg = q.min(V::Zero(nc,1));
const V qpos = q.max(V::Zero(nc,1));
const double gfactor = gravity_*(density[0] - density[1]);
const V gravflux = ndz*transi_*gfactor;
// Block pattern for variables.
// Primary variables:
@ -176,12 +218,14 @@ namespace Opm
// Assemble linear system.
const ADB sw = ADB::variable(0, sw1, bpat);
const std::vector<ADB> pmobc = phaseMobility<ADB>(props_, allcells_, sw.value());
const std::vector<ADB> pmobf = upwind.select(pmobc);
const ADB fw_cell = fluxFunc(pmobc);
// const ADB fw_face = upwind.select(fw_cell);
const std::vector<ADB> pmobf = { upwind_w.select(pmobc[0]),
upwind_o.select(pmobc[1]) };
const ADB fw_face = fluxFunc(pmobf);
const ADB flux1 = fw_face * dflux;
const ADB flux = fw_face * (dflux - pmobf[1]*gravflux);
const ADB qtr_ad = qpos + fw_cell*qneg;
const ADB transport_residual = sw - sw0 + dtpv*(ops_.div*flux1 - qtr_ad);
const ADB transport_residual = sw - sw0 + dtpv*(ops_.div*flux - qtr_ad);
res_norm = transport_residual.value().matrix().norm();
std::cout << "Residual l2-norm = " << res_norm << std::endl;

37
TransportSolverTwophaseAd.hpp
View File

@ -43,10 +43,12 @@ namespace Opm
/// \param[in] grid A 2d or 3d grid.
/// \param[in] props Rock and fluid properties.
/// \param[in] linsolver Linear solver for Newton-Raphson scheme.
/// \param[in] gravity Gravity vector (null for no gravity).
/// \param[in] param Parameters for the solver.
TransportSolverTwophaseAd(const UnstructuredGrid& grid,
const IncompPropertiesInterface& props,
const LinearSolverInterface& linsolver,
const double* gravity,
const parameter::ParameterGroup& param);
// Virtual destructor.
@ -75,42 +77,11 @@ namespace Opm
const IncompPropertiesInterface& props_;
const LinearSolverInterface& linsolver_;
const HelperOps ops_;
double gravity_;
double tol_;
int maxit_;
std::vector<int> allcells_;
#if 0
const double* visc_;
std::vector<double> smin_;
std::vector<double> smax_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one source per cell
double dt_;
std::vector<double> saturation_; // one per cell, only water saturation!
std::vector<double> fractionalflow_; // = m[0]/(m[0] + m[1]) per cell
std::vector<int> reorder_iterations_;
//std::vector<double> reorder_fval_;
// For gravity segregation.
std::vector<double> gravflux_;
std::vector<double> mob_;
std::vector<double> s0_;
std::vector<std::vector<int> > columns_;
// Storing the upwind and downwind graphs for experiments.
std::vector<int> ia_upw_;
std::vector<int> ja_upw_;
std::vector<int> ia_downw_;
std::vector<int> ja_downw_;
struct Residual;
double fracFlow(double s, int cell) const;
struct GravityResidual;
void mobility(double s, int cell, double* mob) const;
#endif
V transi_;
};
} // namespace Opm