OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

merge

Browse Source
This commit is contained in:
Kjetil Olsen Lye 2012-04-10 14:48:35 +02:00
commit b39d5c823f
7 changed files with 254 additions and 232 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Makefile.am
View File

@ -147,6 +147,8 @@ opm/core/utility/StopWatch.hpp \
opm/core/utility/UniformTableLinear.hpp \
opm/core/utility/Units.hpp \
opm/core/utility/buildUniformMonotoneTable.hpp \
opm/core/utility/initState.hpp \
opm/core/utility/initState_impl.hpp \
opm/core/utility/linInt.hpp \
opm/core/utility/linearInterpolation.hpp \
opm/core/utility/miscUtilities.hpp \
@ -173,6 +175,7 @@ opm/core/ProductionSpecification.hpp \
opm/core/ColumnExtract.hpp \
opm/core/GridAdapter.hpp \
opm/core/GridManager.hpp \
opm/core/TwophaseState.hpp \
opm/core/WellsManager.hpp \
opm/core/pressure/fsh.h \
opm/core/pressure/HybridPressureSolver.hpp \

249
examples/spu_2p.cpp
View File

@ -45,6 +45,7 @@
#include <opm/core/newwells.h>
#include <opm/core/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/utility/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/Units.hpp>
@ -73,6 +74,7 @@
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/ColumnExtract.hpp>
#include <opm/core/TwophaseState.hpp>
#include <opm/core/transport/GravityColumnSolver.hpp>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
@ -96,104 +98,8 @@
class ReservoirState
{
public:
ReservoirState(const UnstructuredGrid* g, const double init_sat, const double init_p)
: press_ (g->number_of_cells, 0.0),
fpress_(g->number_of_faces, 0.0),
flux_ (g->number_of_faces, 0.0),
sat_ (2 * g->number_of_cells, 0.0)
{
for (int cell = 0; cell < g->number_of_cells; ++cell) {
sat_[2*cell] = init_sat;
sat_[2*cell + 1] = 1.0 - init_sat;
press_[cell] = init_p;
}
}
enum ExtremalSat { MinSat, MaxSat };
void setToMinimumWaterSat(const Opm::IncompPropertiesInterface& props)
{
const int n = props.numCells();
std::vector<int> cells(n);
for (int i = 0; i < n; ++i) {
cells[i] = i;
}
setWaterSat(cells, props, MinSat);
}
void setWaterSat(const std::vector<int>& cells,
const Opm::IncompPropertiesInterface& props,
ExtremalSat es)
{
const int n = cells.size();
std::vector<double> smin(2*n);
std::vector<double> smax(2*n);
props.satRange(n, &cells[0], &smin[0], &smax[0]);
const double* svals = (es == MinSat) ? &smin[0] : &smax[0];
for (int ci = 0; ci < n; ++ci) {
const int cell = cells[ci];
sat_[2*cell] = svals[2*ci];
sat_[2*cell + 1] = 1.0 - sat_[2*cell];
}
}
// Initialize saturations so that there is water below woc,
// and oil above.
// TODO: add 'anitialiasing', obtaining a more precise woc
// by f. ex. subdividing cells cut by the woc.
void initWaterOilContact(const UnstructuredGrid& grid,
const Opm::IncompPropertiesInterface& props,
const double woc)
{
// Find out which cells should have water and which should have oil.
std::vector<int> oil;
std::vector<int> water;
const int num_cells = grid.number_of_cells;
oil.reserve(num_cells);
water.reserve(num_cells);
const int dim = grid.dimensions;
for (int c = 0; c < num_cells; ++c) {
const double z = grid.cell_centroids[dim*c + dim - 1];
if (z > woc) {
// Z is depth, we put water in the deepest parts
// (even if oil is heavier...).
water.push_back(c);
} else {
oil.push_back(c);
}
}
// Set saturations.
setWaterSat(oil, props, MinSat);
setWaterSat(water, props, MaxSat);
}
int numPhases() const { return sat_.size()/press_.size(); }
std::vector<double>& pressure () { return press_ ; }
std::vector<double>& facepressure() { return fpress_; }
std::vector<double>& faceflux () { return flux_ ; }
std::vector<double>& saturation () { return sat_ ; }
const std::vector<double>& pressure () const { return press_ ; }
const std::vector<double>& facepressure() const { return fpress_; }
const std::vector<double>& faceflux () const { return flux_ ; }
const std::vector<double>& saturation () const { return sat_ ; }
private:
std::vector<double> press_ ;
std::vector<double> fpress_;
std::vector<double> flux_ ;
std::vector<double> sat_ ;
};
static void outputState(const UnstructuredGrid& grid,
const ReservoirState& state,
const Opm::TwophaseState& state,
const int step,
const std::string& output_dir)
{
@ -374,8 +280,8 @@ main(int argc, char** argv)
boost::scoped_ptr<Opm::WellsManager> wells;
boost::scoped_ptr<Opm::RockCompressibility> rock_comp;
Opm::SimulatorTimer simtimer;
double water_oil_contact = 0.0;
bool woc_set = false;
Opm::TwophaseState state;
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename);
@ -393,16 +299,12 @@ main(int argc, char** argv)
} else {
simtimer.init(param);
}
// Water-oil contact.
if (deck.hasField("EQUIL")) {
water_oil_contact = deck.getEQUIL().equil[0].water_oil_contact_depth_;
woc_set = true;
} else if (param.has("water_oil_contact")) {
water_oil_contact = param.get<double>("water_oil_contact");
woc_set = true;
}
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(deck));
// Gravity.
gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
// Init state variables (saturation and pressure).
initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
@ -418,23 +320,20 @@ main(int argc, char** argv)
wells.reset(new Opm::WellsManager());
// Timer init.
simtimer.init(param);
if (param.has("water_oil_contact")) {
water_oil_contact = param.get<double>("water_oil_contact");
woc_set = true;
}
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state);
}
// Extra fluid init for transport solver.
TwophaseFluid fluid(*props);
// Gravity init.
double gravity[3] = { 0.0 };
double g = param.getDefault("gravity", 0.0);
bool use_gravity = g != 0.0;
// Warn if gravity but no density difference.
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
if (use_gravity) {
gravity[grid->c_grid()->dimensions - 1] = g;
if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
@ -468,17 +367,11 @@ main(int argc, char** argv)
nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
}
// State-related and source-related variables init.
// Source-related variables init.
int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> totmob;
std::vector<double> omega; // Will remain empty if no gravity.
std::vector<double> rc; // Will remain empty if no rock compressibility.
double init_sat = param.getDefault("init_sat", 0.0);
double init_p = param.getDefault("init_p_bar", 235)*Opm::unit::barsa;
ReservoirState state(grid->c_grid(), init_sat, init_p);
if (!param.has("init_sat")) {
state.setToMinimumWaterSat(*props);
}
// Extra rock init.
std::vector<double> porevol;
@ -493,106 +386,18 @@ main(int argc, char** argv)
// code expects a scalar sw, not both sw and so.
std::vector<double> reorder_sat(num_cells);
std::vector<double> src(num_cells, 0.0);
int scenario = param.getDefault("scenario", woc_set ? 4 : 0);
switch (scenario) {
case 0:
{
std::cout << "==== Scenario 0: simple wells or single-cell source and sink.\n";
if (wells->c_wells()) {
Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
} else {
double flow_per_sec = 0.1*tot_porevol_init/Opm::unit::day;
if (param.has("injection_rate_per_day")) {
flow_per_sec = param.get<double>("injection_rate_per_day")/Opm::unit::day;
}
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
}
break;
}
case 1:
{
std::cout << "==== Scenario 1: half source, half sink.\n";
double flow_per_sec = 0.1*porevol[0]/Opm::unit::day;
std::fill(src.begin(), src.begin() + src.size()/2, flow_per_sec);
std::fill(src.begin() + src.size()/2, src.end(), -flow_per_sec);
break;
}
case 2:
{
std::cout << "==== Scenario 2: gravity convection.\n";
if (!use_gravity) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (use_deck) {
std::cout << "**** Warning: running gravity convection scenario, which expects a cartesian grid."
<< std::endl;
}
if (grid->c_grid()->cartdims[2] <= 1) {
std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
}
std::vector<int> left_cells;
left_cells.reserve(num_cells/2);
const int *glob_cell = grid->c_grid()->global_cell;
for (int cell = 0; cell < num_cells; ++cell) {
const int* cd = grid->c_grid()->cartdims;
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool left = (gc % cd[0]) < cd[0]/2;
if (left) {
left_cells.push_back(cell);
}
}
state.setWaterSat(left_cells, *props, ReservoirState::MaxSat);
break;
}
case 3:
{
std::cout << "==== Scenario 3: gravity segregation.\n";
if (!use_gravity) {
std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
}
if (use_deck) {
std::cout << "**** Warning: running gravity segregation scenario, which expects a cartesian grid."
<< std::endl;
}
if (grid->c_grid()->cartdims[2] <= 1) {
std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl;
}
std::vector<int> top_cells;
const int *glob_cell = grid->c_grid()->global_cell;
// Water on top
for (int cell = 0; cell < num_cells; ++cell) {
const int* cd = grid->c_grid()->cartdims;
const int gc = glob_cell == 0 ? cell : glob_cell[cell];
bool top = (gc / cd[0] / cd[1]) < cd[2]/2;
if (top) {
top_cells.push_back(cell);
}
}
state.setWaterSat(top_cells, *props, ReservoirState::MaxSat);
break;
}
case 4:
{
std::cout << "==== Scenario 4: water-oil contact and simple wells or sources\n";
if (!use_gravity) {
std::cout << "**** Warning: initializing segregated water and oil zones, but gravity is zero." << std::endl;
}
state.initWaterOilContact(*grid->c_grid(), *props, water_oil_contact);
if (wells->c_wells()) {
Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
} else {
double flow_per_sec = 0.01*tot_porevol_init/Opm::unit::day;
src[0] = flow_per_sec;
src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
}
break;
}
default:
{
THROW("==== Scenario " << scenario << " is unknown.");
}
// Initialising src
if (wells->c_wells()) {
Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
} else {
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
*tot_porevol_init/Opm::unit::day;
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
}
TransportSource* tsrc = create_transport_source(2, 2);
double ssrc[] = { 1.0, 0.0 };
double ssink[] = { 0.0, 1.0 };

89
opm/core/TwophaseState.hpp Normal file
View File

@ -0,0 +1,89 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TWOPHASESTATE_HEADER_INCLUDED
#define OPM_TWOPHASESTATE_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <vector>
namespace Opm
{
/// Simulator state for a two-phase simulator.
class TwophaseState
{
public:
void init(const UnstructuredGrid& g)
{
press_.resize(g.number_of_cells, 0.0);
fpress_.resize(g.number_of_faces, 0.0);
flux_.resize(g.number_of_faces, 0.0);
sat_.resize(2 * g.number_of_cells, 0.0);
for (int cell = 0; cell < g.number_of_cells; ++cell) {
sat_[2*cell + 1] = 1.0; // Defaulting oil saturations to 1.0.
}
}
enum ExtremalSat { MinSat, MaxSat };
void setWaterSat(const std::vector<int>& cells,
const Opm::IncompPropertiesInterface& props,
ExtremalSat es)
{
const int n = cells.size();
std::vector<double> smin(2*n);
std::vector<double> smax(2*n);
props.satRange(n, &cells[0], &smin[0], &smax[0]);
const double* svals = (es == MinSat) ? &smin[0] : &smax[0];
for (int ci = 0; ci < n; ++ci) {
const int cell = cells[ci];
sat_[2*cell] = svals[2*ci];
sat_[2*cell + 1] = 1.0 - sat_[2*cell];
}
}
int numPhases() const
{
return 2;
}
std::vector<double>& pressure () { return press_ ; }
std::vector<double>& facepressure() { return fpress_; }
std::vector<double>& faceflux () { return flux_ ; }
std::vector<double>& saturation () { return sat_ ; }
const std::vector<double>& pressure () const { return press_ ; }
const std::vector<double>& facepressure() const { return fpress_; }
const std::vector<double>& faceflux () const { return flux_ ; }
const std::vector<double>& saturation () const { return sat_ ; }
private:
std::vector<double> press_ ;
std::vector<double> fpress_;
std::vector<double> flux_ ;
std::vector<double> sat_ ;
};
} // namespace Opm
#endif // OPM_TWOPHASESTATE_HEADER_INCLUDED

5
opm/core/transport/reorder/TransportModelTwophase.cpp
View File

@ -154,7 +154,7 @@ namespace Opm
}
double operator()(double s) const
{
return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx) + dtpv*s*comp_term;
return s - s0 + dtpv*(outflux*tm.fracFlow(s, cell) + influx + s*comp_term);
}
};
@ -167,7 +167,8 @@ namespace Opm
// return;
// }
int iters_used;
saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
// saturation_[cell] = modifiedRegulaFalsi(res, smin_[2*cell], smax_[2*cell], maxit_, tol_, iters_used);
saturation_[cell] = modifiedRegulaFalsi(res, saturation_[cell], 0.0, 1.0, maxit_, tol_, iters_used);
fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
}

102
opm/core/utility/RootFinders.hpp
View File

@ -137,6 +137,108 @@ namespace Opm
}
/// Implements a modified regula falsi method as described in
/// "Improved algorithms of Illinois-type for the numerical
/// solution of nonlinear equations"
/// by J. A. Ford.
/// Current variant is the 'Pegasus' method.
/// This version takes an extra parameter for the initial guess.
template <class Functor>
inline double modifiedRegulaFalsi(const Functor& f,
const double initial_guess,
const double a,
const double b,
const int max_iter,
const double tolerance,
int& iterations_used)
{
using namespace std;
const double macheps = numeric_limits<double>::epsilon();
const double eps = tolerance + macheps*max(max(fabs(a), fabs(b)), 1.0);
double f_initial = f(initial_guess);
const double epsF = tolerance + macheps*max(fabs(f_initial), 1.0);
if (fabs(f_initial) < epsF) {
return initial_guess;
}
double x0 = a;
double x1 = b;
double f0 = f_initial;
double f1 = f_initial;
if (x0 != initial_guess) {
f0 = f(x0);
if (fabs(f0) < epsF) {
return x0;
}
}
if (x1 != initial_guess) {
f1 = f(x1);
if (fabs(f1) < epsF) {
return x1;
}
}
if (f0*f_initial < 0.0) {
x1 = initial_guess;
f1 = f_initial;
} else {
x0 = initial_guess;
f0 = f_initial;
}
if (f0*f1 > 0.0) {
THROW("Error in parameters, zero not bracketed: [a, b] = ["
<< a << ", " << b << "] fa = " << f0 << " fb = " << f1);
}
iterations_used = 0;
// In every iteraton, x1 is the last point computed,
// and x0 is the last point computed that makes it a bracket.
while (fabs(x1 - x0) >= 1e-9*eps) {
double xnew = regulaFalsiStep(x0, x1, f0, f1);
double fnew = f(xnew);
// cout << "xnew = " << xnew << " fnew = " << fnew << endl;
++iterations_used;
if (iterations_used > max_iter) {
THROW("Maximum number of iterations exceeded.\n"
<< "Current interval is [" << min(x0, x1) << ", "
<< max(x0, x1) << "]");
}
if (fabs(fnew) < epsF) {
return xnew;
}
// Now we must check which point we must replace.
if ((fnew > 0.0) == (f0 > 0.0)) {
// We must replace x0.
x0 = x1;
f0 = f1;
} else {
// We must replace x1, this is the case where
// the modification to regula falsi kicks in,
// by modifying f0.
// 1. The classic Illinois method
// const double gamma = 0.5;
// @afr: The next two methods do not work??!!?
// 2. The method called 'Method 3' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi1/(1.0 - phi0);
// 3. The method called 'Method 4' in the paper.
// const double phi0 = f1/f0;
// const double phi1 = fnew/f1;
// const double gamma = 1.0 - phi0 - phi1;
// cout << "phi0 = " << phi0 <<" phi1 = " << phi1 <<
// " gamma = " << gamma << endl;
// 4. The 'Pegasus' method
const double gamma = f1/(f1 + fnew);
f0 *= gamma;
}
x1 = xnew;
f1 = fnew;
}
return 0.5*(x0 + x1);
}
template <class Functor>
inline void bracketZero(const Functor& f,
const double x0,

29
opm/core/utility/initState_impl.hpp
View File

@ -159,7 +159,14 @@ namespace Opm
if (state.numPhases() != 2) {
THROW("initStateTwophaseFromDeck(): state must have two phases.");
}
state.init(grid);
const int num_cells = props.numCells();
// By default: initialise water saturation to minimum everywhere.
std::vector<int> all_cells(num_cells);
for (int i = 0; i < num_cells; ++i) {
all_cells[i] = i;
}
state.setWaterSat(all_cells, props, State::MinSat);
const bool convection_testcase = param.getDefault("convection_testcase", false);
const bool segregation_testcase = param.getDefault("segregation_testcase", false);
if (convection_testcase) {
@ -179,6 +186,13 @@ namespace Opm
const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa;
std::fill(state.pressure().begin(), state.pressure().end(), init_p);
} else if (segregation_testcase) {
// Warn against error-prone usage.
if (gravity == 0.0) {
std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
}
if (grid.cartdims[2] <= 1) {
std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl;
}
// Initialise water saturation to max *above* water-oil contact.
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterAbove, state);
@ -187,6 +201,13 @@ namespace Opm
double dens[2] = { props.density()[1], props.density()[0] };
initHydrostaticPressure(grid, dens, woc, gravity, woc, ref_p, state);
} else if (param.has("water_oil_contact")) {
// Warn against error-prone usage.
if (gravity == 0.0) {
std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
}
if (grid.cartdims[2] <= 1) {
std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
}
// Initialise water saturation to max below water-oil contact.
const double woc = param.get<double>("water_oil_contact");
initWaterOilContact(grid, props, woc, WaterBelow, state);
@ -207,12 +228,7 @@ namespace Opm
const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1];
initHydrostaticPressure(grid, dens, ref_z, gravity, ref_z, ref_p, state);
} else {
// By default: initialise water saturation to minimum everywhere.
std::vector<int> all_cells(num_cells);
for (int i = 0; i < num_cells; ++i) {
all_cells[i] = i;
}
state.setWaterSat(all_cells, props, State::MinSat);
// Use default: water saturation is minimum everywhere.
// Initialise pressure to hydrostatic state.
const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa;
const double rho = props.density()[1];
@ -240,6 +256,7 @@ namespace Opm
if (state.numPhases() != 2) {
THROW("initStateTwophaseFromDeck(): state must have two phases.");
}
state.init(grid);
if (deck.hasField("EQUIL")) {
// Set saturations depending on oil-water contact.
const EQUIL& equil= deck.getEQUIL();

9
opm/core/utility/writeVtkData.cpp
View File

@ -283,7 +283,6 @@ namespace Opm
}
Tag celldatatag("CellData", pm, os);
pm.clear();
pm["type"] = "Int32";
pm["NumberOfComponents"] = "1";
pm["format"] = "ascii";
pm["type"] = "Float64";
@ -298,7 +297,13 @@ namespace Opm
if (item % num_per_line == 0) {
Tag::indent(os);
}
os << field[item] << ' ';
double value = field[item];
if (std::fabs(value) < std::numeric_limits<double>::min()) {
// Avoiding denormal numbers to work around
// bug in Paraview.
value = 0.0;
}
os << value << ' ';
if (item % num_per_line == num_per_line - 1
|| item == num_cells - 1) {
os << '\n';