Loading...
Searching...
No Matches
SteadyStateUpscaler_impl.hpp
//===========================================================================
//
// File: SteadyStateUpscaler_impl.hpp
//
// Created: Fri Aug 28 14:07:51 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Brd Skaflestad <bard.skaflestad@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
Copyright 2009, 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_STEADYSTATEUPSCALER_IMPL_HEADER
#define OPM_STEADYSTATEUPSCALER_IMPL_HEADER
#include <opm/porsol/common/MatrixInverse.hpp>
#include <opm/porsol/common/SimulatorUtilities.hpp>
#include <opm/porsol/common/ReservoirPropertyFixedMobility.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <algorithm>
#include <iostream>
namespace Opm
{
template <class Traits>
inline SteadyStateUpscaler<Traits>::SteadyStateUpscaler()
: Super(),
use_gravity_(false),
output_vtk_(false),
print_inoutflows_(false),
simulation_steps_(10),
stepsize_(0.1),
relperm_threshold_(1.0e-8),
maximum_mobility_contrast_(1.0e9),
sat_change_threshold_(0.0)
{
}
template <class Traits>
inline void SteadyStateUpscaler<Traits>::initImpl(const Opm::parameter::ParameterGroup& param)
{
Super::initImpl(param);
use_gravity_ = param.getDefault("use_gravity", use_gravity_);
output_vtk_ = param.getDefault("output_vtk", output_vtk_);
print_inoutflows_ = param.getDefault("print_inoutflows", print_inoutflows_);
simulation_steps_ = param.getDefault("simulation_steps", simulation_steps_);
stepsize_ = Opm::unit::convert::from(param.getDefault("stepsize", stepsize_),
Opm::unit::day);
relperm_threshold_ = param.getDefault("relperm_threshold", relperm_threshold_);
maximum_mobility_contrast_ = param.getDefault("maximum_mobility_contrast", maximum_mobility_contrast_);
sat_change_threshold_ = param.getDefault("sat_change_threshold", sat_change_threshold_);
transport_solver_.init(param);
// Set viscosities and densities if given.
double v1_default = this->res_prop_.viscosityFirstPhase();
double v2_default = this->res_prop_.viscositySecondPhase();
this->res_prop_.setViscosities(param.getDefault("viscosity1", v1_default), param.getDefault("viscosity2", v2_default));
double d1_default = this->res_prop_.densityFirstPhase();
double d2_default = this->res_prop_.densitySecondPhase();
this->res_prop_.setDensities(param.getDefault("density1", d1_default), param.getDefault("density2", d2_default));
}
namespace {
double maxMobility(double m1, double m2)
{
return std::max(m1, m2);
}
// The matrix variant expects diagonal mobilities.
template <class SomeMatrixType>
double maxMobility(double m1, SomeMatrixType& m2)
{
double m = m1;
for (int i = 0; i < std::min(m2.numRows(), m2.numCols()); ++i) {
m = std::max(m, m2(i,i));
}
return m;
}
void thresholdMobility(double& m, double threshold)
{
m = std::max(m, threshold);
}
// The matrix variant expects diagonal mobilities.
template <class SomeMatrixType>
void thresholdMobility(SomeMatrixType& m, double threshold)
{
for (int i = 0; i < std::min(m.numRows(), m.numCols()); ++i) {
m(i, i) = std::max(m(i, i), threshold);
}
}
} // anon namespace
template <class Traits>
inline std::pair<typename SteadyStateUpscaler<Traits>::permtensor_t,
typename SteadyStateUpscaler<Traits>::permtensor_t>
SteadyStateUpscaler<Traits>::
upscaleSteadyState(const int flow_direction,
const std::vector<double>& initial_saturation,
const double boundary_saturation,
const double pressure_drop,
const permtensor_t& upscaled_perm)
{
static int count = 0;
++count;
int num_cells = this->ginterf_.numberOfCells();
// No source or sink.
std::vector<double> src(num_cells, 0.0);
Opm::SparseVector<double> injection(num_cells);
// Gravity.
Dune::FieldVector<double, 3> gravity(0.0);
if (use_gravity_) {
gravity[2] = Opm::unit::gravity;
}
if (gravity.two_norm() > 0.0) {
OPM_MESSAGE("Warning: Gravity is experimental for flow solver.");
}
// Set up initial saturation profile.
std::vector<double> saturation = initial_saturation;
// Set up boundary conditions.
setupUpscalingConditions(this->ginterf_, this->bctype_, flow_direction,
pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
// Set up solvers.
if (flow_direction == 0) {
this->flow_solver_.init(this->ginterf_, this->res_prop_, gravity, this->bcond_);
}
transport_solver_.initObj(this->ginterf_, this->res_prop_, this->bcond_);
// Run pressure solver.
this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
this->residual_tolerance_, this->linsolver_verbosity_,
this->linsolver_type_, false,
this->linsolver_maxit_, this->linsolver_prolongate_factor_,
this->linsolver_smooth_steps_);
double max_mod = this->flow_solver_.postProcessFluxes();
std::cout << "Max mod = " << max_mod << std::endl;
// Do a run till steady state. For now, we just do some pressure and transport steps...
std::vector<double> saturation_old = saturation;
for (int iter = 0; iter < simulation_steps_; ++iter) {
// Run transport solver.
transport_solver_.transportSolve(saturation, stepsize_, gravity, this->flow_solver_.getSolution(), injection);
// Run pressure solver.
this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
this->residual_tolerance_, this->linsolver_verbosity_,
this->linsolver_type_, false,
this->linsolver_maxit_, this->linsolver_prolongate_factor_,
this->linsolver_smooth_steps_);
max_mod = this->flow_solver_.postProcessFluxes();
std::cout << "Max mod = " << max_mod << std::endl;
// Print in-out flows if requested.
if (print_inoutflows_) {
std::pair<double, double> w_io, o_io;
computeInOutFlows(w_io, o_io, this->flow_solver_.getSolution(), saturation);
std::cout << "Pressure step " << iter
<< "\nWater flow [in] " << w_io.first
<< " [out] " << w_io.second
<< "\nOil flow [in] " << o_io.first
<< " [out] " << o_io.second
<< std::endl;
}
// Output.
if (output_vtk_) {
writeVtkOutput(this->ginterf_,
this->res_prop_,
this->flow_solver_.getSolution(),
saturation,
std::string("output-steadystate")
+ '-' + std::to_string(count)
+ '-' + std::to_string(flow_direction)
+ '-' + std::to_string(iter));
}
// Comparing old to new.
int num_cells = saturation.size();
double maxdiff = 0.0;
for (int i = 0; i < num_cells; ++i) {
maxdiff = std::max(maxdiff, std::fabs(saturation[i] - saturation_old[i]));
}
#ifdef VERBOSE
std::cout << "Maximum saturation change: " << maxdiff << std::endl;
#endif
if (maxdiff < sat_change_threshold_) {
#ifdef VERBOSE
std::cout << "Maximum saturation change is under steady state threshold." << std::endl;
#endif
break;
}
// Copy to old.
saturation_old = saturation;
}
// Compute phase mobilities.
// First: compute maximal mobilities.
typedef typename Super::ResProp::Mobility Mob;
Mob m;
double m1max = 0;
double m2max = 0;
for (int c = 0; c < num_cells; ++c) {
this->res_prop_.phaseMobility(0, c, saturation[c], m.mob);
m1max = maxMobility(m1max, m.mob);
this->res_prop_.phaseMobility(1, c, saturation[c], m.mob);
m2max = maxMobility(m2max, m.mob);
}
// Second: set thresholds.
const double mob1_abs_thres = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
const double mob1_rel_thres = m1max / maximum_mobility_contrast_;
const double mob1_threshold = std::max(mob1_abs_thres, mob1_rel_thres);
const double mob2_abs_thres = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
const double mob2_rel_thres = m2max / maximum_mobility_contrast_;
const double mob2_threshold = std::max(mob2_abs_thres, mob2_rel_thres);
// Third: extract and threshold.
std::vector<Mob> mob1(num_cells);
std::vector<Mob> mob2(num_cells);
for (int c = 0; c < num_cells; ++c) {
this->res_prop_.phaseMobility(0, c, saturation[c], mob1[c].mob);
thresholdMobility(mob1[c].mob, mob1_threshold);
this->res_prop_.phaseMobility(1, c, saturation[c], mob2[c].mob);
thresholdMobility(mob2[c].mob, mob2_threshold);
}
// Compute upscaled relperm for each phase.
ReservoirPropertyFixedMobility<Mob> fluid_first(mob1);
permtensor_t eff_Kw = Super::upscaleEffectivePerm(fluid_first);
ReservoirPropertyFixedMobility<Mob> fluid_second(mob2);
permtensor_t eff_Ko = Super::upscaleEffectivePerm(fluid_second);
// Set the steady state saturation fields for eventual outside access.
last_saturation_state_.swap(saturation);
// Compute the (anisotropic) upscaled mobilities.
// eff_Kw := lambda_w*K
// => lambda_w = eff_Kw*inv(K);
permtensor_t lambda_w(matprod(eff_Kw, inverse3x3(upscaled_perm)));
permtensor_t lambda_o(matprod(eff_Ko, inverse3x3(upscaled_perm)));
// Compute (anisotropic) upscaled relative permeabilities.
// lambda = k_r/mu
permtensor_t k_rw(lambda_w);
k_rw *= this->res_prop_.viscosityFirstPhase();
permtensor_t k_ro(lambda_o);
k_ro *= this->res_prop_.viscositySecondPhase();
return std::make_pair(k_rw, k_ro);
}
template <class Traits>
inline const std::vector<double>&
SteadyStateUpscaler<Traits>::lastSaturationState() const
{
return last_saturation_state_;
}
template <class Traits>
double SteadyStateUpscaler<Traits>::lastSaturationUpscaled() const
{
typedef typename GridInterface::CellIterator CellIter;
double pore_vol = 0.0;
double sat_vol = 0.0;
for (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
pore_vol += cell_pore_vol;
sat_vol += cell_pore_vol*last_saturation_state_[c->index()];
}
// Dividing by pore volume gives average saturations.
return sat_vol/pore_vol;
}
static inline double calc_frac_flow (
std::map< int, double >& frac_flow_by_bid,
const BCs& bcond,
const ResProp& res_prop,
GridInterface::CellIterator c,
CellIterator::FaceIterator f ) {
const SatBC& sc = bcond.satCond( f );
if( sc.isPeriodic() ) {
assert( sc.isDirichlet() );
return res_prop.fractionalFlow( c->index(), sc.saturation() );
}
assert( sc.saturationDifference() == 0.0 );
const int partner_bid = bcond.getPeriodicPartner( f->boundaryId() );
const auto it = frac_flow_by_bid.find( partner_bid );
// This is for the periodic case, so that we are sure all fractional flows have
// been set in frac_flow_by_bid.
if( it == frac_flow_by_bid.end() ) {
OPM_THROW(std::runtime_error, "Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
<< " and partner bid = " << partner_bid);
}
return it->second;
}
template <class Traits>
template <class FlowSol>
void SteadyStateUpscaler<Traits>::computeInOutFlows(std::pair<double, double>& water_inout,
std::pair<double, double>& oil_inout,
const FlowSol& flow_solution,
const std::vector<double>& saturations) const
{
struct sideflux {
double water;
double oil;
sideflux& operator+=( const sideflux& other ) {
this->water += other.water;
this->oil += other.oil;
return *this;
}
};
std::map< int, double > frac_flow_by_bid;
// Two passes: First pass, deal with outflow
sideflux outflow;
for( auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
const double frac_flow = this->res_prop_.fractional_flow( c->index(), saturations[ c->index() ] );
for( auto f = c->facebegin(); f != c->faceend(); ++f ) {
if( !f->boundary() ) continue;
const double flux = flow_solution.outflux( f );
if( flux < 0.0 ) continue;
const SatBC& sc = this->bcond_.satCond( f );
if( sc.isPeriodic() ) frac_flow_by_bid[ f->boundaryId() ] = frac_flow;
outflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
}
}
// second pass, deal with inflow.
sideflux inflow;
for( auto c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c ) {
for( auto f = c->facebegin(); f != c->faceend(); ++f ) {
if( !f->boundary() ) continue;
const double flux = flow_solution.outflux( f );
if( flux >= 0.0 ) continue;
const double frac_flow = calc_frac_flow( frac_flow_by_bid, this->bcond_, this->res_prop_, c, f );
inflow += { flux * frac_flow, flux * ( 1.0 - frac_flow ) };
}
}
water_inout = std::make_pair( in.water, out.water );
oil_inout = std::make_pair( in.oil, out.oil );
}
} // namespace Opm
#endif // OPM_STEADYSTATEUPSCALER_IMPL_HEADER
Opm::BCBase::isPeriodic
bool isPeriodic() const
Type query.
Definition BoundaryConditions.hpp:97
Opm::BCBase::isDirichlet
bool isDirichlet() const
Type query.
Definition BoundaryConditions.hpp:85
Opm::ReservoirPropertyFixedMobility
Definition ReservoirPropertyFixedMobility.hpp:46
Opm::SatBC
A class for representing a saturation boundary condition.
Definition BoundaryConditions.hpp:176
Opm::SatBC::saturation
double saturation() const
Query a Dirichlet condition.
Definition BoundaryConditions.hpp:197
Opm::SatBC::saturationDifference
double saturationDifference() const
Query a Periodic condition.
Definition BoundaryConditions.hpp:204
Opm::SteadyStateUpscaler::initImpl
virtual void initImpl(const Opm::parameter::ParameterGroup &param)
Override from superclass.
Definition SteadyStateUpscaler_impl.hpp:70
Opm::SteadyStateUpscaler::lastSaturationState
const std::vector< double > & lastSaturationState() const
Accessor for the steady state saturation field.
Definition SteadyStateUpscaler_impl.hpp:292
Opm::SteadyStateUpscaler::upscaleSteadyState
std::pair< permtensor_t, permtensor_t > upscaleSteadyState(const int flow_direction, const std::vector< double > &initial_saturation, const double boundary_saturation, const double pressure_drop, const permtensor_t &upscaled_perm)
Does a steady-state upscaling.
Definition SteadyStateUpscaler_impl.hpp:132
Opm::SteadyStateUpscaler::lastSaturationUpscaled
double lastSaturationUpscaled() const
Computes the upscaled saturation corresponding to the saturation field returned by lastSaturationStat...
Definition SteadyStateUpscaler_impl.hpp:301
Opm::SteadyStateUpscaler::SteadyStateUpscaler
SteadyStateUpscaler()
Default constructor.
Definition SteadyStateUpscaler_impl.hpp:53
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
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