From 5e78fc1c9f347b309a014c61361503ee6c289837 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 29 Mar 2012 13:05:59 +0200 Subject: [PATCH 1/5] Multiple changes dealing with initialization. Work in progress. - Moved simulator state class into its own file. - Using new initState...() methods in spu_2p.cpp - No longer controlled by 'scenario' parameter. --- opm/core/TwophaseState.hpp | 85 +++++++++++++++++++++++++++++ opm/core/utility/initState_impl.hpp | 16 ++++++ 2 files changed, 101 insertions(+) create mode 100644 opm/core/TwophaseState.hpp diff --git a/opm/core/TwophaseState.hpp b/opm/core/TwophaseState.hpp new file mode 100644 index 000000000..76253d448 --- /dev/null +++ b/opm/core/TwophaseState.hpp @@ -0,0 +1,85 @@ +/* + 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 . +*/ + +#ifndef OPM_TWOPHASESTATE_HEADER_INCLUDED +#define OPM_TWOPHASESTATE_HEADER_INCLUDED + +#include +#include +#include + +namespace Opm +{ + + 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); + } + + enum ExtremalSat { MinSat, MaxSat }; + + void setWaterSat(const std::vector& cells, + const Opm::IncompPropertiesInterface& props, + ExtremalSat es) + { + const int n = cells.size(); + std::vector smin(2*n); + std::vector 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& pressure () { return press_ ; } + std::vector& facepressure() { return fpress_; } + std::vector& faceflux () { return flux_ ; } + std::vector& saturation () { return sat_ ; } + + const std::vector& pressure () const { return press_ ; } + const std::vector& facepressure() const { return fpress_; } + const std::vector& faceflux () const { return flux_ ; } + const std::vector& saturation () const { return sat_ ; } + + private: + std::vector press_ ; + std::vector fpress_; + std::vector flux_ ; + std::vector sat_ ; + }; + +} // namespace Opm + + +#endif // OPM_TWOPHASESTATE_HEADER_INCLUDED diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index 6815bc391..c24e82548 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -159,6 +159,7 @@ namespace Opm if (state.numPhases() != 2) { THROW("initStateTwophaseFromDeck(): state must have two phases."); } + state.init(grid); const int num_cells = props.numCells(); const bool convection_testcase = param.getDefault("convection_testcase", false); const bool segregation_testcase = param.getDefault("segregation_testcase", false); @@ -179,6 +180,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("water_oil_contact"); initWaterOilContact(grid, props, woc, WaterAbove, state); @@ -187,6 +195,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("water_oil_contact"); initWaterOilContact(grid, props, woc, WaterBelow, state); @@ -240,6 +255,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(); From ed5a7802b3e2283e53ef876e7138bdc49bbf9a2a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 29 Mar 2012 21:10:14 +0200 Subject: [PATCH 2/5] Bugfix: ensure all cells have a valid saturation initially. --- opm/core/utility/initState_impl.hpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index c24e82548..4bbd0250b 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -161,6 +161,12 @@ namespace Opm } state.init(grid); const int num_cells = props.numCells(); + // By default: initialise water saturation to minimum everywhere. + std::vector 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) { @@ -222,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 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]; From baf099cd5a562a9f17885c230cab044a0716fe11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 30 Mar 2012 16:11:07 +0200 Subject: [PATCH 3/5] Minor modification of residual formula (equivalent to the old). --- opm/core/transport/reorder/TransportModelTwophase.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index 806ef7c84..dc9cd5316 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -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); } }; From d397fd847ba765bc8b6dea51d2e6959947f9a20b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 2 Apr 2012 13:24:57 +0200 Subject: [PATCH 4/5] Changed single-cell solver call. - Using [0,1] interval instead of [smin, smax] interval to handle compressible case. - Using new version of Regula Falsi function which exploits initial guess. --- opm/core/transport/reorder/TransportModelTwophase.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index dc9cd5316..33a91905a 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -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); } From ee4123c3d8718469763b08be43b4bac3506c2d29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 2 Apr 2012 15:41:13 +0200 Subject: [PATCH 5/5] Made state before init() valid (0.0 water sat, 1.0 oil sat). --- opm/core/TwophaseState.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/core/TwophaseState.hpp b/opm/core/TwophaseState.hpp index 76253d448..74e4b0adf 100644 --- a/opm/core/TwophaseState.hpp +++ b/opm/core/TwophaseState.hpp @@ -27,6 +27,7 @@ namespace Opm { + /// Simulator state for a two-phase simulator. class TwophaseState { public: @@ -37,6 +38,9 @@ namespace Opm 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 };