This commit is contained in:
Kjetil Olsen Lye 2012-04-10 14:48:35 +02:00
commit 65d3ad547e
3 changed files with 115 additions and 8 deletions

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

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);
}

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();