Loading...
Searching...
No Matches
UpscalerBase_impl.hpp
//===========================================================================
//
// File: UpscalerBase_impl.hpp
//
// Created: Thu Apr 29 10:22:06 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
Copyright 2010 Statoil ASA.
This file is part of The Open Reservoir Simulator Project (OpenRS).
OpenRS 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.
OpenRS 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 OpenRS. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_UPSCALERBASE_IMPL_HEADER
#define OPM_UPSCALERBASE_IMPL_HEADER
#include <opm/porsol/common/setupGridAndProps.hpp>
#include <opm/porsol/common/setupBoundaryConditions.hpp>
#include <opm/porsol/common/ReservoirPropertyTracerFluid.hpp>
#include <iostream>
namespace Opm
{
template <class Traits>
inline UpscalerBase<Traits>::UpscalerBase()
: bctype_(Fixed),
twodim_hack_(false),
residual_tolerance_(1e-8),
linsolver_maxit_(0),
linsolver_prolongate_factor_(1.0),
linsolver_verbosity_(0),
linsolver_type_(3),
linsolver_smooth_steps_(1),
gravity_(0.0)
#if HAVE_MPI
, grid_(MPI_COMM_SELF)
#endif
{
}
template <class Traits>
inline void UpscalerBase<Traits>::init(const Opm::ParameterGroup& param)
{
initImpl(param);
initFinal(param);
}
template <class Traits>
inline void UpscalerBase<Traits>::initImpl(const Opm::ParameterGroup& param)
{
// Request the boundary condition type parameter early since,
// depending on the actual type, we may have to manufacture
// and insert other parameters into the ParameterGroup prior
// to constructing the grid and associated properties.
//
int bct = param.get<int>("boundary_condition_type");
bctype_ = static_cast<BoundaryConditionType>(bct);
// Import control parameters pertaining to reduced physical
// dimensionality ("2d_hack = true" precludes periodic
// boundary conditions in the Z direction), and linear solves.
//
twodim_hack_ = param.getDefault("2d_hack", twodim_hack_);
residual_tolerance_ = param.getDefault("residual_tolerance", residual_tolerance_);
linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_);
linsolver_type_ = param.getDefault("linsolver_type", linsolver_type_);
linsolver_maxit_ = param.getDefault("linsolver_max_iterations", linsolver_maxit_);
linsolver_prolongate_factor_ = param.getDefault("linsolver_prolongate_factor", linsolver_prolongate_factor_);
linsolver_smooth_steps_ = param.getDefault("linsolver_smooth_steps", linsolver_smooth_steps_);
gravity_ = param.getDefault("gravity", gravity_);
// Ensure sufficient grid support for requested boundary
// condition type.
//
Opm::ParameterGroup temp_param = param;
if (bctype_ == Linear || bctype_ == Periodic) {
if (!temp_param.has("use_unique_boundary_ids")) {
temp_param.insertParameter("use_unique_boundary_ids", "true");
}
if (!temp_param.has("clip_z")) {
temp_param.insertParameter("clip_z", "true");
}
}
if (bctype_ == Periodic) {
if (!temp_param.has("periodic_extension")) {
temp_param.insertParameter("periodic_extension", "true");
}
}
setupGridAndProps(temp_param, grid_, res_prop_);
ginterf_.init(grid_);
}
template <class Traits>
inline void UpscalerBase<Traits>::initFinal(const Opm::ParameterGroup& param)
{
// Report any unused parameters.
std::cout << "==================== Unused parameters: ====================\n";
param.displayUsage();
std::cout << "================================================================\n";
}
template <class Traits>
inline void UpscalerBase<Traits>::init(const Opm::Deck& deck,
BoundaryConditionType bctype,
double perm_threshold,
double residual_tolerance,
int linsolver_verbosity,
int linsolver_type,
bool twodim_hack,
int linsolver_maxit,
double linsolver_prolongate_factor,
int linsolver_smooth_steps,
const double gravity)
{
bctype_ = bctype;
residual_tolerance_ = residual_tolerance;
linsolver_verbosity_ = linsolver_verbosity;
linsolver_type_ = linsolver_type;
linsolver_maxit_ = linsolver_maxit;
linsolver_prolongate_factor_ = linsolver_prolongate_factor;
linsolver_smooth_steps_ = linsolver_smooth_steps;
twodim_hack_ = twodim_hack;
gravity_ = gravity;
// Faking some parameters depending on bc type.
bool periodic_ext = (bctype_ == Periodic);
bool turn_normals = false;
bool clip_z = (bctype_ == Linear || bctype_ == Periodic);
bool unique_bids = (bctype_ == Linear || bctype_ == Periodic);
std::string rock_list("no_list");
setupGridAndPropsEclipse(deck,
periodic_ext, turn_normals, clip_z, unique_bids,
perm_threshold, rock_list,
useJ<ResProp>(), 1.0, 0.0,
grid_, res_prop_);
ginterf_.init(grid_);
}
template <class Traits>
inline const typename UpscalerBase<Traits>::GridType&
UpscalerBase<Traits>::grid() const
{
return grid_;
}
template <class Traits>
inline void
UpscalerBase<Traits>::setBoundaryConditionType(BoundaryConditionType type)
{
if ((type == Periodic && bctype_ != Periodic)
|| (type != Periodic && bctype_ == Periodic)) {
OPM_THROW(std::runtime_error, "Cannot switch to or from Periodic boundary condition, "
"periodic must be set in init() params.");
} else {
bctype_ = type;
if (type == Periodic || type == Linear) {
grid_.setUniqueBoundaryIds(true);
} else {
grid_.setUniqueBoundaryIds(false);
}
}
}
template <class Traits>
inline void
UpscalerBase<Traits>::setPermeability(const int cell_index, const permtensor_t& k)
{
res_prop_.permeabilityModifiable(cell_index) = k;
}
template <class Traits>
inline typename UpscalerBase<Traits>::permtensor_t
UpscalerBase<Traits>::upscaleSinglePhase()
{
ReservoirPropertyTracerFluid fluid;
return upscaleEffectivePerm(fluid);
}
template <class Traits>
template <class FluidInterface>
inline typename UpscalerBase<Traits>::permtensor_t
UpscalerBase<Traits>::upscaleEffectivePerm(const FluidInterface& fluid)
{
int num_cells = ginterf_.numberOfCells();
// No source or sink.
std::vector<double> src(num_cells, 0.0);
// Just water.
std::vector<double> sat(num_cells, 1.0);
// Gravity.
Dune::FieldVector<double, 3> gravity(0.0);
gravity[2] = gravity_;
permtensor_t upscaled_K(3, 3, (double*)0);
for (int pdd = 0; pdd < Dimension; ++pdd) {
setupUpscalingConditions(ginterf_, bctype_, pdd, 1.0, 1.0, twodim_hack_, bcond_);
if (pdd == 0) {
// Only on first iteration, since we do not change the
// structure of the system, the way the flow solver is
// implemented.
flow_solver_.init(ginterf_, res_prop_, gravity, bcond_);
}
// Run pressure solver.
bool same_matrix = (bctype_ != Fixed) && (pdd != 0);
flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_,
linsolver_verbosity_,
linsolver_type_, same_matrix,
linsolver_maxit_, linsolver_prolongate_factor_,
linsolver_smooth_steps_);
double max_mod = flow_solver_.postProcessFluxes();
std::cout << "Max mod = " << max_mod << std::endl;
// Compute upscaled K.
double Q[Dimension] = { 0 };
switch (bctype_) {
case Fixed:
Q[pdd] = computeAverageVelocity(flow_solver_.getSolution(), pdd, pdd);
break;
case Linear:
case Periodic:
for (int i = 0; i < Dimension; ++i) {
Q[i] = computeAverageVelocity(flow_solver_.getSolution(), i, pdd);
}
break;
default:
OPM_THROW(std::runtime_error, "Unknown boundary type: " + std::to_string(bctype_));
}
double delta = computeDelta(pdd);
for (int i = 0; i < Dimension; ++i) {
upscaled_K(i, pdd) = Q[i] * delta;
}
}
return upscaled_K;
}
template <class Traits>
template <class FlowSol>
double UpscalerBase<Traits>::computeAverageVelocity(const FlowSol& flow_solution,
const int flow_dir,
const int pdrop_dir) const
{
double side1_flux = 0.0;
double side2_flux = 0.0;
double side1_area = 0.0;
double side2_area = 0.0;
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
if (f->boundary()) {
int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
if ((canon_bid - 1)/2 == flow_dir) {
double flux = flow_solution.outflux(f);
double area = f->area();
double norm_comp = f->normal()[flow_dir];
// std::cout << "bid " << f->boundaryId() << " area " << area << " n " << norm_comp << std::endl;
if (canon_bid - 1 == 2*flow_dir) {
if (flow_dir == pdrop_dir && flux > 0.0) {
#ifdef VERBOSE
std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()<<" (canonical: "<<canon_bid
<< ") Magnitude: " << std::fabs(flux) << std::endl;
#endif
// OPM_THROW(std::runtime_error, "Detected outflow at entry face: " << face);
}
side1_flux += flux*norm_comp;
side1_area += area;
} else {
assert(canon_bid - 1 == 2*flow_dir + 1);
if (flow_dir == pdrop_dir && flux < 0.0) {
#ifdef VERBOSE
std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()
<< " Magnitude: " << std::fabs(flux) << std::endl;
#endif
// OPM_THROW(std::runtime_error, "Detected inflow at exit face: " << face);
}
side2_flux += flux*norm_comp;
side2_area += area;
}
}
}
}
}
// q is the average velocity.
return 0.5*(side1_flux/side1_area + side2_flux/side2_area);
}
template <class Traits>
inline double UpscalerBase<Traits>::computeDelta(const int flow_dir) const
{
double side1_pos = 0.0;
double side2_pos = 0.0;
double side1_area = 0.0;
double side2_area = 0.0;
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
if (f->boundary()) {
int canon_bid = bcond_.getCanonicalBoundaryId(f->boundaryId());
if ((canon_bid - 1)/2 == flow_dir) {
double area = f->area();
double pos_comp = f->centroid()[flow_dir];
if (canon_bid - 1 == 2*flow_dir) {
side1_pos += area*pos_comp;
side1_area += area;
} else {
side2_pos += area*pos_comp;
side2_area += area;
}
}
}
}
}
// delta is the average length.
return side2_pos/side2_area - side1_pos/side1_area;
}
template <class Traits>
double UpscalerBase<Traits>::upscalePorosity() const
{
double total_vol = 0.0;
double total_pore_vol = 0.0;
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
total_vol += c->volume();
total_pore_vol += c->volume()*res_prop_.porosity(c->index());
}
return total_pore_vol/total_vol;
}
template <class Traits>
double UpscalerBase<Traits>::upscaleNetPorosity() const
{
double total_net_vol = 0.0;
double total_pore_vol = 0.0;
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
total_net_vol += c->volume()*res_prop_.ntg(c->index());
total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
}
if (total_net_vol>0.0) return total_pore_vol/total_net_vol;
else return 0.0;
}
template <class Traits>
double UpscalerBase<Traits>::upscaleNTG() const
{
double total_vol = 0.0;
double total_net_vol = 0.0;
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
total_vol += c->volume();
total_net_vol += c->volume()*res_prop_.ntg(c->index());
}
return total_net_vol/total_vol;
}
template <class Traits>
double UpscalerBase<Traits>::upscaleSWCR(const bool NTG) const
{
double total_swcr = 0.0;
double total_pore_vol = 0.0;
if (NTG) {
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.swcr(c->index());
total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
}
}
else {
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
total_swcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.swcr(c->index());
total_pore_vol += c->volume()*res_prop_.porosity(c->index());
}
}
return total_swcr/total_pore_vol;
}
template <class Traits>
double UpscalerBase<Traits>::upscaleSOWCR(const bool NTG) const
{
double total_sowcr = 0.0;
double total_pore_vol = 0.0;
if (NTG) {
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index())*res_prop_.sowcr(c->index());
total_pore_vol += c->volume()*res_prop_.porosity(c->index())*res_prop_.ntg(c->index());
}
}
else {
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
total_sowcr += c->volume()*res_prop_.porosity(c->index())*res_prop_.sowcr(c->index());
total_pore_vol += c->volume()*res_prop_.porosity(c->index());
}
}
return total_sowcr/total_pore_vol;
}
} // namespace Opm
#endif // OPM_UPSCALERBASE_IMPL_HEADER
Opm::ReservoirPropertyTracerFluid
Definition ReservoirPropertyTracerFluid.hpp:40
Opm::UpscalerBase::upscaleSinglePhase
permtensor_t upscaleSinglePhase()
Does a single-phase upscaling.
Definition UpscalerBase_impl.hpp:222
Opm::UpscalerBase::setPermeability
void setPermeability(const int cell_index, const permtensor_t &k)
Set the permeability of a cell directly.
Definition UpscalerBase_impl.hpp:212
Opm::UpscalerBase::upscalePorosity
double upscalePorosity() const
Compute upscaled porosity.
Definition UpscalerBase_impl.hpp:377
Opm::UpscalerBase::permtensor_t
ResProp::MutablePermTensor permtensor_t
A type for the upscaled permeability.
Definition UpscalerBase.hpp:66
Opm::UpscalerBase::grid
const GridType & grid() const
Access the grid.
Definition UpscalerBase_impl.hpp:181
Opm::UpscalerBase::upscaleNTG
double upscaleNTG() const
Compute upscaled NTG.
Definition UpscalerBase_impl.hpp:403
Opm::UpscalerBase::upscaleSOWCR
double upscaleSOWCR(const bool NTG) const
Compute upscaled SOWCR.
Definition UpscalerBase_impl.hpp:435
Opm::UpscalerBase::init
void init(const Opm::ParameterGroup &param)
Initializes the upscaler from parameters.
Definition UpscalerBase_impl.hpp:68
Opm::UpscalerBase::setBoundaryConditionType
void setBoundaryConditionType(BoundaryConditionType type)
Set boundary condition type.
Definition UpscalerBase_impl.hpp:191
Opm::UpscalerBase::upscaleNetPorosity
double upscaleNetPorosity() const
Compute upscaled net porosity.
Definition UpscalerBase_impl.hpp:390
Opm::UpscalerBase::upscaleSWCR
double upscaleSWCR(const bool NTG) const
Compute upscaled SWCR.
Definition UpscalerBase_impl.hpp:415
Opm::UpscalerBase::UpscalerBase
UpscalerBase()
Default constructor.
Definition UpscalerBase_impl.hpp:48
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::useJ
bool useJ()
Helper for determining whether we should.
Definition setupGridAndProps.hpp:59
Opm::setupGridAndPropsEclipse
void setupGridAndPropsEclipse(const Opm::Deck &deck, bool periodic_extension, bool turn_normals, bool clip_z, bool unique_bids, double perm_threshold, const std::string &rock_list, bool use_jfunction_scaling, double sigma, double theta, Dune::CpGrid &grid, ResProp< 3 > &res_prop)
Definition setupGridAndProps.hpp:152
Opm::setupUpscalingConditions
void setupUpscalingConditions(const GridInterface &g, int bct, int pddir, double pdrop, double bdy_sat, bool twodim_hack, BCs &bcs)
Definition setupBoundaryConditions.hpp:99
Opm::setupGridAndProps
void setupGridAndProps(const Opm::ParameterGroup &param, Dune::CpGrid &grid, ResProp< 3 > &res_prop)
Definition setupGridAndProps.hpp:71