Loading...
Searching...
No Matches
SteadyStateUpscalerImplicit_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_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
#define OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
#include <boost/date_time/posix_time/posix_time.hpp>
#include <opm/porsol/common/MatrixInverse.hpp>
#include <opm/porsol/common/SimulatorUtilities.hpp>
#include <opm/porsol/common/ReservoirPropertyFixedMobility.hpp>
#include <opm/porsol/euler/MatchSaturatedVolumeFunctor.hpp>
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/upscaling/writeECLData.hpp>
#include <algorithm>
#include <iostream>
namespace Opm
{
template <class Traits>
inline SteadyStateUpscalerImplicit<Traits>::SteadyStateUpscalerImplicit()
: Super(),
use_gravity_(false),
output_vtk_(false),
output_ecl_(false),
print_inoutflows_(false),
simulation_steps_(10),
init_stepsize_(10),
relperm_threshold_(1.0e-8),
maximum_mobility_contrast_(1.0e9),
sat_change_year_(1.0e-5),
max_it_(100),
max_stepsize_(1e4),
dt_sat_tol_(1e-2),
use_maxdiff_(true)
{
}
template <class Traits>
inline void SteadyStateUpscalerImplicit<Traits>::initImpl(const Opm::ParameterGroup& param)
{
Super::initImpl(param);
use_gravity_ = param.getDefault("use_gravity", use_gravity_);
output_vtk_ = param.getDefault("output_vtk", output_vtk_);
output_ecl_ = param.getDefault("output_ecl", output_ecl_);
if (output_ecl_) {
grid_adapter_.init(Super::grid());
}
print_inoutflows_ = param.getDefault("print_inoutflows", print_inoutflows_);
simulation_steps_ = param.getDefault("simulation_steps", simulation_steps_);
init_stepsize_ = Opm::unit::convert::from(param.getDefault("init_stepsize", init_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_year_ = param.getDefault("sat_change_year", sat_change_year_);
dt_sat_tol_ = param.getDefault("dt_sat_tol", dt_sat_tol_);
max_it_ = param.getDefault("max_it", max_it_);
max_stepsize_ = Opm::unit::convert::from(param.getDefault("max_stepsize", max_stepsize_),Opm::unit::year);
use_maxdiff_ = param.getDefault("use_maxdiff", use_maxdiff_);
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);
}
}
// Workaround: duplicate of Opm::destripe in opm-simulators.
inline std::vector< double > destripe( const std::vector< double >& v,
size_t stride,
size_t offset )
{
std::vector< double > dst( v.size() / stride );
size_t di = 0;
for( size_t i = offset; i < v.size(); i += stride ) {
dst[ di++ ] = v[ i ];
}
return dst;
}
void waterSatToBothSat(const std::vector<double>& sw,
std::vector<double>& sboth)
{
int num = sw.size();
sboth.resize(2*num);
for (int i = 0; i < num; ++i) {
sboth[2*i] = sw[i];
sboth[2*i + 1] = 1.0 - sw[i];
}
}
} // anon namespace
template <class Traits>
inline std::pair<typename SteadyStateUpscalerImplicit<Traits>::permtensor_t,
typename SteadyStateUpscalerImplicit<Traits>::permtensor_t>
SteadyStateUpscalerImplicit<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,
bool& success)
{
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.");
}
// Put pore volume in vector.
std::vector<double> pore_vol;
pore_vol.reserve(num_cells);
double tot_pore_vol = 0.0;
typedef typename GridInterface::CellIterator CellIterT;
for (CellIterT c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
double cell_pore_vol = c->volume()*this->res_prop_.porosity(c->index());
pore_vol.push_back(cell_pore_vol);
tot_pore_vol += cell_pore_vol;
}
// 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_);
double max_mod = this->flow_solver_.postProcessFluxes();
std::cout << "Max mod = " << max_mod << std::endl;
// Do a run till steady state.
std::vector<double> saturation_old = saturation;
bool stationary = false;
int it_count = 0;
double stepsize = init_stepsize_;
double ecl_time = 0.0;
std::vector<double> ecl_sat;
std::vector<double> ecl_press;
while ((!stationary) && (it_count < max_it_)) { // && transport_cost < max_transport_cost_)
// Run transport solver.
std::cout << "Running transport step " << it_count << " with stepsize "
<< stepsize/Opm::unit::year << " years." << std::endl;
bool converged = transport_solver_.transportSolve(saturation, stepsize, gravity,
this->flow_solver_.getSolution(), injection);
// Run pressure solver.
if (converged) {
// this->flow_solver_.solve(this->res_prop_, saturation, this->bcond_, src,
// this->residual_tolerance_, this->linsolver_verbosity_, this->linsolver_type_);
// max_mod = this->flow_solver_.postProcessFluxes();
// std::cout << "Max mod of fluxes= " << max_mod << std::endl;
this->flow_solver_.postProcessFluxes();
// 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 " << it_count
<< "\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(it_count));
}
if (output_ecl_) {
const char* fd = "xyz";
std::string basename = std::string("ecldump-steadystate")
+ '-' + std::to_string(boundary_saturation)
+ '-' + std::to_string(fd[flow_direction])
+ '-' + std::to_string(pressure_drop);
waterSatToBothSat(saturation, ecl_sat);
getCellPressure(ecl_press, this->ginterf_, this->flow_solver_.getSolution());
data::Solution solution;
solution.insert( "PRESSURE" , UnitSystem::measure::pressure , ecl_press , data::TargetType::RESTART_SOLUTION );
solution.insert( "SWAT" , UnitSystem::measure::identity , destripe( ecl_sat, 2, 0) , data::TargetType::RESTART_SOLUTION );
ecl_time += stepsize;
boost::posix_time::ptime ecl_startdate( boost::gregorian::date(2012, 1, 1) );
boost::posix_time::ptime ecl_curdate = ecl_startdate + boost::posix_time::seconds(int(ecl_time));
boost::posix_time::ptime epoch( boost::gregorian::date( 1970, 1, 1 ) );
auto ecl_posix_time = ( ecl_curdate - epoch ).total_seconds();
const auto* cgrid = grid_adapter_.c_grid();
Opm::writeECLData(cgrid->cartdims[ 0 ],
cgrid->cartdims[ 1 ],
cgrid->cartdims[ 2 ],
cgrid->number_of_cells,
solution, it_count,
ecl_time, ecl_posix_time,
"./", basename);
}
// Comparing old to new.
double maxdiff = 0.0;
double euclidean_diff = 0.0;
for (int i = 0; i < num_cells; ++i) {
const double sat_diff_cell = saturation[i] - saturation_old[i];
maxdiff = std::max(maxdiff, std::fabs(sat_diff_cell));
euclidean_diff += sat_diff_cell * sat_diff_cell * pore_vol[i];
}
euclidean_diff = std::sqrt(euclidean_diff / tot_pore_vol);
double ds_year;
if (use_maxdiff_) {
ds_year = maxdiff*Opm::unit::year/stepsize;
std::cout << "Maximum saturation change/year: " << ds_year << std::endl;
if (maxdiff < dt_sat_tol_) {
stepsize=std::min(max_stepsize_,2*stepsize);
}
}
else {
ds_year = euclidean_diff*Opm::unit::year/stepsize;
std::cout << "Euclidean saturation change/year: " << ds_year << std::endl;
if (euclidean_diff < dt_sat_tol_) {
stepsize=std::min(max_stepsize_,2*stepsize);
}
}
if (ds_year < sat_change_year_) {
stationary = true;
}
} else {
std::cerr << "Cutting time step\n";
stepsize=stepsize/2.0;
}
++it_count;
// Copy to old.
saturation_old = saturation;
}
success = stationary;
// 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>&
SteadyStateUpscalerImplicit<Traits>::lastSaturationState() const
{
return last_saturation_state_;
}
template <class Traits>
double SteadyStateUpscalerImplicit<Traits>::lastSaturationUpscaled() const
{
typedef typename GridInterface::CellIterator CellIterT;
double pore_vol = 0.0;
double sat_vol = 0.0;
for (CellIterT 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;
}
template <class Traits>
void SteadyStateUpscalerImplicit<Traits>::initSatLimits(std::vector<double>& s) const
{
for (int c = 0; c < int (s.size()); c++ ) {
double s_min = this->res_prop_.s_min(c);
double s_max = this->res_prop_.s_max(c);
s[c] = std::max(s_min+1e-4, s[c]);
s[c] = std::min(s_max-1e-4, s[c]);
}
}
template <class Traits>
void SteadyStateUpscalerImplicit<Traits>::setToCapillaryLimit(double average_s, std::vector<double>& s) const
{
int num_cells = this->ginterf_.numberOfCells();
std::vector<double> s_orig(num_cells, average_s);
initSatLimits(s_orig);
std::vector<double> cap_press(num_cells, 0.0);
typedef typename UpscalerBase<Traits>::ResProp ResPropT;
MatchSaturatedVolumeFunctor<GridInterface, ResPropT> func(this->ginterf_, this->res_prop_, s_orig, cap_press);
double cap_press_range = 1e2;
double mod_low = 1e100;
double mod_high = -1e100;
Opm::bracketZero(func, 0.0, cap_press_range, mod_low, mod_high);
const int max_iter = 40;
const double nonlinear_tolerance = 1e-12;
int iterations_used = -1;
typedef Opm::RegulaFalsi<Opm::ThrowOnError> RootFinder;
double mod_correct = RootFinder::solve(func, mod_low, mod_high, max_iter, nonlinear_tolerance, iterations_used);
std::cout << "Moved capillary pressure solution by " << mod_correct << " after "
<< iterations_used << " iterations." << std::endl;
s = func.lastSaturations();
}
template <class Traits>
template <class FlowSol>
void SteadyStateUpscalerImplicit<Traits>::computeInOutFlows(std::pair<double, double>& water_inout,
std::pair<double, double>& oil_inout,
const FlowSol& flow_solution,
const std::vector<double>& saturations) const
{
typedef typename GridInterface::CellIterator CellIterT;
typedef typename CellIterT::FaceIterator FaceIterT;
double side1_flux = 0.0;
double side2_flux = 0.0;
double side1_flux_oil = 0.0;
double side2_flux_oil = 0.0;
std::map<int, double> frac_flow_by_bid;
// Two passes: First pass, deal with outflow, second pass, deal with inflow.
// This is for the periodic case, so that we are sure all fractional flows have
// been set in frac_flow_by_bid.
for (int pass = 0; pass < 2; ++pass) {
for (CellIterT c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
for (FaceIterT f = c->facebegin(); f != c->faceend(); ++f) {
if (f->boundary()) {
double flux = flow_solution.outflux(f);
const SatBC& sc = this->bcond_.satCond(f);
if (flux < 0.0 && pass == 1) {
// This is an inflow face.
double frac_flow = 0.0;
if (sc.isPeriodic()) {
assert(sc.saturationDifference() == 0.0);
int partner_bid = this->bcond_.getPeriodicPartner(f->boundaryId());
std::map<int, double>::const_iterator it = frac_flow_by_bid.find(partner_bid);
if (it == frac_flow_by_bid.end()) {
OPM_THROW(std::runtime_error,
"Could not find periodic partner fractional flow. "
"Face bid = " + std::to_string(f->boundaryId()) +
" and partner bid = " + std::to_string(partner_bid));
}
frac_flow = it->second;
} else {
assert(sc.isDirichlet());
frac_flow = this->res_prop_.fractionalFlow(c->index(), sc.saturation());
}
side1_flux += flux*frac_flow;
side1_flux_oil += flux*(1.0 - frac_flow);
} else if (flux >= 0.0 && pass == 0) {
// This is an outflow face.
double frac_flow = this->res_prop_.fractionalFlow(c->index(), saturations[c->index()]);
if (sc.isPeriodic()) {
frac_flow_by_bid[f->boundaryId()] = frac_flow;
// std::cout << "Inserted bid " << f->boundaryId() << std::endl;
}
side2_flux += flux*frac_flow;
side2_flux_oil += flux*(1.0 - frac_flow);
}
}
}
}
}
water_inout = std::make_pair(side1_flux, side2_flux);
oil_inout = std::make_pair(side1_flux_oil, side2_flux_oil);
}
} // namespace Opm
#endif // OPM_STEADYSTATEUPSCALERIMPLICIT_IMPL_HEADER
Opm::ReservoirPropertyFixedMobility
Definition ReservoirPropertyFixedMobility.hpp:46
Opm::SteadyStateUpscalerImplicit::initImpl
void initImpl(const Opm::ParameterGroup &param) override
Override from superclass.
Definition SteadyStateUpscalerImplicit_impl.hpp:78
Opm::SteadyStateUpscalerImplicit::initSatLimits
void initSatLimits(std::vector< double > &s) const
Ensure saturations are not outside table.
Definition SteadyStateUpscalerImplicit_impl.hpp:416
Opm::SteadyStateUpscalerImplicit::lastSaturationUpscaled
double lastSaturationUpscaled() const
Computes the upscaled saturation corresponding to the saturation field returned by lastSaturationStat...
Definition SteadyStateUpscalerImplicit_impl.hpp:400
Opm::SteadyStateUpscalerImplicit::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, bool &success)
Does a steady-state upscaling.
Definition SteadyStateUpscalerImplicit_impl.hpp:175
Opm::SteadyStateUpscalerImplicit::SteadyStateUpscalerImplicit
SteadyStateUpscalerImplicit()
Default constructor.
Definition SteadyStateUpscalerImplicit_impl.hpp:56
Opm::SteadyStateUpscalerImplicit::lastSaturationState
const std::vector< double > & lastSaturationState() const
Accessor for the steady state saturation field.
Definition SteadyStateUpscalerImplicit_impl.hpp:391
Opm::UpscalerBase::grid
const GridType & grid() const
Access the grid.
Definition UpscalerBase_impl.hpp:181
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
Opm::getCellPressure
void getCellPressure(std::vector< double > &cell_pressure, const GridInterface &ginterf, const FlowSol &flow_solution)
Definition SimulatorUtilities.hpp:182
Opm::MatchSaturatedVolumeFunctor
Definition MatchSaturatedVolumeFunctor.hpp:68