From 29b2bcad454524b43865aeae9bc558eb62fee9cd 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 1/6] 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/compute_tof_from_files.cpp | 173 ++++++++++++++++++++++++++++ 1 file changed, 173 insertions(+) create mode 100644 examples/compute_tof_from_files.cpp diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp new file mode 100644 index 000000000..30b76f7b8 --- /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 43162d7d99255d66390d4fd82387d8ee6c78bc43 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 2/6] 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 30b76f7b8..012def4e6 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 bd013a7d4d6bdf0b8ec959ce86249fe3f612716d 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 3/6] Added parameters for controlling use of multidim upwinding. For now, you will simply get SPU even with use_multidim_upwind=true. --- .../reorder/TransportModelTracerTof.cpp | 17 +++++++++++++---- .../reorder/TransportModelTracerTof.hpp | 7 ++++++- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 8d93f2480..e5202a6d3 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 270ac328a..39b280d4b 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 1272575f1f405d8f7f648af309d081b68055fb0c 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 4/6] 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 012def4e6..17ad7b85b 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 a124d2e3be9942d9bfbb43338bdbe138939be528 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 5/6] 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 e5202a6d3..f4c2e72bb 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 39b280d4b..7a27e0ca0 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 7deba2cce0f53ac69706aab5da109c0387e9b02b 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 6/6] 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 f4c2e72bb..91b20ef01 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 7a27e0ca0..d9540abd2 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_;