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

Move upscaling code to separate module .

Browse Source
This commit is contained in:
Bård Skaflestad
2010-08-23 09:34:14 +00:00
parent cef76ee3e6
commit 2a2f2a95fa
17 changed files with 2223 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
dune/upscaling/Makefile.am
View File

@@ -1,4 +1,10 @@
upscalingincludedir = $(includedir)/dune/upscaling
upscalinginclude_HEADERS = upscaling.hh
# $Date: 2010-04-29 11:11:10 +0200 (Thu, 29 Apr 2010) $
# $Revision: 836 $
SUBDIRS = test
upscalingdir = $(includedir)/dune/upscaling
upscaling_HEADERS = SinglePhaseUpscaler.hpp \
SteadyStateUpscaler.hpp SteadyStateUpscaler_impl.hpp
include $(top_srcdir)/am/global-rules

61
dune/upscaling/SinglePhaseUpscaler.hpp Normal file
View File

@@ -0,0 +1,61 @@
//===========================================================================
//
// File: SinglePhaseUpscaler.hpp
//
// Created: Fri Aug 28 13:38:54 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd 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 OPENRS_SINGLEPHASEUPSCALER_HEADER
#define OPENRS_SINGLEPHASEUPSCALER_HEADER
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/upscaling/UpscalingTraits.hpp>
#include <dune/upscaling/UpscalerBase.hpp>
namespace Dune
{
/**
@brief A class for doing single phase (permeability) upscaling.
@author Atgeirr F. Rasmussen <atgeirr@sintef.no>
*/
class SinglePhaseUpscaler : public UpscalerBase<UpscalingTraitsBasic>
{
};
} // namespace Dune
#endif // OPENRS_SINGLEPHASEUPSCALER_HEADER

127
dune/upscaling/SteadyStateUpscaler.hpp Normal file
View File

@@ -0,0 +1,127 @@
//===========================================================================
//
// File: SteadyStateUpscaler.hpp
//
// Created: Fri Aug 28 14:01:19 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd 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 OPENRS_STEADYSTATEUPSCALER_HEADER
#define OPENRS_STEADYSTATEUPSCALER_HEADER
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/upscaling/UpscalerBase.hpp>
#include <dune/solvers/euler/EulerUpstream.hpp>
#include <dune/solvers/euler/ImplicitCapillarity.hpp>
#include <boost/array.hpp>
namespace Dune
{
/**
@brief A class for doing steady state upscaling.
@author Atgeirr F. Rasmussen <atgeirr@sintef.no>
*/
template <class Traits>
class SteadyStateUpscaler : public UpscalerBase<Traits>
{
public:
// ------- Typedefs and enums -------
typedef UpscalerBase<Traits> Super;
typedef typename Super::permtensor_t permtensor_t;
typedef typename UpscalerBase<Traits>::GridInterface GridInterface;
enum { Dimension = UpscalerBase<Traits>::Dimension };
// ------- Methods -------
/// Default constructor.
SteadyStateUpscaler();
/// Does a steady-state upscaling.
/// @param flow_direction The cardinal direction in which to impose a pressure gradient for the purpose of converging to steady state.
/// @param initial_saturation the initial saturation profile for the steady-state computation.
/// The vector must have size equal to the number of cells in the grid.
/// @param boundary_saturation the saturation of fluid flowing in across the boundary,
/// only needed for nonperiodic upscaling.
/// @param pressure_drop the pressure drop in Pascal over the domain.
/// @param upscaled_perm typically the output of upscaleSinglePhase().
/// @return the upscaled relative permeability matrices of both phases.
/// The relative permeability matrix, call it k, is such that if K_w is the phase
/// permeability and K the absolute permeability, K_w = k*K.
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);
/// Accessor for the steady state saturation field. This is empty until
/// upscaleSteadyState() is called, at which point it will
/// contain the last computed (steady) saturation state.
const std::vector<double>& lastSaturationState() const;
/// 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;
protected:
// ------- Typedefs -------
typedef typename Traits::template TransportSolver<GridInterface, typename Super::BCs>::Type TransportSolver;
// ------- Methods -------
template <class FlowSol>
void computeInOutFlows(std::pair<double, double>& water_inout,
std::pair<double, double>& oil_inout,
const FlowSol& flow_solution,
const std::vector<double>& saturations) const;
/// Override from superclass.
virtual void initImpl(const parameter::ParameterGroup& param);
// ------- Data members -------
std::vector<double> last_saturation_state_;
bool output_vtk_;
bool print_inoutflows_;
int simulation_steps_;
double stepsize_;
double relperm_threshold_;
double sat_change_threshold_;
TransportSolver transport_solver_;
};
} // namespace Dune
#include "SteadyStateUpscaler_impl.hpp"
#endif // OPENRS_STEADYSTATEUPSCALER_HEADER

226
dune/upscaling/SteadyStateUpscalerManager.hpp Normal file
View File

@@ -0,0 +1,226 @@
//===========================================================================
//
// File: SteadyStateUpscalerManager.hpp
//
// Created: Thu Apr 29 11:12:56 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Jostein R Natvig <jostein.r.natvig@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 OPENRS_STEADYSTATEUPSCALERMANAGER_HEADER
#define OPENRS_STEADYSTATEUPSCALERMANAGER_HEADER
#include <dune/upscaling/SteadyStateUpscaler.hpp>
#include <dune/upscaling/UpscalingTraits.hpp>
#include <dune/common/Units.hpp>
#include <dune/common/SparseTable.hpp>
#include <cmath>
#include <fstream>
#include <iostream>
namespace Dune
{
/// Reads saturation and pressure drop data from an input stream.
template <class Istream>
void readControl(Istream& is, std::vector<double>& saturations, SparseTable<double>& all_pdrops)
{
int num_sats;
is >> num_sats;
std::vector<double> sat(num_sats);
SparseTable<double> ap;
std::vector<double> tmp_pd;
for (int i = 0; i < num_sats; ++i) {
is >> sat[i];
int num_pdrops;
is >> num_pdrops;
tmp_pd.resize(num_pdrops);
for (int j = 0; j < num_pdrops; ++j) {
is >> tmp_pd[j];
}
ap.appendRow(tmp_pd.begin(), tmp_pd.end());
}
all_pdrops.swap(ap);
saturations.swap(sat);
}
/// Writes saturation and pressure drop data to an output stream.
template <class Ostream>
void writeControl(Ostream& os, const std::vector<double>& saturations, const SparseTable<double>& all_pdrops)
{
int num_sats = saturations.size();
os << num_sats << '\n';
for (int i = 0; i < num_sats; ++i) {
os << saturations[i] << '\n';
int num_pdrops = all_pdrops[i].size();
os << num_pdrops;
for (int j = 0; j < num_pdrops; ++j) {
os << ' ' << all_pdrops[i][j];
}
os << '\n';
}
}
template<class Ostream, class Tensor>
void writeRelPerm(Ostream& os, const Tensor& K, double sat, double pdrop)
{
const int num_rows = K.numRows();
const int num_cols = K.numCols();
os << sat << ' ';
for (int i = 0; i < num_rows; ++i) {
for (int j = 0; j < num_cols; ++j) {
os << K(i,j) << ' ';
}
}
os << pdrop << std::endl;
}
template <class Traits>
class SteadyStateUpscalerManager
{
public:
void upscale(const parameter::ParameterGroup& param)
{
// Control structure.
std::vector<double> saturations;
SparseTable<double> all_pdrops;
bool from_file = param.has("sat_pdrop_filename");
if (from_file) {
std::string filename = param.get<std::string>("sat_pdrop_filename");
std::ifstream file(filename.c_str());
if (!file) {
THROW("Could not open file " << filename);
}
readControl(file, saturations, all_pdrops);
} else {
// Get a linear range of saturations.
int num_sats = param.getDefault("num_sats", 4);
double min_sat = param.getDefault("min_sat", 0.2);
double max_sat = param.getDefault("max_sat", 0.8);
saturations.resize(num_sats);
for (int i = 0; i < num_sats; ++i) {
double factor = num_sats == 1 ? 0 : double(i)/double(num_sats - 1);
saturations[i] = (1.0 - factor)*min_sat + factor*max_sat;
}
// Get a logarithmic range of pressure drops.
int num_pdrops = param.getDefault("num_pdrops", 5);
double log_min_pdrop = std::log(param.getDefault("min_pdrop", 1e2));
double log_max_pdrop = std::log(param.getDefault("max_pdrop", 1e6));
std::vector<double> pdrops;
pdrops.resize(num_pdrops);
for (int i = 0; i < num_pdrops; ++i) {
double factor = num_pdrops == 1 ? 0 : double(i)/double(num_pdrops - 1);
pdrops[i] = std::exp((1.0 - factor)*log_min_pdrop + factor*log_max_pdrop);
}
// Assign the same pressure drops to all saturations.
for (int i = 0; i < num_sats; ++i) {
all_pdrops.appendRow(pdrops.begin(), pdrops.end());
}
}
int flow_direction = param.getDefault("flow_direction", 0);
// Print the saturations and pressure drops.
// writeControl(std::cout, saturations, all_pdrops);
// Initialize upscaler.
typedef SteadyStateUpscaler<Traits> Upscaler;
typedef typename Upscaler::permtensor_t permtensor_t;
Upscaler upscaler;
upscaler.init(param);
// First, compute an upscaled permeability.
permtensor_t upscaled_K = upscaler.upscaleSinglePhase();
permtensor_t upscaled_K_copy = upscaled_K;
upscaled_K_copy *= (1.0/(prefix::milli*unit::darcy));
std::cout.precision(15);
std::cout << "Upscaled K in millidarcy:\n" << upscaled_K_copy << std::endl;
std::cout << "Upscaled porosity: " << upscaler.upscalePorosity() << std::endl;
#if WRITE_RELPERM_TO_FILE
// Create output streams for upscaled relative permeabilities
std::string kr_filename = param.getDefault<std::string>("kr_filename", "upscaled_relperm");
std::string krw_filename = kr_filename + "_wat";
std::string kro_filename = kr_filename + "_oil";
std::ofstream krw_out(krw_filename.c_str());
std::ofstream kro_out(kro_filename.c_str());
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);
#endif
// Then, compute some upscaled relative permeabilities.
int num_cells = upscaler.grid().size(0);
int num_sats = saturations.size();
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]);
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];
std::pair<permtensor_t, permtensor_t> lambda
= upscaler.upscaleSteadyState(flow_direction, init_sat, saturations[i], pdrop, upscaled_K);
double usat = upscaler.lastSaturationUpscaled();
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();
#if WRITE_RELPERM_TO_FILE
writeRelPerm(krw_out, lambda.first , usat, pdrop);
writeRelPerm(kro_out, lambda.second, usat, pdrop);
#endif
}
}
}
};
} // namespace Dune
#endif // OPENRS_STEADYSTATEUPSCALERMANAGER_HEADER

348
dune/upscaling/SteadyStateUpscaler_impl.hpp Normal file
View File

@@ -0,0 +1,348 @@
//===========================================================================
//
// File: SteadyStateUpscaler_impl.hpp
//
// Created: Fri Aug 28 14:07:51 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd 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 OPENRS_STEADYSTATEUPSCALER_IMPL_HEADER
#define OPENRS_STEADYSTATEUPSCALER_IMPL_HEADER
#include <boost/lexical_cast.hpp>
#include <dune/solvers/common/MatrixInverse.hpp>
#include <dune/solvers/common/SimulatorUtilities.hpp>
#include <dune/solvers/common/ReservoirPropertyFixedMobility.hpp>
#include <algorithm>
namespace Dune
{
template <class Traits>
inline SteadyStateUpscaler<Traits>::SteadyStateUpscaler()
: Super(),
output_vtk_(false),
print_inoutflows_(false),
simulation_steps_(10),
stepsize_(0.1),
relperm_threshold_(1.0e-4),
sat_change_threshold_(0.0)
{
}
template <class Traits>
inline void SteadyStateUpscaler<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_);
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 {
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);
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;
// 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. 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_);
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")
+ '-' + 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;
}
// Copy to old.
saturation_old = saturation;
}
// Compute phase mobilities.
typedef typename Super::ResProp::Mobility Mob;
std::vector<Mob> mob1(num_cells);
std::vector<Mob> mob2(num_cells);
const double mob1_threshold = relperm_threshold_ / this->res_prop_.viscosityFirstPhase();
const double mob2_threshold = relperm_threshold_ / this->res_prop_.viscositySecondPhase();
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;
}
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
{
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;
std::map<int, double> frac_flow_by_bid;
int num_cells = this->ginterf_.numberOfCells();
std::vector<double> cell_inflows_w(num_cells, 0.0);
std::vector<double> cell_outflows_w(num_cells, 0.0);
// 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 (CellIter c = this->ginterf_.cellbegin(); c != this->ginterf_.cellend(); ++c) {
for (FaceIter 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()) {
THROW("Could not find periodic partner fractional flow. Face bid = " << f->boundaryId()
<< " and partner bid = " << partner_bid);
}
frac_flow = it->second;
} else {
ASSERT(sc.isDirichlet());
frac_flow = this->res_prop_.fractionalFlow(c->index(), sc.saturation());
}
cell_inflows_w[c->index()] += flux*frac_flow;
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;
}
cell_outflows_w[c->index()] += flux*frac_flow;
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 Dune
#endif // OPENRS_STEADYSTATEUPSCALER_IMPL_HEADER

153
dune/upscaling/UpscalerBase.hpp Normal file
View File

@@ -0,0 +1,153 @@
//===========================================================================
//
// File: UpscalerBase.hpp
//
// Created: Thu Apr 29 10:20:22 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 OPENRS_UPSCALERBASE_HEADER
#define OPENRS_UPSCALERBASE_HEADER
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/param/ParameterGroup.hpp>
#include <dune/grid/CpGrid.hpp>
#include <dune/common/EclipseGridParser.hpp>
#include <dune/solvers/common/GridInterfaceEuler.hpp>
#include <dune/solvers/common/BoundaryConditions.hpp>
namespace Dune
{
/**
@brief A base class for upscaling.
@author Atgeirr F. Rasmussen <atgeirr@sintef.no>
*/
template <class Traits>
class UpscalerBase
{
protected:
public:
// ------- Typedefs -------
typedef CpGrid GridType;
enum { Dimension = GridType::dimension };
typedef GridInterfaceEuler<GridType> GridInterface;
typedef typename Traits::template ResProp<Dimension>::Type ResProp;
/// A type for the upscaled permeability.
typedef typename ResProp::MutablePermTensor permtensor_t;
enum BoundaryConditionType { Fixed = 0, Linear = 1, Periodic = 2 };
// ------- Methods -------
/// Default constructor.
UpscalerBase();
/// Initializes the upscaler from parameters.
void init(const parameter::ParameterGroup& param);
/// Initializes the upscaler from given arguments.
void init(const EclipseGridParser& parser,
BoundaryConditionType bctype,
double perm_threshold,
double z_tolerance = 0.0,
double residual_tolerance = 1e-8,
int linsolver_verbosity = 0,
int linsolver_type = 1,
bool twodim_hack = false);
/// 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 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;
protected:
// ------- Typedefs and enums -------
typedef GridInterface::CellIterator CellIter;
typedef CellIter::FaceIterator FaceIter;
typedef BasicBoundaryConditions<true, true> BCs;
typedef typename Traits::template FlowSolver<GridInterface, BCs>::Type FlowSolver;
// ------- Methods -------
template <class FlowSol>
double computeAverageVelocity(const FlowSol& flow_solution,
const int flow_dir,
const int pdrop_dir) const;
double computeDelta(const int flow_dir) const;
template <class FluidInterface>
permtensor_t upscaleEffectivePerm(const FluidInterface& fluid);
virtual void initImpl(const parameter::ParameterGroup& param);
virtual void initFinal(const parameter::ParameterGroup& param);
// ------- Data members -------
BoundaryConditionType bctype_;
bool twodim_hack_;
double residual_tolerance_;
int linsolver_verbosity_;
int linsolver_type_;
GridType grid_;
GridInterface ginterf_;
ResProp res_prop_;
BCs bcond_;
FlowSolver flow_solver_;
};
} // namespace Dune
#include "UpscalerBase_impl.hpp"
#endif // OPENRS_UPSCALERBASE_HEADER

373
dune/upscaling/UpscalerBase_impl.hpp Normal file
View File

@@ -0,0 +1,373 @@
//===========================================================================
//
// 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 OPENRS_UPSCALERBASE_IMPL_HEADER
#define OPENRS_UPSCALERBASE_IMPL_HEADER
#include <dune/solvers/common/setupGridAndProps.hpp>
#include <dune/solvers/common/setupBoundaryConditions.hpp>
#include <dune/solvers/common/ReservoirPropertyTracerFluid.hpp>
namespace Dune
{
template <class Traits>
inline UpscalerBase<Traits>::UpscalerBase()
: bctype_(Fixed),
twodim_hack_(false),
residual_tolerance_(1e-8),
linsolver_verbosity_(0),
linsolver_type_(1)
{
}
template <class Traits>
inline void UpscalerBase<Traits>::init(const parameter::ParameterGroup& param)
{
initImpl(param);
initFinal(param);
}
template <class Traits>
inline void UpscalerBase<Traits>::initImpl(const parameter::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_);
// Ensure sufficient grid support for requested boundary
// condition type.
//
parameter::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 (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 parameter::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 EclipseGridParser& parser,
BoundaryConditionType bctype,
double perm_threshold,
double z_tolerance,
double residual_tolerance,
int linsolver_verbosity,
int linsolver_type,
bool twodim_hack)
{
bctype_ = bctype;
residual_tolerance_ = residual_tolerance;
linsolver_verbosity_ = linsolver_verbosity;
linsolver_type_ = linsolver_type;
twodim_hack_ = twodim_hack;
// Faking some parameters depending on bc type.
bool periodic_ext = (bctype_ == Periodic);
bool turn_normals = false;
bool clip_z = (bctype_ == Periodic);
bool unique_bids = (bctype_ == Linear || bctype_ == Periodic);
std::string rock_list("no_list");
setupGridAndPropsEclipse(parser, z_tolerance,
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)) {
THROW("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.
FieldVector<double, 3> gravity(0.0);
// gravity[2] = -Dune::unit::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.
flow_solver_.solve(fluid, sat, bcond_, src, residual_tolerance_, linsolver_verbosity_, linsolver_type_);
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:
THROW("Unknown boundary type: " << 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;
int num_faces = 0;
int num_bdyfaces = 0;
int num_side1 = 0;
int num_side2 = 0;
for (CellIter c = ginterf_.cellbegin(); c != ginterf_.cellend(); ++c) {
for (FaceIter f = c->facebegin(); f != c->faceend(); ++f) {
++num_faces;
if (f->boundary()) {
++num_bdyfaces;
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) {
++num_side1;
if (flow_dir == pdrop_dir && flux > 0.0) {
std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()
<< " Magnitude: " << std::fabs(flux) << std::endl;
// THROW("Detected outflow at entry face: " << face);
}
side1_flux += flux*norm_comp;
side1_area += area;
} else {
ASSERT(canon_bid - 1 == 2*flow_dir + 1);
++num_side2;
if (flow_dir == pdrop_dir && flux < 0.0) {
std::cerr << "Flow may be in wrong direction at bid: " << f->boundaryId()
<< " Magnitude: " << std::fabs(flux) << std::endl;
// THROW("Detected inflow at exit face: " << face);
}
side2_flux += flux*norm_comp;
side2_area += area;
}
}
}
}
}
// std::cout << "Faces: " << num_faces << " Boundary faces: " << num_bdyfaces
// << " Side 1 faces: " << num_side1 << " Side 2 faces: " << num_side2 << std::endl;
// 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;
}
} // namespace Dune
#endif // OPENRS_UPSCALERBASE_IMPL_HEADER

49
dune/upscaling/UpscalingTraits.hpp Normal file
View File

@@ -0,0 +1,49 @@
//===========================================================================
//
// File: UpscalingTraits.hpp
//
// Created: Wed Apr 28 10:36:42 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 OPENRS_UPSCALINGTRAITS_HEADER
#define OPENRS_UPSCALINGTRAITS_HEADER
#include <dune/solvers/common/SimulatorTraits.hpp>
namespace Dune
{
typedef SimulatorTraits<Isotropic, Explicit> UpscalingTraitsBasic;
typedef SimulatorTraits<Anisotropic, Explicit> UpscalingTraitsAnisoRelperm;
} // namespace Dune
#endif // OPENRS_UPSCALINGTRAITS_HEADER

130
dune/upscaling/periodic_testresults.txt Normal file
View File

@@ -0,0 +1,130 @@
Old old code (the code first delivered december 2006).
=====================================================
Command:
--------
atgeirr@elrond:~/gotools/mimetic_statoil$ time app/upscale_test ~/workdir/data/opm_data/grids/wave_zlim.grdecl 2
Computing tensor using PERIODIC boundary conditions.
--------
1) EPS = 1e-16
154.083438502815 0.0892604992569009 0.212772595133234
0.146359826095459 156.999934343510 0.0177304108534503
0.213914312649649 0.00862349720761078 28.8258061024766
real 0m49.700s
user 0m46.868s
sys 0m2.480s
2) EPS = 1e-10
154.083438504296 0.0892604955824901 0.212772484789619
0.146359823701169 156.999934345602 0.0177303899771650
0.213914312547322 0.00862349620390634 28.8258060039328
real 0m42.343s
user 0m39.686s
sys 0m2.417s
Old code (delivered december 2008).
==================================
Command:
-------
atgeirr@elrond:~/workdir/trunk/zzz_work_here/upscaling/test_fall_2009$ time ../../../upscaling/upscale_singlephase config.xml
-------
[ 154.269327967824 0.086113779584908 0.214772628024014 |
| 0.0861138229423957 157.499004400429 0.00961100572800983 |
| 0.214773290593839 0.00961022430314095 28.8157255228606 ]
real 2m8.481s
user 1m16.656s
sys 0m51.289s
New code (the dune-cornerpoint module).
======================================
Command:
-------
atgeirr@elrond:~/dune_release/dune-cornerpoint/upscaling/test$ ./upscaling_test fileformat=eclipse filename=~/workdir/data/opm_data/grids/wave_zlim.grdecl use_unique_boundary_ids=true boundary_condition_type=2 periodic_extension=true residual_tolerance=1e-10
-------
1) Debug mode (all others are in release mode), residual_tolerance = 1e-10
[done to verify results with optimization, result identical => ok]
154.14299689765 -0.0951137484162823 0.213990389076227
-0.0951137479212948 157.371580359195 0.00652322307363517
0.213990127124522 0.00652273083131193 28.4615206338715
2) residual_tolerance = 1e-16
154.142996896361 -0.09511374418474 0.213990133897884
-0.0951137442634729 157.371580348364 0.00652274146037189
0.213990133899447 0.0065227414798937 28.4615196615018
real 1m24.065s
user 1m19.297s
sys 0m4.258s
3) residual_tolerance = 1e-10
154.14299689765 -0.0951137484162823 0.213990389076227
-0.0951137479212948 157.371580359195 0.00652322307363517
0.213990127124522 0.00652273083131193 28.4615206338715
real 0m58.475s
user 0m54.116s
sys 0m4.134s
3.5) residual_tolerance = 1e-8
154.142996830925 -0.0951135283476298 0.213993638248287
-0.0951134787818674 157.371579632703 0.00651798355914985
0.213991129269592 0.00652141149859271 28.4615086537167
real 0m47.659s
user 0m42.728s
sys 0m4.113s
4) No linear solver.
real 0m11.577s
user 0m7.962s
sys 0m3.551s
5) boundary_condition_type = 0, periodic_extension = true
[after bugfix, we have a difference due to extension, but no longer a
large one]
154.746173362754 0 0
0 158.269808171916 0
0 0 29.0714824825093
real 0m46.266s
user 0m42.007s
sys 0m3.867s
6) boundary_condition_type = 0, periodic_extension = false
154.750896840188 0 0
0 158.27098254722 0
0 0 29.0774483359038
real 0m41.786s
user 0m37.747s
sys 0m3.847s

12
dune/upscaling/solver_tests.txt Normal file
View File

@@ -0,0 +1,12 @@
SYMMETRIC = 1 for all cases
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 0 0 1 (== SPE10 setting) : 90, 91, 61 iters, 46783 level 1 z-system
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 0 1 1 : 26, 27, 16 iters, 46783 level 1 z-system
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 1 1 0 (== old setting) : 19, 19, 9 iters, 62497 level 1 z-system
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 1 0 0 : 80, 89, 61 iters, 62497 level 1 z-system
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 0 1 0 : 21, 21, 11 iters, 74787 level 1 z-system
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 0 0 0 : 85, 84, 58 iters, 74787 level 1 z-system
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 1 0 1 : 92, 93, 61 iters, 38486 level 1 z-system
FIRST_DIAGONAL, SMOOTHER_ILU, ANISOTROPIC_3D = 1 1 1 : 49, 41, 17 iters, 38486 level 1 z-system

31
dune/upscaling/test/Makefile.am Normal file
View File

@@ -0,0 +1,31 @@
# $Date$
# $Revision$
check_PROGRAMS =
noinst_PROGRAMS = upscaling_test steadystate_test aniso_steadystate_test implicit_steadystate_test aniso_implicit_steadystate_test #upscale_perm upscale_relperm
AM_CPPFLAGS += $(DUNEMPICPPFLAGS) $(BOOST_CPPFLAGS) $(SUPERLU_CPPFLAGS)
AM_LDFLAGS += $(DUNEMPILDFLAGS) $(BOOST_LDFLAGS) $(SUPERLU_LDFLAGS)
LDADD = $(DUNE_LIBS) $(DUNEMPILIBS) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
$(BOOST_FILESYSTEM_LIB) $(BOOST_SYSTEM_LIB) \
../../lib/libdunecornerpoint.la \
$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SUPERLU_LIBS)
upscaling_test_SOURCES = upscaling_test.cpp
steadystate_test_SOURCES = steadystate_test.cpp
aniso_steadystate_test_SOURCES = aniso_steadystate_test.cpp
implicit_steadystate_test_SOURCES = implicit_steadystate_test.cpp
aniso_implicit_steadystate_test_SOURCES = aniso_implicit_steadystate_test.cpp
#upscale_perm_SOURCES = upscale_perm.C
#upscale_relperm_SOURCES = upscale_relperm.C
TESTS = $(check_PROGRAMS)
include $(top_srcdir)/am/global-rules

49
dune/upscaling/test/aniso_implicit_steadystate_test.cpp Normal file
View File

@@ -0,0 +1,49 @@
//===========================================================================
//
// File: aniso_implicit_steadystate_test.cpp
//
// Created: Thu Jul 22 15:25:07 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd Skaflestad <bard.skaflestad@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/>.
*/
#define VERBOSE
#include <dune/upscaling/SteadyStateUpscalerManager.hpp>
#include <dune/upscaling/UpscalingTraits.hpp>
using namespace Dune;
int main(int argc, char** argv)
{
// Initialize.
parameter::ParameterGroup param(argc, argv);
SteadyStateUpscalerManager<SimulatorTraits<Anisotropic, ImplicitCap> > mgr;
mgr.upscale(param);
}

51
dune/upscaling/test/aniso_steadystate_test.cpp Normal file
View File

@@ -0,0 +1,51 @@
//===========================================================================
//
// File: aniso_steadystate_test.cpp
//
// Created: Thu Apr 29 10:29:35 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// Jostein R Natvig <jostein.r.natvig@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/>.
*/
//#define VERBOSE
#include <dune/upscaling/SteadyStateUpscalerManager.hpp>
#include <dune/upscaling/UpscalingTraits.hpp>
using namespace Dune;
int main(int argc, char** argv)
{
// Initialize.
parameter::ParameterGroup param(argc, argv);
// MPIHelper::instance(argc,argv);
SteadyStateUpscalerManager<UpscalingTraitsAnisoRelperm> mgr;
mgr.upscale(param);
}

49
dune/upscaling/test/implicit_steadystate_test.cpp Normal file
View File

@@ -0,0 +1,49 @@
//===========================================================================
//
// File: implicit_steadystate_test.cpp
//
// Created: Thu Jul 22 15:36:52 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd Skaflestad <bard.skaflestad@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/>.
*/
#define VERBOSE
#include <dune/upscaling/SteadyStateUpscalerManager.hpp>
#include <dune/upscaling/UpscalingTraits.hpp>
using namespace Dune;
int main(int argc, char** argv)
{
// Initialize.
parameter::ParameterGroup param(argc, argv);
SteadyStateUpscalerManager<SimulatorTraits<Isotropic, ImplicitCap> > mgr;
mgr.upscale(param);
}

52
dune/upscaling/test/steadystate_test.cpp Normal file
View File

@@ -0,0 +1,52 @@
//===========================================================================
//
// File: steadystate_test.cpp
//
// Created: Fri Aug 28 14:11:03 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd 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/>.
*/
//#define VERBOSE
#include <dune/upscaling/SteadyStateUpscalerManager.hpp>
#include <dune/upscaling/UpscalingTraits.hpp>
using namespace Dune;
int main(int argc, char** argv)
{
// Initialize.
parameter::ParameterGroup param(argc, argv);
// MPIHelper::instance(argc,argv);
SteadyStateUpscalerManager<UpscalingTraitsBasic> mgr;
mgr.upscale(param);
}

57
dune/upscaling/test/upscaling_test.cpp Normal file
View File

@@ -0,0 +1,57 @@
//===========================================================================
//
// File: upscaling_test.cpp
//
// Created: Mon Aug 10 08:21:44 2009
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
// B<>rd 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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/upscaling/SinglePhaseUpscaler.hpp>
#include <dune/common/Units.hpp>
using namespace Dune;
using namespace Dune::prefix;
using namespace Dune::unit;
int main(int argc, char** argv)
{
parameter::ParameterGroup param(argc, argv);
// MPIHelper::instance(argc,argv);
SinglePhaseUpscaler upscaler;
upscaler.init(param);
SinglePhaseUpscaler::permtensor_t upscaled_K = upscaler.upscaleSinglePhase();
upscaled_K *= (1.0/(milli*darcy));
std::cout.precision(15);
std::cout << "Upscaled K in millidarcy:\n" << upscaled_K << std::endl;
}

447
dune/upscaling/testresults.txt Normal file
View File

@@ -0,0 +1,447 @@
All results are from my laptop. Your results may differ.
Conclusions first: New code (without running the linear solver) is
faster, by almost 4 seconds. When actually running the solver, the new
code keeps its speed edge until the residual tolerance becomes low
(less than 1e-10). This means that the Dune AMG is a little slower
than SAMG, which is not all that surprising. The surprise is that it
(mostly) keeps pace with SAMG, which was not quite expected.
The codes seem to converge to different solutions, which is
unfortunate, but explainable: The codes use
a) different approximations to areas, volumes and normals,
b) different linear systems (related to handling of bdy cond).
The first difference is zero for cartesian grids, testing with a
cartesian grid and heterogenous permeability should shed some light on
this.
The SAMG solver seems to give a little more precision (compared to its
own 1e-16 solution) than the Dune solver. I offer no explanation, but
do not think this necessarily means anything problematic. Will keep an
eye on this, though.
--Atgeirr
Old old code (the code first delivered december 2006).
=====================================================
Command:
--------
atgeirr@elrond:~/gotools/mimetic_statoil$ time app/upscale_test ~/workdir/data/opm_data/grids/wave_orig.grdecl 0
Computing tensor using FIXED boundary conditions.
--------
1) EPS = -1e-16 (negative means absolute residual, not relative!)
156.231510034364 0.00000000000000 0.00000000000000
0.00000000000000 161.928040253848 0.00000000000000
0.00000000000000 0.00000000000000 31.6524093079307
real 0m43.129s
user 0m40.558s
sys 0m2.447s
2) EPS = -1e-8
156.231510034373 0.00000000000000 0.00000000000000
0.00000000000000 161.928040253847 0.00000000000000
0.00000000000000 0.00000000000000 31.6524093079313
real 0m40.681s
user 0m38.129s
sys 0m2.456s
3) EPS = 1e-16 (relative residual)
156.231510034364 0.00000000000000 0.00000000000000
0.00000000000000 161.928040253848 0.00000000000000
0.00000000000000 0.00000000000000 31.6524093079307
real 0m43.110s
user 0m40.551s
sys 0m2.451s
4) EPS = 1e-10
156.231510032888 0.00000000000000 0.00000000000000
0.00000000000000 161.928040251955 0.00000000000000
0.00000000000000 0.00000000000000 31.6524093060648
real 0m37.090s
user 0m34.545s
sys 0m2.426s
5) EPS = 1e-8
156.231510035841 0.00000000000000 0.00000000000000
0.00000000000000 161.928040509462 0.00000000000000
0.00000000000000 0.00000000000000 31.6524097342785
real 0m33.974s
user 0m31.551s
sys 0m2.395s
6) EPS = 1e-5
156.232835862799 0.00000000000000 0.00000000000000
0.00000000000000 161.927447816447 0.00000000000000
0.00000000000000 0.00000000000000 31.6513325481661
real 0m29.270s
user 0m26.752s
sys 0m2.432s
7) EPS = 1e-3
156.289741449408 0.00000000000000 0.00000000000000
0.00000000000000 161.963502056854 0.00000000000000
0.00000000000000 0.00000000000000 31.7244015876773
real 0m26.270s
user 0m23.765s
sys 0m2.426s
8) EPS = 1e-2
155.645697059578 0.00000000000000 0.00000000000000
0.00000000000000 161.717893480317 0.00000000000000
0.00000000000000 0.00000000000000 38.3150616056072
real 0m24.720s
user 0m22.243s
sys 0m2.414s
9) Not running the linear solver (commented out call to SAMG).
real 0m14.999s
user 0m12.867s
sys 0m2.071s
real 0m14.911s
user 0m12.829s
sys 0m2.039s
Old code (delivered december 2008).
==================================
Command:
-------
atgeirr@elrond:~/workdir/trunk/zzz_work_here/upscaling/test_fall_2009$ time ../../../upscaling/upscale_singlephase config.xml
-------
1) tolerance = 1e-10
[ 156.192703205124 0 0 |
| 0 159.455859912012 0 |
| 0 0 30.0168798433113 ]
real 1m57.181s
user 1m6.030s
sys 0m50.704s
New code (the dune-cornerpoint module).
======================================
Command:
-------
atgeirr@elrond:~/dune_release/dune-cornerpoint/upscaling/test$ time ./upscaling_test fileformat=eclipse filename=~/workdir/data/opm_data/grids/wave_orig.grdecl use_unique_boundary_ids=false boundary_condition_type=0 residual_tolerance=1e-8
-------
1) Debug mode (all others are in release mode), residual_tolerance = 1e-8
[done to verify results with optimization, result identical => ok]
156.108070493624 0 0
0 161.812986340526 0
0 0 31.2464464394367
real 4m58.007s
user 4m49.991s
sys 0m7.027s
2) residual_tolerance = 1e-16
156.108072515849 0 0
0 161.812984342818 0
0 0 31.2464597239113
real 0m47.023s
user 0m42.571s
sys 0m4.231s
[again, Dec. 9 2009 - changes in AMG affect solution, I guess I should have run this case in October, too]
156.108072515843 0 0
0 161.812984342796 0
0 0 31.2464597239132
real 0m42.342s
user 0m41.632s
sys 0m0.638s
[again, Apr. 16 2010 - all the changes since last have modified the solution a little]
156.108072515789 0 0
0 161.812984342748 0
0 0 31.2464597239119
real 0m38.425s
user 0m37.649s
sys 0m0.602s
3) residual_tolerance = 1e-10
156.108072556463 0 0
0 161.812984353691 0
0 0 31.2464593728902
real 0m35.906s
user 0m31.487s
sys 0m4.179s
4) residual_tolerance = 1e-8
156.108070493624 0 0
0 161.812986340526 0
0 0 31.2464464394367
real 0m32.599s
user 0m28.111s
sys 0m4.231s
[again, Oct. 8 2009 - changes in AMG affect solution]
156.108070306758 0 0
0 161.812982946198 0
0 0 31.2464373309417
real 0m31.417s
user 0m27.212s
sys 0m4.058s
[again, Oct. 8 2009 - this time with branch version of Incom* and Mimetic*]
156.108070306758 0 0
0 161.812982946198 0
0 0 31.2464373309417
real 0m31.228s
user 0m27.018s
sys 0m4.067s
[again, Dec. 9 2009 - AMG changes again]
156.108068958263 0 0
0 161.812985066158 0
0 0 31.2464395640301
real 0m29.376s
user 0m28.700s
sys 0m0.617s
[again, Mar. 4 2010 - AMG changes again I think, maybe other changes as well. Intersection iterator improvements?]
156.108068958234 0 0
0 161.812985066163 0
0 0 31.24643956403
real 0m26.003s
user 0m24.991s
sys 0m0.597s
[again, Mar. 26 2010 - After changing boundary condition object to deal with Faces instead of boundary ids. Too slow?]
156.108068958234 0 0
0 161.812985066163 0
0 0 31.24643956403
real 0m27.530s
user 0m25.903s
sys 0m0.699s
real 0m25.990s [repeating: faster now, maybe error in experiment before?]
user 0m25.052s
sys 0m0.631s
real 0m25.536s
user 0m24.881s
sys 0m0.602s
[again, Apr. 16 2010 - after refactoring the upscaler to use ...TracerFluid. Apparently no big change.]
156.108068958209 0 0
0 161.81298506611 0
0 0 31.2464395640289
real 0m25.737s
user 0m24.583s
sys 0m0.581s
real 0m25.272s
user 0m24.638s
sys 0m0.584s
-------------------8<------------------------------8<-------------------
[Again, 2010-05-31 -- Five runs. Different computer (Ubuntu Linux 10.04)]
{ for i in $(seq 1 5); do
time ./test/upscaling_test $XML_DATA_FILE;
done; } 2>&1 | perl -ne 'print if /^Upscaled K/ .. /^sys/'
156.108069013645 0 0
0 161.812984922129 0
0 0 31.2464395657857
real 0m13.860s
user 0m13.740s
sys 0m0.120s
---
156.108069013645 0 0
0 161.812984922129 0
0 0 31.2464395657857
real 0m13.902s
user 0m13.800s
sys 0m0.090s
---
156.108069013645 0 0
0 161.812984922129 0
0 0 31.2464395657857
real 0m13.881s
user 0m13.750s
sys 0m0.130s
---
156.108069013645 0 0
0 161.812984922129 0
0 0 31.2464395657857
real 0m13.939s
user 0m13.750s
sys 0m0.190s
---
156.108069013645 0 0
0 161.812984922129 0
0 0 31.2464395657857
real 0m13.912s
user 0m13.800s
sys 0m0.110s
-------------------8<------------------------------8<-------------------
5) residual_tolerance = 1e-5
156.105632072812 0 0
0 161.811906903242 0
0 0 31.2401740405732
real 0m26.004s
user 0m21.783s
sys 0m4.110s
[again, Mar. 4 2010 - see comments above]
156.109771894298 0 0
0 161.809594220408 0
0 0 31.2442748077186
real 0m20.700s
user 0m20.018s
sys 0m0.593s
[again, Apr. 16 2010]
156.109771894273 0 0
0 161.809594220355 0
0 0 31.2442748077174
real 0m20.157s
user 0m19.561s
sys 0m0.555s
6) residual_tolerance = 1e-3
157.903332537278 0 0
0 159.271030166896 0
0 0 32.3461234723618
real 0m21.946s
user 0m17.519s
sys 0m4.194s
[again, Dec. 9 2009 - AMG changes modify solution a lot]
152.545493422819 0 0
0 158.642855651543 0
0 0 27.2903819981729
real 0m20.066s
user 0m19.440s
sys 0m0.603s
7) residual_tolerance = 1e-2
229.190446701104 0 0
0 169.760196673778 0
0 0 56.7072225968723
real 0m21.105s
user 0m16.524s
sys 0m4.207s
8) Not running the linear solver.
real 0m11.291s
user 0m7.524s
sys 0m3.723s
real 0m11.242s
user 0m7.507s
sys 0m3.710s
[again - after a little optimization in flow solver initialization]
real 0m10.426s
user 0m7.031s
sys 0m3.372s
real 0m10.543s
user 0m7.065s
sys 0m3.428s
[again Dec. 9 2009 - maybe Snow Leopard has some of the credit?]
real 0m5.828s
user 0m5.482s
sys 0m0.334s
9) Mar. 26 2010, residual_tolerance = 1e-8, boundary_condition_type = 2 (periodic)
155.314828184991 -0.0871907368904636 0.165154074915333
-0.0871902640179007 158.392114140835 0.00515244278494333
0.165214873343881 0.00515360711329919 29.0905335339699
real 0m36.793s
user 0m35.756s
sys 0m0.655s