From 72f8da074e61a901c9633bc8ad05c5a3e35cee82 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 25 Oct 2012 13:23:50 +0200 Subject: [PATCH 1/8] Add constructor taking file name to GridManager. Calls read_grid(). Grid fileformat not yet documented. --- opm/core/GridManager.cpp | 14 ++++++++++++++ opm/core/GridManager.hpp | 6 ++++++ 2 files changed, 20 insertions(+) diff --git a/opm/core/GridManager.cpp b/opm/core/GridManager.cpp index d286e504..528e2e12 100644 --- a/opm/core/GridManager.cpp +++ b/opm/core/GridManager.cpp @@ -97,6 +97,20 @@ namespace Opm + /// Construct a grid from an input file. + /// The file format used is currently undocumented, + /// and is therefore only suited for internal use. + GridManager::GridManager(const std::string& input_filename) + { + ug_ = read_grid(input_filename.c_str()); + if (!ug_) { + THROW("Failed to read grid from file " << input_filename); + } + } + + + + /// Destructor. GridManager::~GridManager() { diff --git a/opm/core/GridManager.hpp b/opm/core/GridManager.hpp index 26e335f2..36846f2c 100644 --- a/opm/core/GridManager.hpp +++ b/opm/core/GridManager.hpp @@ -20,6 +20,7 @@ #ifndef OPM_GRIDMANAGER_HEADER_INCLUDED #define OPM_GRIDMANAGER_HEADER_INCLUDED +#include struct UnstructuredGrid; @@ -53,6 +54,11 @@ namespace Opm GridManager(int nx, int ny, int nz, double dx, double dy, double dz); + /// Construct a grid from an input file. + /// The file format used is currently undocumented, + /// and is therefore only suited for internal use. + GridManager(const std::string& input_filename); + /// Destructor. ~GridManager(); From 9152d78405d6f06795a2d0a4b43f0560fc4bd2f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 25 Oct 2012 14:47:26 +0200 Subject: [PATCH 2/8] Added new tof computation utility. This program assumes that grid, porosity, flux field and sources are available to read from files. It does not compute any flux field itself. --- examples/Makefile.am | 4 + examples/compute_tof_from_files.cpp | 173 ++++++++++++++++++++++++++++ 2 files changed, 177 insertions(+) create mode 100644 examples/compute_tof_from_files.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index 043a44c3..85dce93a 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -31,6 +31,7 @@ $(BOOST_SYSTEM_LIB) noinst_PROGRAMS = \ compute_tof \ +compute_tof_from_files \ refine_wells \ scaneclipsedeck \ sim_2p_comp_reorder \ @@ -54,6 +55,9 @@ endif compute_tof_SOURCES = compute_tof.cpp compute_tof_LDADD = $(LDADD) $(LINK_BOOST_FILESYSTEM) +compute_tof_from_files_SOURCES = compute_tof_from_files.cpp +compute_tof_from_files_LDADD = $(LDADD) $(LINK_BOOST_FILESYSTEM) + refine_wells_SOURCES = refine_wells.cpp if HAVE_ERT diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp new file mode 100644 index 00000000..30b76f7b --- /dev/null +++ b/examples/compute_tof_from_files.cpp @@ -0,0 +1,173 @@ +/* + 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 . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + + +namespace +{ + void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) + { + if (param.anyUnused()) { + std::cout << "-------------------- Warning: unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "-------------------------------------------------------------------------" << std::endl; + } + } +} // anon namespace + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + parameter::ParameterGroup param(argc, argv, false); + + // Read grid. + GridManager grid_manager(param.get("grid_filename")); + const UnstructuredGrid& grid = *grid_manager.c_grid(); + + // Read porosity, compute pore volume. + std::vector porevol; + { + std::ifstream poro_stream(param.get("poro_filename").c_str()); + std::istream_iterator beg(poro_stream); + std::istream_iterator end; + porevol.assign(beg, end); // Now contains poro. + if (int(porevol.size()) != grid.number_of_cells) { + THROW("Size of porosity field differs from number of cells."); + } + for (int i = 0; i < grid.number_of_cells; ++i) { + porevol[i] *= grid.cell_volumes[i]; + } + } + + // Read flux. + std::vector flux; + { + std::ifstream flux_stream(param.get("flux_filename").c_str()); + std::istream_iterator beg(flux_stream); + std::istream_iterator end; + flux.assign(beg, end); + if (int(flux.size()) != grid.number_of_faces) { + THROW("Size of flux field differs from number of faces."); + } + } + + // Read source terms. + std::vector src; + { + std::ifstream src_stream(param.get("src_filename").c_str()); + std::istream_iterator beg(src_stream); + std::istream_iterator end; + src.assign(beg, end); + if (int(src.size()) != grid.number_of_cells) { + THROW("Size of source term field differs from number of cells."); + } + } + + // Choice of tof solver. + bool use_dg = param.getDefault("use_dg", false); + int dg_degree = -1; + bool use_cvi = false; + if (use_dg) { + dg_degree = param.getDefault("dg_degree", 0); + use_cvi = param.getDefault("use_cvi", false); + } + + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::ofstream epoch_os; + std::string output_dir; + if (output) { + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + param.writeParam(output_dir + "/simulation.param"); + } + + // Issue a warning if any parameters were unused. + warnIfUnusedParams(param); + + // Solve time-of-flight. + Opm::time::StopWatch transport_timer; + transport_timer.start(); + std::vector tof; + if (use_dg) { + Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi); + tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof); + } else { + Opm::TransportModelTracerTof tofsolver(grid); + tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof); + } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + + // Output. + if (output) { + std::string tof_filename = output_dir + "/tof.txt"; + std::ofstream tof_stream(tof_filename.c_str()); + std::copy(tof.begin(), tof.end(), std::ostream_iterator(tof_stream, "\n")); + } +} From 20c1567eca40ceb41da65e16c3bbb96f3cda3cbb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 29 Oct 2012 08:27:59 +0100 Subject: [PATCH 3/8] Increased number of digits in output. --- examples/compute_tof_from_files.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index 30b76f7b..012def4e 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -168,6 +168,7 @@ main(int argc, char** argv) if (output) { std::string tof_filename = output_dir + "/tof.txt"; std::ofstream tof_stream(tof_filename.c_str()); + tof_stream.precision(16); std::copy(tof.begin(), tof.end(), std::ostream_iterator(tof_stream, "\n")); } } From 020f452a28dd731c0e38ccdf132a462fb76acc73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 29 Oct 2012 17:23:17 +0100 Subject: [PATCH 4/8] Added parameters for controlling use of multidim upwinding. For now, you will simply get SPU even with use_multidim_upwind=true. --- examples/compute_tof.cpp | 5 ++++- .../reorder/TransportModelTracerTof.cpp | 17 +++++++++++++---- .../reorder/TransportModelTracerTof.hpp | 7 ++++++- 3 files changed, 23 insertions(+), 6 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index ce04aa59..afe89321 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -169,8 +169,11 @@ main(int argc, char** argv) // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); int dg_degree = -1; + bool use_multidim_upwind = false; if (use_dg) { dg_degree = param.getDefault("dg_degree", 0); + } else { + use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } // Write parameters used for later reference. @@ -231,7 +234,7 @@ main(int argc, char** argv) tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof); transport_timer.stop(); } else { - Opm::TransportModelTracerTof tofsolver(*grid->c_grid()); + Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind); transport_timer.start(); tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); transport_timer.stop(); diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 8d93f248..e5202a6d 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -30,8 +30,9 @@ namespace Opm /// Construct solver. /// \param[in] grid A 2d or 3d grid. - TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid) - : grid_(grid) + TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid, + const bool use_multidim_upwind) + : grid_(grid), use_multidim_upwind_(use_multidim_upwind) { } @@ -94,7 +95,11 @@ namespace Opm // Add flux to upwind_term or downwind_flux, if interior. if (other != -1) { if (flux < 0.0) { - upwind_term += flux*tof_[other]; + if (use_multidim_upwind_) { + upwind_term += flux*multidimUpwindTof(f, other); + } else { + upwind_term += flux*tof_[other]; + } } else { downwind_flux += flux; } @@ -117,6 +122,10 @@ namespace Opm } - + double TransportModelTracerTof::multidimUpwindTof(const int face, const int upwind_cell) const + { + // Just SPU for now... + return tof_[upwind_cell]; + } } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 270ac328..39b280d4 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -43,7 +43,9 @@ namespace Opm public: /// Construct solver. /// \param[in] grid A 2d or 3d grid. - TransportModelTracerTof(const UnstructuredGrid& grid); + /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. + TransportModelTracerTof(const UnstructuredGrid& grid, + const bool use_multidim_upwind = false); /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. @@ -61,12 +63,15 @@ namespace Opm virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); + double multidimUpwindTof(const int face, const int upwind_cell) const; + private: const UnstructuredGrid& grid_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell double* tof_; + bool use_multidim_upwind_; }; } // namespace Opm From ca351e716dce5e6cb63941091d4c74a625ee14aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 30 Oct 2012 10:31:30 +0100 Subject: [PATCH 5/8] Change transport timer use to include construction time. Also minor mods to be more similar to compute_tof_from_files.cpp --- examples/compute_tof.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index afe89321..259bd253 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -169,9 +169,11 @@ main(int argc, char** argv) // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); int dg_degree = -1; + bool use_cvi = false; bool use_multidim_upwind = false; if (use_dg) { dg_degree = param.getDefault("dg_degree", 0); + use_cvi = param.getDefault("use_cvi", false); } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } @@ -226,19 +228,16 @@ main(int argc, char** argv) wells->c_wells(), well_state.perfRates(), transport_src); // Solve time-of-flight. + transport_timer.start(); std::vector tof; if (use_dg) { - bool use_cvi = param.getDefault("use_cvi", false); Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid(), use_cvi); - transport_timer.start(); tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof); - transport_timer.stop(); } else { Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind); - transport_timer.start(); tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); - transport_timer.stop(); } + transport_timer.stop(); double tt = transport_timer.secsSinceStart(); std::cout << "Transport solver took: " << tt << " seconds." << std::endl; ttime += tt; From 386b4e4942b6655b73ec9ef2df359664f82d694f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 30 Oct 2012 10:32:29 +0100 Subject: [PATCH 6/8] Add use_multidim_upwind parameter. --- examples/compute_tof_from_files.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index 012def4e..17ad7b85 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -124,9 +124,12 @@ main(int argc, char** argv) bool use_dg = param.getDefault("use_dg", false); int dg_degree = -1; bool use_cvi = false; + bool use_multidim_upwind = false; if (use_dg) { dg_degree = param.getDefault("dg_degree", 0); use_cvi = param.getDefault("use_cvi", false); + } else { + use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } // Write parameters used for later reference. @@ -157,7 +160,7 @@ main(int argc, char** argv) Opm::TransportModelTracerTofDiscGal tofsolver(grid, use_cvi); tofsolver.solveTof(&flux[0], &porevol[0], &src[0], dg_degree, tof); } else { - Opm::TransportModelTracerTof tofsolver(grid); + Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind); tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof); } transport_timer.stop(); From 3325d4ef39faa17d9ed560b76fc8e2e11a6a02e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 30 Oct 2012 13:10:50 +0100 Subject: [PATCH 7/8] Implemented multidimensional upwinding. The 'SMU' variant is chosen for its smoothness. --- .../reorder/TransportModelTracerTof.cpp | 67 ++++++++++++++++++- .../reorder/TransportModelTracerTof.hpp | 2 + 2 files changed, 66 insertions(+), 3 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index e5202a6d..f4c2e72b 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -64,6 +64,10 @@ namespace Opm tof.resize(grid_.number_of_cells); std::fill(tof.begin(), tof.end(), 0.0); tof_ = &tof[0]; + if (use_multidim_upwind_) { + face_tof_.resize(grid_.number_of_faces); + std::fill(face_tof_.begin(), face_tof_.end(), 0.0); + } reorderAndTransport(grid_, darcyflux); } @@ -96,7 +100,9 @@ namespace Opm if (other != -1) { if (flux < 0.0) { if (use_multidim_upwind_) { - upwind_term += flux*multidimUpwindTof(f, other); + const double ftof = multidimUpwindTof(f, other); + face_tof_[f] = ftof; + upwind_term += flux*ftof; } else { upwind_term += flux*tof_[other]; } @@ -124,8 +130,63 @@ namespace Opm double TransportModelTracerTof::multidimUpwindTof(const int face, const int upwind_cell) const { - // Just SPU for now... - return tof_[upwind_cell]; + // Implements multidim upwind according to + // "Multidimensional upstream weighting for multiphase transport on general grids" + // by Keilegavlen, Kozdon, Mallison. + // However, that article does not give a 3d extension other than noting that using + // multidimensional upwinding in the XY-plane and not in the Z-direction may be + // a good idea. We have here attempted some generalization, by looking at all + // face-neighbours across edges as upwind candidates, and giving them all uniform weight. + // This will over-weight the immediate upstream cell value in an extruded 2d grid with + // one layer (top and bottom no-flow faces will enter the computation) compared to the + // original 2d case. Improvements are welcome. + + // Identify the adjacent faces of the upwind cell. + const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[face]; + const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[face + 1]; + ASSERT(face_nodes_end - face_nodes_beg == 2 || grid_.dimensions != 2); + adj_faces_.clear(); + for (int hf = grid_.cell_facepos[upwind_cell]; hf < grid_.cell_facepos[upwind_cell + 1]; ++hf) { + const int f = grid_.cell_faces[hf]; + if (f != face) { + const int* f_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f]; + const int* f_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1]; + // Find out how many vertices they have in common. + // Using simple linear searches since sets are small. + int num_common = 0; + for (const int* f_iter = f_nodes_beg; f_iter < f_nodes_end; ++f_iter) { + num_common += std::count(face_nodes_beg, face_nodes_end, *f_iter); + } + if (num_common == grid_.dimensions - 1) { + // Neighbours over an edge (3d) or vertex (2d). + adj_faces_.push_back(f); + } else { + ASSERT(num_common == 0); + } + } + } + + // Indentify adjacent faces with inflows, compute omega_star, omega, + // add up contributions. + const int num_adj = adj_faces_.size(); + ASSERT(num_adj == face_nodes_end - face_nodes_beg); + const double flux_face = std::fabs(darcyflux_[face]); + double tof = 0.0; + for (int ii = 0; ii < num_adj; ++ii) { + const int f = adj_faces_[ii]; + const double influx_f = (grid_.face_cells[2*f] == upwind_cell) ? -darcyflux_[f] : darcyflux_[f]; + const double omega_star = influx_f/flux_face; + // SPU + // const double omega = 0.0; + // TMU + // const double omega = omega_star > 0.0 ? std::min(omega_star, 1.0) : 0.0; + // SMU + const double omega = omega_star > 0.0 ? omega_star/(1.0 + omega_star) : 0.0; + tof += (1.0 - omega)*tof_[upwind_cell] + omega*face_tof_[f]; + } + + // For now taking a simple average. + return tof/double(num_adj); } } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 39b280d4..7a27e0ca 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -72,6 +72,8 @@ namespace Opm const double* source_; // one volumetric source term per cell double* tof_; bool use_multidim_upwind_; + std::vector face_tof_; // For multidim upwind face tofs. + mutable std::vector adj_faces_; // For multidim upwind logic. }; } // namespace Opm From 51bbc85a265f74746877d8f6fc7e6d8d02d7f050 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 5 Nov 2012 14:26:00 +0100 Subject: [PATCH 8/8] Bugfix tof computations with multidimensional upwinding. Cell tof depends on downwind face tof in a more complicated way with multidim upwinding, this was not done correctly. --- .../reorder/TransportModelTracerTof.cpp | 92 ++++++++++++++++--- .../reorder/TransportModelTracerTof.hpp | 4 +- 2 files changed, 82 insertions(+), 14 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index f4c2e72b..91b20ef0 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -76,6 +76,10 @@ namespace Opm void TransportModelTracerTof::solveSingleCell(const int cell) { + if (use_multidim_upwind_) { + solveSingleCellMultidimUpwind(cell); + return; + } // Compute flux terms. // Sources have zero tof, and therefore do not contribute // to upwind_term. Sinks on the other hand, must be added @@ -99,13 +103,7 @@ namespace Opm // Add flux to upwind_term or downwind_flux, if interior. if (other != -1) { if (flux < 0.0) { - if (use_multidim_upwind_) { - const double ftof = multidimUpwindTof(f, other); - face_tof_[f] = ftof; - upwind_term += flux*ftof; - } else { - upwind_term += flux*tof_[other]; - } + upwind_term += flux*tof_[other]; } else { downwind_flux += flux; } @@ -119,6 +117,60 @@ namespace Opm + void TransportModelTracerTof::solveSingleCellMultidimUpwind(const int cell) + { + // Compute flux terms. + // Sources have zero tof, and therefore do not contribute + // to upwind_term. Sinks on the other hand, must be added + // to the downwind terms (note sign change resulting from + // different sign conventions: pos. source is injection, + // pos. flux is outflow). + double upwind_term = 0.0; + double downwind_term_cell_factor = std::max(-source_[cell], 0.0); + double downwind_term_face = 0.0; + for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { + int f = grid_.cell_faces[i]; + double flux; + int other; + // Compute cell flux + if (cell == grid_.face_cells[2*f]) { + flux = darcyflux_[f]; + other = grid_.face_cells[2*f+1]; + } else { + flux =-darcyflux_[f]; + other = grid_.face_cells[2*f]; + } + // Add flux to upwind_term or downwind_flux, if interior. + if (other != -1) { + if (flux < 0.0) { + upwind_term += flux*face_tof_[f]; + } else { + double fterm, cterm_factor; + multidimUpwindTerms(f, cell, fterm, cterm_factor); + downwind_term_face += fterm*flux; + downwind_term_cell_factor += cterm_factor*flux; + } + } + } + + // Compute tof for cell. + tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor; // } + + // Compute tof for downwind faces. + for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { + int f = grid_.cell_faces[i]; + const double outflux_f = (grid_.face_cells[2*f] == cell) ? darcyflux_[f] : -darcyflux_[f]; + if (outflux_f > 0.0) { + double fterm, cterm_factor; + multidimUpwindTerms(f, cell, fterm, cterm_factor); + face_tof_[f] = fterm + cterm_factor*tof_[cell]; + } + } + } + + + + void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells) { std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; @@ -128,7 +180,18 @@ namespace Opm } - double TransportModelTracerTof::multidimUpwindTof(const int face, const int upwind_cell) const + + + // Assumes that face_tof_[f] is known for all upstream faces f of upwind_cell. + // Assumes that darcyflux_[face] is != 0.0. + // This function returns factors to compute the tof for 'face': + // tof(face) = face_term + cell_term_factor*tof(upwind_cell). + // It is not computed here, since these factors are needed to + // compute the tof(upwind_cell) itself. + void TransportModelTracerTof::multidimUpwindTerms(const int face, + const int upwind_cell, + double& face_term, + double& cell_term_factor) const { // Implements multidim upwind according to // "Multidimensional upstream weighting for multiphase transport on general grids" @@ -171,7 +234,8 @@ namespace Opm const int num_adj = adj_faces_.size(); ASSERT(num_adj == face_nodes_end - face_nodes_beg); const double flux_face = std::fabs(darcyflux_[face]); - double tof = 0.0; + face_term = 0.0; + cell_term_factor = 0.0; for (int ii = 0; ii < num_adj; ++ii) { const int f = adj_faces_[ii]; const double influx_f = (grid_.face_cells[2*f] == upwind_cell) ? -darcyflux_[f] : darcyflux_[f]; @@ -182,11 +246,13 @@ namespace Opm // const double omega = omega_star > 0.0 ? std::min(omega_star, 1.0) : 0.0; // SMU const double omega = omega_star > 0.0 ? omega_star/(1.0 + omega_star) : 0.0; - tof += (1.0 - omega)*tof_[upwind_cell] + omega*face_tof_[f]; + face_term += omega*face_tof_[f]; + cell_term_factor += (1.0 - omega); } - - // For now taking a simple average. - return tof/double(num_adj); + face_term /= double(num_adj); + cell_term_factor /= double(num_adj); } + + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 7a27e0ca..d9540abd 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -61,9 +61,11 @@ namespace Opm private: virtual void solveSingleCell(const int cell); + void solveSingleCellMultidimUpwind(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); - double multidimUpwindTof(const int face, const int upwind_cell) const; + void multidimUpwindTerms(const int face, const int upwind_cell, + double& face_term, double& cell_term_factor) const; private: const UnstructuredGrid& grid_;