OilfieldToolsNet/opm-upscaling
4
0
Fork 0
mirror of https://github.com/OPM/opm-upscaling.git synced 2025-02-25 18:45:23 -06:00
Code Issues Packages Projects Releases Wiki Activity

added upscale_relpermvisc.cpp - relperm upscaling in viscous

Browse Source 
limit. Not tested yet, committed due to structural changes in repo
This commit is contained in:
Kari B. Skjerve
2011-12-12 12:14:08 +01:00
parent 1810dfc754 ae6cd9723f
commit bb9c657d16
8 changed files with 194 additions and 129 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
.hgignore
View File

@@ -39,3 +39,4 @@ examples/upscale_*
examples/grdecldips
data
steadystate_test_implicit
steadystate_test_explicit

14
dune/upscaling/SteadyStateUpscalerImplicit.hpp
View File

@@ -83,7 +83,7 @@ namespace Dune
const std::vector<double>& initial_saturation,
const double boundary_saturation,
const double pressure_drop,
const permtensor_t& upscaled_perm);
const permtensor_t& upscaled_perm,bool& success);
/// Accessor for the steady state saturation field. This is empty until
/// upscaleSteadyState() is called, at which point it will
@@ -93,8 +93,11 @@ namespace Dune
/// Computes the upscaled saturation corresponding to the saturation field
/// returned by lastSaturationState(). Does this by computing total saturated
/// volume divided by total pore volume.
double lastSaturationUpscaled() const;
void setToCapillaryLimit(double average_s, std::vector<double>& s) const;
protected:
// ------- Typedefs -------
typedef typename Traits::template TransportSolver<GridInterface, typename Super::BCs>::Type TransportSolver;
@@ -114,10 +117,13 @@ namespace Dune
bool output_vtk_;
bool print_inoutflows_;
int simulation_steps_;
double stepsize_;
double init_stepsize_;
double relperm_threshold_;
double maximum_mobility_contrast_;
double sat_change_threshold_;
double sat_change_year_;
int max_it_;
double max_stepsize_;
double dt_sat_tol_;
TransportSolver transport_solver_;
};

240
dune/upscaling/SteadyStateUpscalerImplicit_impl.hpp
View File

@@ -41,7 +41,8 @@
#include <dune/porsol/common/MatrixInverse.hpp>
#include <dune/porsol/common/SimulatorUtilities.hpp>
#include <dune/porsol/common/ReservoirPropertyFixedMobility.hpp>
#include <dune/porsol/euler/MatchSaturatedVolumeFunctor.hpp>
#include <dune/common/Units.hpp>
#include <algorithm>
namespace Dune
@@ -51,14 +52,17 @@ namespace Dune
template <class Traits>
inline SteadyStateUpscalerImplicit<Traits>::SteadyStateUpscalerImplicit()
: Super(),
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)
: Super(),
output_vtk_(false),
print_inoutflows_(false),
simulation_steps_(10),
init_stepsize_(0.1),
relperm_threshold_(1.0e-8),
maximum_mobility_contrast_(1.0e9),
sat_change_year_(0.0),
max_it_(20),
max_stepsize_(10),
dt_sat_tol_(1e-2)
{
}
@@ -68,17 +72,19 @@ namespace Dune
template <class Traits>
inline void SteadyStateUpscalerImplicit<Traits>::initImpl(const parameter::ParameterGroup& param)
{
Super::initImpl(param);
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_ = Dune::unit::convert::from(param.getDefault("stepsize", stepsize_),
Dune::unit::day);
relperm_threshold_ = param.getDefault("relperm_threshold", relperm_threshold_);
Super::initImpl(param);
output_vtk_ = param.getDefault("output_vtk", output_vtk_);
print_inoutflows_ = param.getDefault("print_inoutflows", print_inoutflows_);
simulation_steps_ = param.getDefault("simulation_steps", simulation_steps_);
init_stepsize_ = Dune::unit::convert::from(param.getDefault("init_stepsize", init_stepsize_),
Dune::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);
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_ = Dune::unit::convert::from(param.getDefault("max_stepsize", max_stepsize_),Dune::unit::year);
transport_solver_.init(param);
// Set viscosities and densities if given.
double v1_default = this->res_prop_.viscosityFirstPhase();
double v2_default = this->res_prop_.viscositySecondPhase();
@@ -129,26 +135,33 @@ namespace Dune
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)
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);
SparseVector<double> injection(num_cells);
// Gravity.
FieldVector<double, 3> gravity(0.0);
// gravity[2] = -Dune::unit::gravity;
if (gravity.two_norm() > 0.0) {
MESSAGE("Warning: Gravity not yet handled by flow solver.");
}
static int count = 0;
++count;
int num_cells = this->ginterf_.numberOfCells();
// No source or sink.
std::vector<double> src(num_cells, 0.0);
SparseVector<double> injection(num_cells);
// Gravity.
FieldVector<double, 3> gravity(0.0);
// gravity[2] = -Dune::unit::gravity;
if (gravity.two_norm() > 0.0) {
MESSAGE("Warning: Gravity not yet handled by flow solver.");
}
// Set up initial saturation profile.
std::vector<double> saturation = initial_saturation;
for (int c = 0; c < int(saturation.size()); c++ ) {
double s_min = this->res_prop_.s_min(c);
double s_max = this->res_prop_.s_max(c);
saturation[c] = std::max(s_min+1e-4, saturation[c]);
saturation[c] = std::min(s_max-1e-4, saturation[c]);
}
// Set up boundary conditions.
setupUpscalingConditions(this->ginterf_, this->bctype_, flow_direction,
pressure_drop, boundary_saturation, this->twodim_hack_, this->bcond_);
@@ -167,60 +180,75 @@ namespace Dune
// Do a run till steady state. For now, we just do some pressure and transport steps...
std::vector<double> saturation_old = saturation;
bool stationary = false;
int it_count=0;
double stepsize=init_stepsize_;
for (int iter = 0; iter < simulation_steps_; ++iter) {
// Run transport solver.
transport_solver_.transportSolve(saturation, stepsize_, gravity, this->flow_solver_.getSolution(), injection);
std::vector<double> init_saturation(saturation);
while((!stationary) & (it_count < max_it_)){// & transport_cost < max_transport_cost_)
//for (int iter = 0; iter < simulation_steps_; ++iter) {
// Run transport solver.
std::cout << "Running transport step " << it_count <<" with stepsize " << stepsize/Dune::unit::year << " year \n";
bool converged=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_);
max_mod = this->flow_solver_.postProcessFluxes();
std::cout << "Max mod = " << max_mod << std::endl;
if(converged){
init_saturation=saturation;
/*
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;
// 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")
+ '-' + boost::lexical_cast<std::string>(count)
+ '-' + boost::lexical_cast<std::string>(flow_direction)
+ '-' + boost::lexical_cast<std::string>(it_count));
}
// 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]));
}
double ds_year=maxdiff*Dune::unit::year/stepsize;
std::cout << "Maximum saturation change/year: " << ds_year << 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;
if( ds_year < sat_change_year_){
stationary=true;
}
if(maxdiff< dt_sat_tol_){
stepsize=std::min(max_stepsize_,2*stepsize);
}
}else{
std::cerr << "Cutting time step\n";
init_saturation = saturation_old;
stepsize=stepsize/2.0;
}
// Output.
if (output_vtk_) {
writeVtkOutput(this->ginterf_,
this->res_prop_,
this->flow_solver_.getSolution(),
saturation,
std::string("output-steadystate")
+ '-' + boost::lexical_cast<std::string>(count)
+ '-' + boost::lexical_cast<std::string>(flow_direction)
+ '-' + boost::lexical_cast<std::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;
}
it_count+=1;
// Copy to old.
saturation_old = saturation;
}
success=stationary;
// Compute phase mobilities.
// First: compute maximal mobilities.
@@ -260,11 +288,11 @@ namespace Dune
// Set the steady state saturation fields for eventual outside access.
last_saturation_state_.swap(saturation);
// Compute the (anisotropic) upscaled mobilities.
// 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)));
// => 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
@@ -272,7 +300,7 @@ namespace Dune
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);
return std::make_pair(k_rw, k_ro);
}
@@ -282,7 +310,7 @@ namespace Dune
inline const std::vector<double>&
SteadyStateUpscalerImplicit<Traits>::lastSaturationState() const
{
return last_saturation_state_;
return last_saturation_state_;
}
@@ -304,22 +332,42 @@ namespace Dune
}
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);
std::vector<double> cap_press(num_cells, 0.0);
typedef typename UpscalerBase<Traits>::ResProp ResProp;
MatchSaturatedVolumeFunctor<GridInterface, ResProp> func(this->ginterf_, this->res_prop_, s_orig, cap_press);
double cap_press_range = 1e2;
double mod_low = 1e100;
double mod_high = -1e100;
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;
double mod_correct = modifiedRegulaFalsi(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
std::pair<double, double>& oil_inout,
const FlowSol& flow_solution,
const std::vector<double>& saturations) const
{
typedef typename GridInterface::CellIterator CellIter;
typedef typename CellIter::FaceIterator FaceIter;
double side1_flux = 0.0;
double side2_flux = 0.0;
double side1_flux_oil = 0.0;
double side2_flux_oil = 0.0;
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;
int num_cells = this->ginterf_.numberOfCells();
std::vector<double> cell_inflows_w(num_cells, 0.0);
@@ -358,7 +406,7 @@ namespace Dune
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;
// std::cout << "Inserted bid " << f->boundaryId() << std::endl;
}
cell_outflows_w[c->index()] += flux*frac_flow;
side2_flux += flux*frac_flow;
@@ -368,8 +416,8 @@ namespace Dune
}
}
}
water_inout = std::make_pair(side1_flux, side2_flux);
oil_inout = std::make_pair(side1_flux_oil, side2_flux_oil);
water_inout = std::make_pair(side1_flux, side2_flux);
oil_inout = std::make_pair(side1_flux_oil, side2_flux_oil);
}

34
dune/upscaling/SteadyStateUpscalerManagerImplicit.hpp
View File

@@ -97,7 +97,7 @@ namespace Dune
template<class Ostream, class Tensor>
void writeRelPerm(Ostream& os, const Tensor& K, double sat, double pdrop)
void writeRelPerm(Ostream& os, const Tensor& K, double sat, double pdrop,bool success)
{
/* const int num_rows = K.numRows(); */
/* const int num_cols = K.numCols(); */
@@ -124,8 +124,8 @@ namespace Dune
os << K(0,1) << '\t'; /* xy */
os << K(2,1) << '\t'; /* zy */
os << K(2,0) << '\t'; /* zx */
os << K(1,2); /* yz */
os << K(1,2) << '\t'; /* yz */
os << success;
os << std::endl;
@@ -176,6 +176,7 @@ namespace Dune
}
}
int flow_direction = param.getDefault("flow_direction", 0);
bool start_from_cl = param.getDefault("start_from_cl", false);
// Print the saturations and pressure drops.
// writeControl(std::cout, saturations, all_pdrops);
@@ -206,34 +207,43 @@ namespace Dune
kro_out << "# Pressuredrop Sw Krxx Kryy Krzz" << std::endl;
krw_out.precision(15); krw_out.setf(std::ios::scientific | std::ios::showpoint);
kro_out.precision(15); kro_out.setf(std::ios::scientific | std::ios::showpoint);
krw_out.precision(4); krw_out.setf(std::ios::scientific | std::ios::showpoint);
kro_out.precision(4); kro_out.setf(std::ios::scientific | std::ios::showpoint);
//#endif
// Then, compute some upscaled relative permeabilities.
int num_cells = upscaler.grid().size(0);
int num_sats = saturations.size();
//std::vector<bool> success_state;
for (int i = 0; i < num_sats; ++i) {
// Starting every computation with a trio of uniform profiles.
std::vector<double> init_sat(num_cells, saturations[i]);
std::vector<double> init_sat;
if (start_from_cl) {
upscaler.setToCapillaryLimit(saturations[i], init_sat);
} {
init_sat.resize(num_cells, saturations[i]);
}
const SparseTable<double>::row_type pdrops = all_pdrops[i];
int num_pdrops = pdrops.size();
for (int j = 0; j < num_pdrops; ++j) {
double pdrop = pdrops[j];
bool success = false;
std::pair<permtensor_t, permtensor_t> lambda
= upscaler.upscaleSteadyState(flow_direction, init_sat, saturations[i], pdrop, upscaled_K);
= upscaler.upscaleSteadyState(flow_direction, init_sat, saturations[i], pdrop, upscaled_K, success);
double usat = upscaler.lastSaturationUpscaled();
if(! success){
std::cout << "Upscaling failed for " << usat << std::endl;
}else{
init_sat = upscaler.lastSaturationState();
}
std::cout << "\n\nTensor of upscaled relperms for initial saturation " << saturations[i]
<< ", real steady-state saturation " << usat
<< " and pressure drop " << pdrop
<< ":\n\n[water]\n" << lambda.first
<< "\n[oil]\n" << lambda.second << std::endl;
// Changing initial saturations for next pressure drop to equal the steady state of the last
init_sat = upscaler.lastSaturationState();
writeRelPerm(krw_out, lambda.first , usat, pdrop);
writeRelPerm(kro_out, lambda.second, usat, pdrop);
writeRelPerm(krw_out, lambda.first , usat, pdrop, success);
writeRelPerm(kro_out, lambda.second, usat, pdrop, success);
}
}

20
dune/upscaling/UpscalerBase.hpp
View File

@@ -92,22 +92,22 @@ namespace Dune
/// Access the grid.
const GridType& grid() const;
/// Set boundary condition type. This may not be used to swicth
/// between Periodic and the other types, since the grid is
/// modified for Periodic conditions.
void setBoundaryConditionType(BoundaryConditionType type);
/// Set boundary condition type. This may not be used to swicth
/// between Periodic and the other types, since the grid is
/// modified for Periodic conditions.
void setBoundaryConditionType(BoundaryConditionType type);
/// Set the permeability of a cell directly. This will override
/// the permeability that was read from the eclipse file.
void setPermeability(const int cell_index, const permtensor_t& k);
/// Set the permeability of a cell directly. This will override
/// the permeability that was read from the eclipse file.
void setPermeability(const int cell_index, const permtensor_t& k);
/// Does a single-phase upscaling.
/// @return an upscaled permeability tensor.
permtensor_t upscaleSinglePhase();
/// Compute upscaled porosity.
/// @return total pore volume of all cells divided by total volume.
double upscalePorosity() const;
/// Compute upscaled porosity.
/// @return total pore volume of all cells divided by total volume.
double upscalePorosity() const;
protected:
// ------- Typedefs and enums -------

8
dune/upscaling/test/Makefile.am
View File

@@ -1,13 +1,13 @@
# $Date$
# $Revision$
check_PROGRAMS =
check_PROGRAMS =
noinst_PROGRAMS = \
aniso_implicit_steadystate_test \
aniso_steadystate_test \
implicit_steadystate_test \
opmpressure_test \
steadystate_test \
steadystate_test_explicit \
steadystate_test_implicit \
upscaling_test
@@ -32,8 +32,8 @@ implicit_steadystate_test.cpp
opmpressure_test_SOURCES = \
opmpressure_test.cpp
steadystate_test_SOURCES = \
steadystate_test.cpp
steadystate_test_explicit_SOURCES = \
steadystate_test_explicit.cpp
steadystate_test_implicit_SOURCES = \
steadystate_test_implicit.cpp

0
dune/upscaling/test/steadystate_test.cpp → dune/upscaling/test/steadystate_test_explicit.cpp
View File

6
examples/cpchop.cpp
View File

@@ -37,7 +37,7 @@
#include <fstream>
#include <iostream>
#include <cfloat>
#include <cmath>
int main(int argc, char** argv)
{
@@ -406,7 +406,7 @@ int main(int argc, char** argv)
largestSaturationInterval = SatDiff.second;
}
// Check for saneness of Ptestvalue:
if (isnan(Ptestvalue) | isinf(Ptestvalue)) {
if (std::isnan(Ptestvalue) || std::isinf(Ptestvalue)) {
std::cerr << "ERROR: Ptestvalue was inf or nan" << std::endl;
break; // Jump out of while-loop, just print out the results
// up to now and exit the program
@@ -440,7 +440,7 @@ int main(int argc, char** argv)
MonotCubicInterpolator CapPressureVsWaterSaturation(WaterSaturationVsCapPressure.get_fVector(),
WaterSaturationVsCapPressure.get_xVector());
std::vector<double> pcs;
for (uint satp=0; satp<nsatpoints; ++satp) {
for (int satp=0; satp<nsatpoints; ++satp) {
pcs.push_back(CapPressureVsWaterSaturation.evaluate(Swir+(Swor-Swir)/(nsatpoints-1)*satp));
}
pcvalues.push_back(pcs);