From 88f4cd894d8eb4e179ff97c894c35c5be84d6814 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 17 Apr 2013 10:38:08 +0200 Subject: [PATCH 01/21] Bugfix: do not compute MDU terms for noflow faces. --- opm/core/tof/TofReorder.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 618b5417..9f7c289f 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -219,7 +219,7 @@ namespace Opm // Add flux to upwind_term or downwind_term_[face|cell_factor]. if (flux < 0.0) { upwind_term += flux*face_tof_[f]; - } else { + } else if (flux > 0.0) { double fterm, cterm_factor; multidimUpwindTerms(f, cell, fterm, cterm_factor); downwind_term_face += fterm*flux; From b7acc70ab76fb4c03b838ba90a18630821501ad8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 17 Apr 2013 12:58:15 +0200 Subject: [PATCH 02/21] Implement solveMultiCell() properly. Interface change: solver now requires a linear solver (for the multi-cell blocks only). Implementation uses new private method assembleSingleCell(), that is a modified copy of solveSingleCell(). Should refactor. --- examples/compute_tof.cpp | 2 +- examples/compute_tof_from_files.cpp | 5 +- opm/core/tof/TofReorder.cpp | 106 +++++++++++++++++++++++++++- opm/core/tof/TofReorder.hpp | 12 +++- 4 files changed, 121 insertions(+), 4 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 7abbe439..2efe285c 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -237,7 +237,7 @@ main(int argc, char** argv) if (use_dg) { dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); } else { - Opm::TofReorder tofsolver(*grid->c_grid(), use_multidim_upwind); + Opm::TofReorder tofsolver(*grid->c_grid(), linsolver, use_multidim_upwind); if (compute_tracer) { tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tof, tracer); } else { diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index deba61ea..2d716777 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -120,6 +120,9 @@ main(int argc, char** argv) } } + // Linear solver. + LinearSolverFactory linsolver(param); + // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); bool use_multidim_upwind = false; @@ -163,7 +166,7 @@ main(int argc, char** argv) if (use_dg) { dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof); } else { - Opm::TofReorder tofsolver(grid, use_multidim_upwind); + Opm::TofReorder tofsolver(grid, linsolver, use_multidim_upwind); if (compute_tracer) { tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tof, tracer); } else { diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 9f7c289f..092b8166 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -20,6 +20,7 @@ #include "config.h" #include #include +#include #include #include #include @@ -30,17 +31,21 @@ namespace Opm /// Construct solver. - /// \param[in] grid A 2d or 3d grid. + /// \param[in] gri d A 2d or 3d grid. + /// \param[in] linsolver Linear solver used for multi-cell blocks. /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. TofReorder::TofReorder(const UnstructuredGrid& grid, + const LinearSolverInterface& linsolver, const bool use_multidim_upwind) : grid_(grid), + linsolver_(linsolver), darcyflux_(0), porevolume_(0), source_(0), tof_(0), tracer_(0), num_tracers_(0), + block_index_(grid.number_of_cells, OutsideBlock), use_multidim_upwind_(use_multidim_upwind) { } @@ -195,6 +200,55 @@ namespace Opm + void TofReorder::assembleSingleCell(const int cell, + std::vector& local_column, + std::vector& local_coefficient, + double& rhs) + { + // 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_flux (note sign change resulting from + // different sign conventions: pos. source is injection, + // pos. flux is outflow). + double downwind_flux = std::max(-source_[cell], 0.0); + double upwind_term = 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 (flux < 0.0) { + // Using tof == 0 on inflow, so we only add a + // nonzero contribution if we are on an internal + // face. + if (other != -1) { + if (block_index_[other] == OutsideBlock) { + upwind_term += flux*tof_[other]; + } else { + local_column.push_back(block_index_[other]); + local_coefficient.push_back(flux); + } + } + } else { + downwind_flux += flux; + } + } + local_column.push_back(block_index_[cell]); + local_coefficient.push_back(downwind_flux); + rhs = porevolume_[cell] - upwind_term; + } + + + void TofReorder::solveSingleCellMultidimUpwind(const int cell) { @@ -247,10 +301,60 @@ namespace Opm void TofReorder::solveMultiCell(const int num_cells, const int* cells) { +#if 0 std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; for (int i = 0; i < num_cells; ++i) { solveSingleCell(cells[i]); } +#else + std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + // Give each cell a local index in this multi-cell block. + for (int ci = 0; ci < num_cells; ++ci) { + block_index_[cells[ci]] = ci; + } + + // Create sparse matrix and rhs vector for block equation system. + std::vector ia(num_cells+1); + std::vector ja; + std::vector sa; + ja.reserve(2*num_cells); + sa.reserve(2*num_cells); + std::vector rhs(num_cells); + std::vector< std::pair > rowdata; + for (int ci = 0; ci < num_cells; ++ci) { + ia[ci] = ja.size(); + assembleSingleCell(cells[ci], ja, sa, rhs[ci]); + // We standardize sparse row format: must sort row content by column index. + int num_row_elem = ja.size() - ia[ci]; + rowdata.resize(num_row_elem); + for (int i = 0; i < num_row_elem; ++i) { + rowdata[i].first = ja[ia[ci] + i]; + rowdata[i].second = sa[ia[ci] + i]; + } + std::sort(rowdata.begin(), rowdata.end()); + for (int i = 0; i < num_row_elem; ++i) { + ja[ia[ci] + i] = rowdata[i].first; + sa[ia[ci] + i] = rowdata[i].second; + } + } + ia.back() = ja.size(); + ASSERT(ja.size() == sa.size()); + + // Solve system. + std::vector tof_block(num_cells); + LinearSolverInterface::LinearSolverReport rep + = linsolver_.solve(num_cells, ja.size(), &ia[0], &ja[0], &sa[0], &rhs[0], &tof_block[0]); + if (!rep.converged) { + THROW("Multicell system with " << num_cells << " failed to converge."); + } + + // Write to global tof vector, and reset block indices to make + // it usable for next solveMultiCell() call. + for (int ci = 0; ci < num_cells; ++ci) { + tof_[cells[ci]] = tof_block[ci]; + block_index_[cells[ci]] = OutsideBlock; + } +#endif } diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index 2069767d..673e0ea1 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -30,6 +30,7 @@ namespace Opm { class IncompPropertiesInterface; + class LinearSolverInterface; /// Implements a first-order finite volume solver for /// (single-phase) time-of-flight using reordering. @@ -42,9 +43,11 @@ namespace Opm { public: /// Construct solver. - /// \param[in] grid A 2d or 3d grid. + /// \param[in] grid A 2d or 3d grid. + /// \param[in] linsolver Linear solver used for multi-cell blocks. /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. TofReorder(const UnstructuredGrid& grid, + const LinearSolverInterface& linsolver, const bool use_multidim_upwind = false); /// Solve for time-of-flight. @@ -79,6 +82,10 @@ namespace Opm private: virtual void solveSingleCell(const int cell); void solveSingleCellMultidimUpwind(const int cell); + void assembleSingleCell(const int cell, + std::vector& local_column, + std::vector& local_coefficient, + double& rhs); virtual void solveMultiCell(const int num_cells, const int* cells); void multidimUpwindTerms(const int face, const int upwind_cell, @@ -86,12 +93,15 @@ namespace Opm private: const UnstructuredGrid& grid_; + const LinearSolverInterface& linsolver_; 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_; double* tracer_; int num_tracers_; + enum { OutsideBlock = -1 }; + std::vector block_index_; // For solveMultiCell(). bool use_multidim_upwind_; std::vector face_tof_; // For multidim upwind face tofs. mutable std::vector adj_faces_; // For multidim upwind logic. From 3d3173d04f5132596d8692a8d6f13fcfb57b159d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 17 Apr 2013 13:18:22 +0200 Subject: [PATCH 03/21] More refined and concise output of block info. --- opm/core/tof/TofReorder.cpp | 8 +++++++- opm/core/tof/TofReorder.hpp | 2 ++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 092b8166..49d57a83 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -83,7 +83,11 @@ namespace Opm std::fill(face_tof_.begin(), face_tof_.end(), 0.0); } num_tracers_ = 0; + num_multicell_ = 0; + max_size_multicell_ = 0; reorderAndTransport(grid_, darcyflux); + std::cout << num_multicell_ << " multicell blocks with max size " + << max_size_multicell_ << " cells." << std::endl; } @@ -307,7 +311,9 @@ namespace Opm solveSingleCell(cells[i]); } #else - std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + // std::cout << "Solving multi-cell dependent equation with " << num_cells << " cells." << std::endl; + ++num_multicell_; + max_size_multicell_ = std::max(max_size_multicell_, num_cells); // Give each cell a local index in this multi-cell block. for (int ci = 0; ci < num_cells; ++ci) { block_index_[cells[ci]] = ci; diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index 673e0ea1..fe4705f7 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -102,6 +102,8 @@ namespace Opm int num_tracers_; enum { OutsideBlock = -1 }; std::vector block_index_; // For solveMultiCell(). + int num_multicell_; + int max_size_multicell_; bool use_multidim_upwind_; std::vector face_tof_; // For multidim upwind face tofs. mutable std::vector adj_faces_; // For multidim upwind logic. From 7c172fe1928349fff4db0c3d22019c72ec5c6621 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 22 Apr 2013 11:07:41 +0200 Subject: [PATCH 04/21] Make multidim upwind more forgiving of bad grids. Now it will no longer trigger assertation failure when grids are not edge-conformal (faulted cornerpoint grids processed by our code will usually not be). Minor algorithm change to handle this. --- opm/core/tof/TofReorder.cpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 49d57a83..88621661 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -387,6 +387,9 @@ namespace Opm // 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. + // Note: Modified algorithm to consider faces that share even a single vertex with + // the input face. This reduces the problem of non-edge-conformal grids, but does not + // eliminate it entirely. // Identify the adjacent faces of the upwind cell. const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[face]; @@ -404,11 +407,11 @@ namespace Opm 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). + // Before: neighbours over an edge (3d) or vertex (2d). + // Now: neighbours across a vertex. + // if (num_common == grid_.dimensions - 1) { + if (num_common > 0) { adj_faces_.push_back(f); - } else { - ASSERT(num_common == 0); } } } @@ -416,7 +419,9 @@ namespace Opm // 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); + // The assertion below only holds if the grid is edge-conformal. + // No longer testing, since method no longer requires it. + // ASSERT(num_adj == face_nodes_end - face_nodes_beg); const double flux_face = std::fabs(darcyflux_[face]); face_term = 0.0; cell_term_factor = 0.0; From 489c41d60dc3aa30d1c98f102f5882f77fdbeac3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 22 Apr 2013 11:11:55 +0200 Subject: [PATCH 05/21] Optimize multi-cell solve and add new Gauss-Seidel variant. --- opm/core/tof/TofReorder.cpp | 80 +++++++++++++++++++++++-------------- opm/core/tof/TofReorder.hpp | 10 ++++- 2 files changed, 60 insertions(+), 30 deletions(-) diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 88621661..82862a75 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -150,6 +150,7 @@ namespace Opm void TofReorder::solveSingleCell(const int cell) { +#if 1 if (use_multidim_upwind_) { solveSingleCellMultidimUpwind(cell); return; @@ -189,7 +190,6 @@ namespace Opm downwind_flux += flux; } } - // Compute tof. tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux; @@ -200,6 +200,18 @@ namespace Opm tracer_[num_tracers_*cell + tr] *= -1.0/downwind_flux; } } +#else + ja_.clear(); + sa_.clear(); + block_index_[cell] = 0; + double rhs = 0.0; + assembleSingleCell(cell, ja_, sa_, rhs); + ASSERT(sa_.size() == 1); + ASSERT(ja_.size() == sa_.size()); + ASSERT(ja_[0] == 0); + tof_[cell] = rhs/sa_[0]; + block_index_[cell] = OutsideBlock; +#endif } @@ -305,51 +317,48 @@ namespace Opm void TofReorder::solveMultiCell(const int num_cells, const int* cells) { -#if 0 - std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; - for (int i = 0; i < num_cells; ++i) { - solveSingleCell(cells[i]); - } -#else - // std::cout << "Solving multi-cell dependent equation with " << num_cells << " cells." << std::endl; ++num_multicell_; max_size_multicell_ = std::max(max_size_multicell_, num_cells); +#if 0 + if (use_multidim_upwind_) { + THROW("Encountered multi-cell dependent block. " + "Block solve not implemented for multidimensional upwind method"); + } + // std::cout << "Solving multi-cell dependent equation with " << num_cells << " cells." << std::endl; // Give each cell a local index in this multi-cell block. for (int ci = 0; ci < num_cells; ++ci) { block_index_[cells[ci]] = ci; } // Create sparse matrix and rhs vector for block equation system. - std::vector ia(num_cells+1); - std::vector ja; - std::vector sa; - ja.reserve(2*num_cells); - sa.reserve(2*num_cells); - std::vector rhs(num_cells); - std::vector< std::pair > rowdata; + ia_.clear(); + ja_.clear(); + sa_.clear(); + rhs_.resize(num_cells); + rowdata_.clear(); for (int ci = 0; ci < num_cells; ++ci) { - ia[ci] = ja.size(); - assembleSingleCell(cells[ci], ja, sa, rhs[ci]); + ia_.push_back(ja_.size()); + assembleSingleCell(cells[ci], ja_, sa_, rhs_[ci]); // We standardize sparse row format: must sort row content by column index. - int num_row_elem = ja.size() - ia[ci]; - rowdata.resize(num_row_elem); + int num_row_elem = ja_.size() - ia_[ci]; + rowdata_.resize(num_row_elem); for (int i = 0; i < num_row_elem; ++i) { - rowdata[i].first = ja[ia[ci] + i]; - rowdata[i].second = sa[ia[ci] + i]; + rowdata_[i].first = ja_[ia_[ci] + i]; + rowdata_[i].second = sa_[ia_[ci] + i]; } - std::sort(rowdata.begin(), rowdata.end()); + std::sort(rowdata_.begin(), rowdata_.end()); for (int i = 0; i < num_row_elem; ++i) { - ja[ia[ci] + i] = rowdata[i].first; - sa[ia[ci] + i] = rowdata[i].second; + ja_[ia_[ci] + i] = rowdata_[i].first; + sa_[ia_[ci] + i] = rowdata_[i].second; } } - ia.back() = ja.size(); - ASSERT(ja.size() == sa.size()); + ia_.push_back(ja_.size()); + ASSERT(ja_.size() == sa_.size()); // Solve system. - std::vector tof_block(num_cells); + tof_block_.resize(num_cells); LinearSolverInterface::LinearSolverReport rep - = linsolver_.solve(num_cells, ja.size(), &ia[0], &ja[0], &sa[0], &rhs[0], &tof_block[0]); + = linsolver_.solve(num_cells, ja_.size(), &ia_[0], &ja_[0], &sa_[0], &rhs_[0], &tof_block_[0]); if (!rep.converged) { THROW("Multicell system with " << num_cells << " failed to converge."); } @@ -357,9 +366,22 @@ namespace Opm // Write to global tof vector, and reset block indices to make // it usable for next solveMultiCell() call. for (int ci = 0; ci < num_cells; ++ci) { - tof_[cells[ci]] = tof_block[ci]; + tof_[cells[ci]] = tof_block_[ci]; block_index_[cells[ci]] = OutsideBlock; } +#else + // std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl; + double max_delta = 1e100; + while (max_delta > 1e-3) { + max_delta = 0.0; + for (int ci = 0; ci < num_cells; ++ci) { + const int cell = cells[ci]; + const double tof_before = tof_[cell]; + solveSingleCell(cell); + max_delta = std::max(max_delta, std::fabs(tof_[cell] - tof_before)); + } + // std::cout << "Max delta = " << max_delta << std::endl; + } #endif } diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index fe4705f7..c88121d2 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -100,10 +100,18 @@ namespace Opm double* tof_; double* tracer_; int num_tracers_; + // For solveMultiCell(): enum { OutsideBlock = -1 }; - std::vector block_index_; // For solveMultiCell(). + std::vector block_index_; int num_multicell_; int max_size_multicell_; + std::vector ia_; + std::vector ja_; + std::vector sa_; + std::vector< std::pair > rowdata_; + std::vector rhs_; + std::vector tof_block_; + // For multidim upwinding: bool use_multidim_upwind_; std::vector face_tof_; // For multidim upwind face tofs. mutable std::vector adj_faces_; // For multidim upwind logic. From cf38c91f742fa6e91b194219739b605c7a981b3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 22 Apr 2013 11:22:23 +0200 Subject: [PATCH 06/21] Removed experimental multi-cell solver code. Since the Gauss-Seidel approach seems to be both simplest and fastest, all parts dealing with assembling multicell systems have been removed. --- examples/compute_tof.cpp | 2 +- examples/compute_tof_from_files.cpp | 5 +- opm/core/tof/TofReorder.cpp | 138 +++------------------------- opm/core/tof/TofReorder.hpp | 14 +-- 4 files changed, 18 insertions(+), 141 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 2efe285c..7abbe439 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -237,7 +237,7 @@ main(int argc, char** argv) if (use_dg) { dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); } else { - Opm::TofReorder tofsolver(*grid->c_grid(), linsolver, use_multidim_upwind); + Opm::TofReorder tofsolver(*grid->c_grid(), use_multidim_upwind); if (compute_tracer) { tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tof, tracer); } else { diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index 2d716777..deba61ea 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -120,9 +120,6 @@ main(int argc, char** argv) } } - // Linear solver. - LinearSolverFactory linsolver(param); - // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); bool use_multidim_upwind = false; @@ -166,7 +163,7 @@ main(int argc, char** argv) if (use_dg) { dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof); } else { - Opm::TofReorder tofsolver(grid, linsolver, use_multidim_upwind); + Opm::TofReorder tofsolver(grid, use_multidim_upwind); if (compute_tracer) { tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tof, tracer); } else { diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 82862a75..09f754e9 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -20,7 +20,6 @@ #include "config.h" #include #include -#include #include #include #include @@ -32,20 +31,17 @@ namespace Opm /// Construct solver. /// \param[in] gri d A 2d or 3d grid. - /// \param[in] linsolver Linear solver used for multi-cell blocks. /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. TofReorder::TofReorder(const UnstructuredGrid& grid, - const LinearSolverInterface& linsolver, const bool use_multidim_upwind) : grid_(grid), - linsolver_(linsolver), darcyflux_(0), porevolume_(0), source_(0), tof_(0), tracer_(0), num_tracers_(0), - block_index_(grid.number_of_cells, OutsideBlock), + gauss_seidel_tol_(1e-3), use_multidim_upwind_(use_multidim_upwind) { } @@ -85,9 +81,13 @@ namespace Opm num_tracers_ = 0; num_multicell_ = 0; max_size_multicell_ = 0; + max_iter_multicell_ = 0; reorderAndTransport(grid_, darcyflux); - std::cout << num_multicell_ << " multicell blocks with max size " - << max_size_multicell_ << " cells." << std::endl; + if (num_multicell_ > 0) { + std::cout << num_multicell_ << " multicell blocks with max size " + << max_size_multicell_ << " cells in upto " + << max_iter_multicell_ << " iterations." << std::endl; + } } @@ -150,7 +150,6 @@ namespace Opm void TofReorder::solveSingleCell(const int cell) { -#if 1 if (use_multidim_upwind_) { solveSingleCellMultidimUpwind(cell); return; @@ -200,72 +199,10 @@ namespace Opm tracer_[num_tracers_*cell + tr] *= -1.0/downwind_flux; } } -#else - ja_.clear(); - sa_.clear(); - block_index_[cell] = 0; - double rhs = 0.0; - assembleSingleCell(cell, ja_, sa_, rhs); - ASSERT(sa_.size() == 1); - ASSERT(ja_.size() == sa_.size()); - ASSERT(ja_[0] == 0); - tof_[cell] = rhs/sa_[0]; - block_index_[cell] = OutsideBlock; -#endif } - void TofReorder::assembleSingleCell(const int cell, - std::vector& local_column, - std::vector& local_coefficient, - double& rhs) - { - // 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_flux (note sign change resulting from - // different sign conventions: pos. source is injection, - // pos. flux is outflow). - double downwind_flux = std::max(-source_[cell], 0.0); - double upwind_term = 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 (flux < 0.0) { - // Using tof == 0 on inflow, so we only add a - // nonzero contribution if we are on an internal - // face. - if (other != -1) { - if (block_index_[other] == OutsideBlock) { - upwind_term += flux*tof_[other]; - } else { - local_column.push_back(block_index_[other]); - local_coefficient.push_back(flux); - } - } - } else { - downwind_flux += flux; - } - } - local_column.push_back(block_index_[cell]); - local_coefficient.push_back(downwind_flux); - rhs = porevolume_[cell] - upwind_term; - } - - - - void TofReorder::solveSingleCellMultidimUpwind(const int cell) { // Compute flux terms. @@ -298,7 +235,7 @@ namespace Opm } // Compute tof for cell. - tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor; // } + 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) { @@ -319,61 +256,14 @@ namespace Opm { ++num_multicell_; max_size_multicell_ = std::max(max_size_multicell_, num_cells); -#if 0 - if (use_multidim_upwind_) { - THROW("Encountered multi-cell dependent block. " - "Block solve not implemented for multidimensional upwind method"); - } - // std::cout << "Solving multi-cell dependent equation with " << num_cells << " cells." << std::endl; - // Give each cell a local index in this multi-cell block. - for (int ci = 0; ci < num_cells; ++ci) { - block_index_[cells[ci]] = ci; - } - - // Create sparse matrix and rhs vector for block equation system. - ia_.clear(); - ja_.clear(); - sa_.clear(); - rhs_.resize(num_cells); - rowdata_.clear(); - for (int ci = 0; ci < num_cells; ++ci) { - ia_.push_back(ja_.size()); - assembleSingleCell(cells[ci], ja_, sa_, rhs_[ci]); - // We standardize sparse row format: must sort row content by column index. - int num_row_elem = ja_.size() - ia_[ci]; - rowdata_.resize(num_row_elem); - for (int i = 0; i < num_row_elem; ++i) { - rowdata_[i].first = ja_[ia_[ci] + i]; - rowdata_[i].second = sa_[ia_[ci] + i]; - } - std::sort(rowdata_.begin(), rowdata_.end()); - for (int i = 0; i < num_row_elem; ++i) { - ja_[ia_[ci] + i] = rowdata_[i].first; - sa_[ia_[ci] + i] = rowdata_[i].second; - } - } - ia_.push_back(ja_.size()); - ASSERT(ja_.size() == sa_.size()); - - // Solve system. - tof_block_.resize(num_cells); - LinearSolverInterface::LinearSolverReport rep - = linsolver_.solve(num_cells, ja_.size(), &ia_[0], &ja_[0], &sa_[0], &rhs_[0], &tof_block_[0]); - if (!rep.converged) { - THROW("Multicell system with " << num_cells << " failed to converge."); - } - - // Write to global tof vector, and reset block indices to make - // it usable for next solveMultiCell() call. - for (int ci = 0; ci < num_cells; ++ci) { - tof_[cells[ci]] = tof_block_[ci]; - block_index_[cells[ci]] = OutsideBlock; - } -#else // std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl; + + // Using a Gauss-Seidel approach. double max_delta = 1e100; - while (max_delta > 1e-3) { + int num_iter = 0; + while (max_delta > gauss_seidel_tol_) { max_delta = 0.0; + ++num_iter; for (int ci = 0; ci < num_cells; ++ci) { const int cell = cells[ci]; const double tof_before = tof_[cell]; @@ -382,7 +272,7 @@ namespace Opm } // std::cout << "Max delta = " << max_delta << std::endl; } -#endif + max_iter_multicell_ = std::max(max_iter_multicell_, num_iter); } diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index c88121d2..57a82cdb 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -30,7 +30,6 @@ namespace Opm { class IncompPropertiesInterface; - class LinearSolverInterface; /// Implements a first-order finite volume solver for /// (single-phase) time-of-flight using reordering. @@ -44,10 +43,8 @@ namespace Opm public: /// Construct solver. /// \param[in] grid A 2d or 3d grid. - /// \param[in] linsolver Linear solver used for multi-cell blocks. /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. TofReorder(const UnstructuredGrid& grid, - const LinearSolverInterface& linsolver, const bool use_multidim_upwind = false); /// Solve for time-of-flight. @@ -93,7 +90,6 @@ namespace Opm private: const UnstructuredGrid& grid_; - const LinearSolverInterface& linsolver_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell @@ -101,16 +97,10 @@ namespace Opm double* tracer_; int num_tracers_; // For solveMultiCell(): - enum { OutsideBlock = -1 }; - std::vector block_index_; + double gauss_seidel_tol_; int num_multicell_; int max_size_multicell_; - std::vector ia_; - std::vector ja_; - std::vector sa_; - std::vector< std::pair > rowdata_; - std::vector rhs_; - std::vector tof_block_; + int max_iter_multicell_; // For multidim upwinding: bool use_multidim_upwind_; std::vector face_tof_; // For multidim upwind face tofs. From a2125fee09bbd35ab14b52d0dcfcaa1824f87d99 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 22 Apr 2013 11:29:59 +0200 Subject: [PATCH 07/21] Fix typo and minimize formatting changes. --- opm/core/tof/TofReorder.cpp | 4 +++- opm/core/tof/TofReorder.hpp | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 09f754e9..1091fbdd 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -30,7 +30,7 @@ namespace Opm /// Construct solver. - /// \param[in] gri d A 2d or 3d grid. + /// \param[in] grid A 2d or 3d grid. /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. TofReorder::TofReorder(const UnstructuredGrid& grid, const bool use_multidim_upwind) @@ -189,6 +189,7 @@ namespace Opm downwind_flux += flux; } } + // Compute tof. tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux; @@ -203,6 +204,7 @@ namespace Opm + void TofReorder::solveSingleCellMultidimUpwind(const int cell) { // Compute flux terms. diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index 57a82cdb..56ccd112 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -42,7 +42,7 @@ namespace Opm { public: /// Construct solver. - /// \param[in] grid A 2d or 3d grid. + /// \param[in] grid A 2d or 3d grid. /// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding. TofReorder(const UnstructuredGrid& grid, const bool use_multidim_upwind = false); From 8b6faffe008b077adfdae5a803702d6df19ed485 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 22 Apr 2013 14:02:45 +0200 Subject: [PATCH 08/21] Work in progress on tracers. - Changed interface. - Read tracerheads (tracer start locations) from file in compute_tof_from_files. - Initialize tracerheads from wells in compute_tof. --- examples/compute_tof.cpp | 24 +++++++++++++++++++++- examples/compute_tof_from_files.cpp | 21 +++++++++++++++++-- opm/core/tof/TofReorder.cpp | 31 +++++++++++++++++++---------- opm/core/tof/TofReorder.hpp | 5 +++++ 4 files changed, 68 insertions(+), 13 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 7abbe439..1d71628a 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -30,6 +30,7 @@ #include #include #include +#include #include #include #include @@ -65,6 +66,25 @@ namespace std::cout << "----------------------------------------------------------------" << std::endl; } } + + void buildTracerheadsFromWells(const Wells* wells, + const Opm::WellState& well_state, + Opm::SparseTable& tracerheads) + { + if (wells == 0) { + return; + } + tracerheads.clear(); + const int num_wells = wells->number_of_wells; + for (int w = 0; w < num_wells; ++w) { + if (wells->type[w] != INJECTOR) { + continue; + } + tracerheads.appendRow(wells->well_cells + wells->well_connpos[w], + wells->well_cells + wells->well_connpos[w + 1]); + } + } + } // anon namespace @@ -239,7 +259,9 @@ main(int argc, char** argv) } else { Opm::TofReorder tofsolver(*grid->c_grid(), use_multidim_upwind); if (compute_tracer) { - tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tof, tracer); + Opm::SparseTable tracerheads; + buildTracerheadsFromWells(wells->c_wells(), well_state, tracerheads); + tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tracerheads, tof, tracer); } else { tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); } diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index deba61ea..a186e0f4 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -120,6 +121,23 @@ main(int argc, char** argv) } } + const bool compute_tracer = param.getDefault("compute_tracer", false); + Opm::SparseTable tracerheads; + { + std::ifstream tr_stream(param.get("tracerheads_filename").c_str()); + int num_rows; + tr_stream >> num_rows; + for (int row = 0; row < num_rows; ++row) { + int row_size; + tr_stream >> row_size; + std::vector rowdata(row_size); + for (int elem = 0; elem < row_size; ++elem) { + tr_stream >> rowdata[elem]; + } + tracerheads.appendRow(rowdata.begin(), rowdata.end()); + } + } + // Choice of tof solver. bool use_dg = param.getDefault("use_dg", false); bool use_multidim_upwind = false; @@ -130,7 +148,6 @@ main(int argc, char** argv) } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } - bool compute_tracer = param.getDefault("compute_tracer", false); if (use_dg && compute_tracer) { THROW("DG for tracer not yet implemented."); } @@ -165,7 +182,7 @@ main(int argc, char** argv) } else { Opm::TofReorder tofsolver(grid, use_multidim_upwind); if (compute_tracer) { - tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tof, tracer); + tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tracerheads, tof, tracer); } else { tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof); } diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 1091fbdd..ccb35ed4 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #include #include @@ -57,9 +58,9 @@ namespace Opm /// (-) outflow flux. /// \param[out] tof Array of time-of-flight values. void TofReorder::solveTof(const double* darcyflux, - const double* porevolume, - const double* source, - std::vector& tof) + const double* porevolume, + const double* source, + std::vector& tof) { darcyflux_ = darcyflux; porevolume_ = porevolume; @@ -101,12 +102,15 @@ namespace Opm /// \param[in] source Source term. Sign convention is: /// (+) inflow flux, /// (-) outflow flux. + /// \param[in] tracerheads Table containing one row per tracer, and each + /// row contains the source cells for that tracer. /// \param[out] tof Array of time-of-flight values (1 per cell). /// \param[out] tracer Array of tracer values (N per cell, where N is /// the number of cells c for which source[c] > 0.0). void TofReorder::solveTofTracer(const double* darcyflux, const double* porevolume, const double* source, + const SparseTable& tracerheads, std::vector& tof, std::vector& tracer) { @@ -123,26 +127,33 @@ namespace Opm tof.resize(grid_.number_of_cells); std::fill(tof.begin(), tof.end(), 0.0); tof_ = &tof[0]; + // Find the tracer heads (injectors). - std::vector tracerheads; - for (int c = 0; c < grid_.number_of_cells; ++c) { - if (source[c] > 0.0) { - tracerheads.push_back(c); - } - } num_tracers_ = tracerheads.size(); tracer.resize(grid_.number_of_cells*num_tracers_); std::fill(tracer.begin(), tracer.end(), 0.0); for (int tr = 0; tr < num_tracers_; ++tr) { - tracer[tracerheads[tr]*num_tracers_ + tr] = 1.0; + for (int i = 0; i < tracerheads[tr].size(); ++i) { + const int cell = tracerheads[tr][i]; + tracer[cell*num_tracers_ + tr] = 1.0; + } } + tracer_ = &tracer[0]; if (use_multidim_upwind_) { face_tof_.resize(grid_.number_of_faces); std::fill(face_tof_.begin(), face_tof_.end(), 0.0); THROW("Multidimensional upwind not yet implemented for tracer."); } + num_multicell_ = 0; + max_size_multicell_ = 0; + max_iter_multicell_ = 0; reorderAndTransport(grid_, darcyflux); + if (num_multicell_ > 0) { + std::cout << num_multicell_ << " multicell blocks with max size " + << max_size_multicell_ << " cells in upto " + << max_iter_multicell_ << " iterations." << std::endl; + } } diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index 56ccd112..7394ad59 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -24,12 +24,14 @@ #include #include #include + struct UnstructuredGrid; namespace Opm { class IncompPropertiesInterface; + template class SparseTable; /// Implements a first-order finite volume solver for /// (single-phase) time-of-flight using reordering. @@ -67,12 +69,15 @@ namespace Opm /// \param[in] source Source term. Sign convention is: /// (+) inflow flux, /// (-) outflow flux. + /// \param[in] tracerheads Table containing one row per tracer, and each + /// row contains the source cells for that tracer. /// \param[out] tof Array of time-of-flight values (1 per cell). /// \param[out] tracer Array of tracer values (N per cell, where N is /// the number of cells c for which source[c] > 0.0). void solveTofTracer(const double* darcyflux, const double* porevolume, const double* source, + const SparseTable& tracerheads, std::vector& tof, std::vector& tracer); From fd4ac88334c1caf35b9198019f15609b332fc3c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 23 Apr 2013 09:54:52 +0200 Subject: [PATCH 09/21] Minor fix: do not read tracer heads unless compute_tracer is true. --- examples/compute_tof_from_files.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index a186e0f4..723df378 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -123,7 +123,7 @@ main(int argc, char** argv) const bool compute_tracer = param.getDefault("compute_tracer", false); Opm::SparseTable tracerheads; - { + if (compute_tracer) { std::ifstream tr_stream(param.get("tracerheads_filename").c_str()); int num_rows; tr_stream >> num_rows; From 2fbf3bd7cb6b378215e249ec64f9d26f5a81956a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 23 Apr 2013 13:38:47 +0200 Subject: [PATCH 10/21] Implement solveMultiCell() with Gauss-Seidel. --- opm/core/tof/TofDiscGalReorder.cpp | 40 ++++++++++++++++++++++++++---- opm/core/tof/TofDiscGalReorder.hpp | 6 +++++ 2 files changed, 41 insertions(+), 5 deletions(-) diff --git a/opm/core/tof/TofDiscGalReorder.cpp b/opm/core/tof/TofDiscGalReorder.cpp index ef50b665..cc348f58 100644 --- a/opm/core/tof/TofDiscGalReorder.cpp +++ b/opm/core/tof/TofDiscGalReorder.cpp @@ -45,7 +45,8 @@ namespace Opm limiter_method_(MinUpwindAverage), limiter_usage_(DuringComputations), coord_(grid.dimensions), - velocity_(grid.dimensions) + velocity_(grid.dimensions), + gauss_seidel_tol_(1e-3) { const int dg_degree = param.getDefault("dg_degree", 0); const bool use_tensorial_basis = param.getDefault("use_tensorial_basis", false); @@ -123,6 +124,10 @@ namespace Opm basis_nb_.resize(num_basis); grad_basis_.resize(num_basis*grid_.dimensions); velocity_interpolation_->setupFluxes(darcyflux); + num_multicell_ = 0; + max_size_multicell_ = 0; + max_iter_multicell_ = 0; + num_singlesolves_ = 0; reorderAndTransport(grid_, darcyflux); switch (limiter_usage_) { case AsPostProcess: @@ -137,6 +142,13 @@ namespace Opm default: THROW("Unknown limiter usage choice: " << limiter_usage_); } + if (num_multicell_ > 0) { + std::cout << num_multicell_ << " multicell blocks with max size " + << max_size_multicell_ << " cells in upto " + << max_iter_multicell_ << " iterations." << std::endl; + std::cout << "Average solves per cell (for all cells) was " + << double(num_singlesolves_)/double(grid_.number_of_cells) << std::endl; + } } @@ -156,6 +168,7 @@ namespace Opm const int dim = grid_.dimensions; const int num_basis = basis_func_->numBasisFunc(); + ++num_singlesolves_; std::fill(rhs_.begin(), rhs_.end(), 0.0); std::fill(jac_.begin(), jac_.end(), 0.0); @@ -367,10 +380,27 @@ namespace Opm void TofDiscGalReorder::solveMultiCell(const int num_cells, const int* cells) { - std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; - for (int i = 0; i < num_cells; ++i) { - solveSingleCell(cells[i]); + ++num_multicell_; + max_size_multicell_ = std::max(max_size_multicell_, num_cells); + // std::cout << "Multiblock solve with " << num_cells << " cells." << std::endl; + + // Using a Gauss-Seidel approach. + const int nb = basis_func_->numBasisFunc(); + double max_delta = 1e100; + int num_iter = 0; + while (max_delta > gauss_seidel_tol_) { + max_delta = 0.0; + ++num_iter; + for (int ci = 0; ci < num_cells; ++ci) { + const int cell = cells[ci]; + const double tof_before = basis_func_->functionAverage(&tof_coeff_[nb*cell]); + solveSingleCell(cell); + const double tof_after = basis_func_->functionAverage(&tof_coeff_[nb*cell]); + max_delta = std::max(max_delta, std::fabs(tof_after - tof_before)); + } + // std::cout << "Max delta = " << max_delta << std::endl; } + max_iter_multicell_ = std::max(max_iter_multicell_, num_iter); } @@ -462,7 +492,7 @@ namespace Opm double limiter = (tof_c - min_upstream_tof)/(tof_c - min_here_tof); if (tof_c < min_upstream_tof) { // Handle by setting a flat solution. - std::cout << "Trouble in cell " << cell << std::endl; + // std::cout << "Trouble in cell " << cell << std::endl; limiter = 0.0; basis_func_->addConstant(min_upstream_tof - tof_c, tof + num_basis*cell); } diff --git a/opm/core/tof/TofDiscGalReorder.hpp b/opm/core/tof/TofDiscGalReorder.hpp index 079051a2..b92d2cce 100644 --- a/opm/core/tof/TofDiscGalReorder.hpp +++ b/opm/core/tof/TofDiscGalReorder.hpp @@ -125,6 +125,12 @@ namespace Opm mutable std::vector basis_nb_; std::vector grad_basis_; std::vector velocity_; + int num_singlesolves_; + // Used by solveMultiCell(): + double gauss_seidel_tol_; + int num_multicell_; + int max_size_multicell_; + int max_iter_multicell_; // Private methods From 8552347080b25d5ebedbe7fd61cb14145d2c6c0a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 23 Apr 2013 15:35:47 +0200 Subject: [PATCH 11/21] Tracer computations are now fixed and robust. - Handles noisy source terms. - Works with repeated solve calls (multi-cell block solves). --- opm/core/tof/TofReorder.cpp | 16 +++++++++++++--- opm/core/tof/TofReorder.hpp | 2 ++ 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index ccb35ed4..0393dda1 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -132,10 +132,13 @@ namespace Opm num_tracers_ = tracerheads.size(); tracer.resize(grid_.number_of_cells*num_tracers_); std::fill(tracer.begin(), tracer.end(), 0.0); + tracerhead_by_cell_.clear(); + tracerhead_by_cell_.resize(grid_.number_of_cells, NoTracerHead); for (int tr = 0; tr < num_tracers_; ++tr) { for (int i = 0; i < tracerheads[tr].size(); ++i) { const int cell = tracerheads[tr][i]; tracer[cell*num_tracers_ + tr] = 1.0; + tracerhead_by_cell_[cell] = tr; } } @@ -171,6 +174,11 @@ namespace Opm // to the downwind_flux (note sign change resulting from // different sign conventions: pos. source is injection, // pos. flux is outflow). + if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) { + for (int tr = 0; tr < num_tracers_; ++tr) { + tracer_[num_tracers_*cell + tr] = 0.0; + } + } double upwind_term = 0.0; double downwind_flux = std::max(-source_[cell], 0.0); for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { @@ -192,8 +200,10 @@ namespace Opm // face. if (other != -1) { upwind_term += flux*tof_[other]; - for (int tr = 0; tr < num_tracers_; ++tr) { - tracer_[num_tracers_*cell + tr] += flux*tracer_[num_tracers_*other + tr]; + if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) { + for (int tr = 0; tr < num_tracers_; ++tr) { + tracer_[num_tracers_*cell + tr] += flux*tracer_[num_tracers_*other + tr]; + } } } } else { @@ -206,7 +216,7 @@ namespace Opm // Compute tracers (if any). // Do not change tracer solution in source cells. - if (source_[cell] <= 0.0) { + if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) { for (int tr = 0; tr < num_tracers_; ++tr) { tracer_[num_tracers_*cell + tr] *= -1.0/downwind_flux; } diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index 7394ad59..42b8f7b9 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -101,6 +101,8 @@ namespace Opm double* tof_; double* tracer_; int num_tracers_; + enum { NoTracerHead = -1 }; + std::vector tracerhead_by_cell_; // For solveMultiCell(): double gauss_seidel_tol_; int num_multicell_; From 3086d5d815ad58c83b73a658fc934b5e01a7092e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 24 Apr 2013 10:36:52 +0200 Subject: [PATCH 12/21] Update doc, only initialize tracerhead_by_cell_ if needed. --- opm/core/tof/TofReorder.cpp | 12 ++++++------ opm/core/tof/TofReorder.hpp | 6 ++---- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 0393dda1..8907d308 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -95,8 +95,6 @@ namespace Opm /// Solve for time-of-flight and a number of tracers. - /// One tracer will be used for each inflow flux specified in - /// the source parameter. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. /// \param[in] source Source term. Sign convention is: @@ -105,8 +103,8 @@ namespace Opm /// \param[in] tracerheads Table containing one row per tracer, and each /// row contains the source cells for that tracer. /// \param[out] tof Array of time-of-flight values (1 per cell). - /// \param[out] tracer Array of tracer values (N per cell, where N is - /// the number of cells c for which source[c] > 0.0). + /// \param[out] tracer Array of tracer values. N per cell, where N is + /// equalt to tracerheads.size(). void TofReorder::solveTofTracer(const double* darcyflux, const double* porevolume, const double* source, @@ -132,8 +130,10 @@ namespace Opm num_tracers_ = tracerheads.size(); tracer.resize(grid_.number_of_cells*num_tracers_); std::fill(tracer.begin(), tracer.end(), 0.0); - tracerhead_by_cell_.clear(); - tracerhead_by_cell_.resize(grid_.number_of_cells, NoTracerHead); + if (num_tracers_ > 0) { + tracerhead_by_cell_.clear(); + tracerhead_by_cell_.resize(grid_.number_of_cells, NoTracerHead); + } for (int tr = 0; tr < num_tracers_; ++tr) { for (int i = 0; i < tracerheads[tr].size(); ++i) { const int cell = tracerheads[tr][i]; diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index 42b8f7b9..4abd640a 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -62,8 +62,6 @@ namespace Opm std::vector& tof); /// Solve for time-of-flight and a number of tracers. - /// One tracer will be used for each inflow flux specified in - /// the source parameter. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. /// \param[in] source Source term. Sign convention is: @@ -72,8 +70,8 @@ namespace Opm /// \param[in] tracerheads Table containing one row per tracer, and each /// row contains the source cells for that tracer. /// \param[out] tof Array of time-of-flight values (1 per cell). - /// \param[out] tracer Array of tracer values (N per cell, where N is - /// the number of cells c for which source[c] > 0.0). + /// \param[out] tracer Array of tracer values. N per cell, where N is + /// equalt to tracerheads.size(). void solveTofTracer(const double* darcyflux, const double* porevolume, const double* source, From 1a84b4fe7d9383f67bb0e40ea50d1e64cf768e83 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 24 Apr 2013 10:39:50 +0200 Subject: [PATCH 13/21] Add tracer computations (method solveTofTracer()). Same interface as in class TofReorder. --- opm/core/tof/TofDiscGalReorder.cpp | 141 +++++++++++++++++++++++++++-- opm/core/tof/TofDiscGalReorder.hpp | 39 +++++++- 2 files changed, 170 insertions(+), 10 deletions(-) diff --git a/opm/core/tof/TofDiscGalReorder.cpp b/opm/core/tof/TofDiscGalReorder.cpp index cc348f58..e2e4839a 100644 --- a/opm/core/tof/TofDiscGalReorder.cpp +++ b/opm/core/tof/TofDiscGalReorder.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -98,9 +99,9 @@ namespace Opm /// Solve for time-of-flight. void TofDiscGalReorder::solveTof(const double* darcyflux, - const double* porevolume, - const double* source, - std::vector& tof_coeff) + const double* porevolume, + const double* source, + std::vector& tof_coeff) { darcyflux_ = darcyflux; porevolume_ = porevolume; @@ -124,6 +125,103 @@ namespace Opm basis_nb_.resize(num_basis); grad_basis_.resize(num_basis*grid_.dimensions); velocity_interpolation_->setupFluxes(darcyflux); + num_tracers_ = 0; + num_multicell_ = 0; + max_size_multicell_ = 0; + max_iter_multicell_ = 0; + num_singlesolves_ = 0; + reorderAndTransport(grid_, darcyflux); + switch (limiter_usage_) { + case AsPostProcess: + applyLimiterAsPostProcess(); + break; + case AsSimultaneousPostProcess: + applyLimiterAsSimultaneousPostProcess(); + break; + case DuringComputations: + // Do nothing. + break; + default: + THROW("Unknown limiter usage choice: " << limiter_usage_); + } + if (num_multicell_ > 0) { + std::cout << num_multicell_ << " multicell blocks with max size " + << max_size_multicell_ << " cells in upto " + << max_iter_multicell_ << " iterations." << std::endl; + std::cout << "Average solves per cell (for all cells) was " + << double(num_singlesolves_)/double(grid_.number_of_cells) << std::endl; + } + } + + + + + /// Solve for time-of-flight and a number of tracers. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. + /// \param[in] tracerheads Table containing one row per tracer, and each + /// row contains the source cells for that tracer. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + /// \param[out] tracer_coeff Array of tracer solution coefficients. N*K per cell, + /// where N is equal to tracerheads.size(). All K coefs + /// for a tracer are consecutive, and all tracers' coefs + /// for a cell come before those for the next cell. + void TofDiscGalReorder::solveTofTracer(const double* darcyflux, + const double* porevolume, + const double* source, + const SparseTable& tracerheads, + std::vector& tof_coeff, + std::vector& tracer_coeff) + { + darcyflux_ = darcyflux; + porevolume_ = porevolume; + source_ = source; +#ifndef NDEBUG + // Sanity check for sources. + const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { + // THROW("Sources do not sum to zero: " << cum_src); + MESSAGE("Warning: sources do not sum to zero: " << cum_src); + } +#endif + const int num_basis = basis_func_->numBasisFunc(); + num_tracers_ = tracerheads.size(); + tof_coeff.resize(num_basis*grid_.number_of_cells); + std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0); + tof_coeff_ = &tof_coeff[0]; + rhs_.resize(num_basis*(num_tracers_ + 1)); + jac_.resize(num_basis*num_basis); + orig_jac_.resize(num_basis*num_basis); + basis_.resize(num_basis); + basis_nb_.resize(num_basis); + grad_basis_.resize(num_basis*grid_.dimensions); + velocity_interpolation_->setupFluxes(darcyflux); + + // Set up tracer + tracer_coeff.resize(grid_.number_of_cells*num_tracers_*num_basis); + std::fill(tracer_coeff.begin(), tracer_coeff.end(), 0.0); + if (num_tracers_ > 0) { + tracerhead_by_cell_.clear(); + tracerhead_by_cell_.resize(grid_.number_of_cells, NoTracerHead); + } + for (int tr = 0; tr < num_tracers_; ++tr) { + for (int i = 0; i < tracerheads[tr].size(); ++i) { + const int cell = tracerheads[tr][i]; + basis_func_->addConstant(1.0, &tracer_coeff[cell*num_tracers_*num_basis + tr*num_basis]); + tracer_coeff[cell*num_tracers_ + tr] = 1.0; + tracerhead_by_cell_[cell] = tr; + } + } + + tracer_coeff_ = &tracer_coeff[0]; num_multicell_ = 0; max_size_multicell_ = 0; max_iter_multicell_ = 0; @@ -165,6 +263,13 @@ namespace Opm // This is linear in c_i, so we do not need any nonlinear iterations. // We assemble the jacobian and the right-hand side. The residual is // equal to Res = Jac*c - rhs, and we compute rhs directly. + // + // For tracers, the equation is the same, except for the last + // term being zero (the one with \phi). + // + // The rhs_ vector contains a (Fortran ordering) matrix of all + // right-hand-sides, first for tof and then (optionally) for + // all tracers. const int dim = grid_.dimensions; const int num_basis = basis_func_->numBasisFunc(); @@ -183,6 +288,7 @@ namespace Opm basis_func_->eval(cell, &coord_[0], &basis_[0]); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { + // Only adding to the tof rhs. rhs_[j] += w * basis_[j] * porevolume_[cell] / grid_.cell_volumes[cell]; } } @@ -206,6 +312,8 @@ namespace Opm } if (upstream_cell < 0) { // This is an outer boundary. Assumed tof = 0 on inflow, so no contribution. + // For tracers, a cell with inflow should be marked as a tracer head cell, + // and not be modified. continue; } // Do quadrature over the face to compute @@ -222,12 +330,23 @@ namespace Opm quad.quadPtCoord(quad_pt, &coord_[0]); basis_func_->eval(cell, &coord_[0], &basis_[0]); basis_func_->eval(upstream_cell, &coord_[0], &basis_nb_[0]); + const double w = quad.quadPtWeight(quad_pt); + // Modify tof rhs const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(), tof_coeff_ + num_basis*upstream_cell, 0.0); - const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j]; } + // Modify tracer rhs + if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) { + for (int tr = 0; tr < num_tracers_; ++tr) { + const double* up_tr_co = tracer_coeff_ + num_tracers_*num_basis*upstream_cell + num_basis*tr; + const double tracer_up = std::inner_product(basis_nb_.begin(), basis_nb_.end(), up_tr_co, 0.0); + for (int j = 0; j < num_basis; ++j) { + rhs_[num_basis*(tr + 1) + j] -= w * tracer_up * normal_velocity * basis_[j]; + } + } + } } } @@ -318,7 +437,13 @@ namespace Opm // Solve linear equation. MAT_SIZE_T n = num_basis; - MAT_SIZE_T nrhs = 1; + int num_tracer_to_compute = num_tracers_; + if (num_tracers_) { + if (tracerhead_by_cell_[cell] != NoTracerHead) { + num_tracer_to_compute = 0; + } + } + MAT_SIZE_T nrhs = 1 + num_tracer_to_compute; MAT_SIZE_T lda = num_basis; std::vector piv(num_basis); MAT_SIZE_T ldb = num_basis; @@ -344,7 +469,10 @@ namespace Opm } // The solution ends up in rhs_, so we must copy it. - std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); + std::copy(rhs_.begin(), rhs_.begin() + num_basis, tof_coeff_ + num_basis*cell); + if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) { + std::copy(rhs_.begin() + num_basis, rhs_.end(), tracer_coeff_ + num_tracers_*num_basis*cell); + } // Apply limiter. if (basis_func_->degree() > 0 && use_limiter_ && limiter_usage_ == DuringComputations) { @@ -372,6 +500,7 @@ namespace Opm std::cout << std::endl; #endif applyLimiter(cell, tof_coeff_); + // We do not (yet) apply a limiter to the tracer solution. } } diff --git a/opm/core/tof/TofDiscGalReorder.hpp b/opm/core/tof/TofDiscGalReorder.hpp index b92d2cce..c4e67f2d 100644 --- a/opm/core/tof/TofDiscGalReorder.hpp +++ b/opm/core/tof/TofDiscGalReorder.hpp @@ -35,6 +35,7 @@ namespace Opm class VelocityInterpolationInterface; class DGBasisInterface; namespace parameter { class ParameterGroup; } + template class SparseTable; /// Implements a discontinuous Galerkin solver for /// (single-phase) time-of-flight using reordering. @@ -83,7 +84,7 @@ namespace Opm /// \param[out] tof_coeff Array of time-of-flight solution coefficients. /// The values are ordered by cell, meaning that /// the K coefficients corresponding to the first - /// cell comes before the K coefficients corresponding + /// cell come before the K coefficients corresponding /// to the second cell etc. /// K depends on degree and grid dimension. void solveTof(const double* darcyflux, @@ -91,6 +92,31 @@ namespace Opm const double* source, std::vector& tof_coeff); + /// Solve for time-of-flight and a number of tracers. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. + /// \param[in] tracerheads Table containing one row per tracer, and each + /// row contains the source cells for that tracer. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + /// \param[out] tracer_coeff Array of tracer solution coefficients. N*K per cell, + /// where N is equal to tracerheads.size(). All K coefs + /// for a tracer are consecutive, and all tracers' coefs + /// for a cell come before those for the next cell. + void solveTofTracer(const double* darcyflux, + const double* porevolume, + const double* source, + const SparseTable& tracerheads, + std::vector& tof_coeff, + std::vector& tracer_coeff); + private: virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); @@ -115,11 +141,16 @@ namespace Opm const double* source_; // one volumetric source term per cell boost::shared_ptr basis_func_; double* tof_coeff_; - std::vector rhs_; // single-cell right-hand-side + // For tracers. + double* tracer_coeff_; + int num_tracers_; + enum { NoTracerHead = -1 }; + std::vector tracerhead_by_cell_; + // Used by solveSingleCell(). + std::vector rhs_; // single-cell right-hand-sides std::vector jac_; // single-cell jacobian - std::vector orig_rhs_; // single-cell right-hand-side (copy) + std::vector orig_rhs_; // single-cell right-hand-sides (copy) std::vector orig_jac_; // single-cell jacobian (copy) - // Below: storage for quantities needed by solveSingleCell(). std::vector coord_; mutable std::vector basis_; mutable std::vector basis_nb_; From da5ea0ccae92a33bfb625abf886f111ee61dd297 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 24 Apr 2013 10:40:56 +0200 Subject: [PATCH 14/21] Add call to DG tracer computations if user requests. Used to throw with an 'unimplemented' message. --- examples/compute_tof.cpp | 15 +++++++++------ examples/compute_tof_from_files.cpp | 9 +++++---- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 1d71628a..7f5adc36 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -197,9 +197,6 @@ main(int argc, char** argv) use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } bool compute_tracer = param.getDefault("compute_tracer", false); - if (use_dg && compute_tracer) { - THROW("DG for tracer not yet implemented."); - } // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -254,13 +251,19 @@ main(int argc, char** argv) transport_timer.start(); std::vector tof; std::vector tracer; + Opm::SparseTable tracerheads; + if (compute_tracer) { + buildTracerheadsFromWells(wells->c_wells(), well_state, tracerheads); + } if (use_dg) { - dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); + if (compute_tracer) { + dg_solver->solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tracerheads, tof, tracer); + } else { + dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); + } } else { Opm::TofReorder tofsolver(*grid->c_grid(), use_multidim_upwind); if (compute_tracer) { - Opm::SparseTable tracerheads; - buildTracerheadsFromWells(wells->c_wells(), well_state, tracerheads); tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tracerheads, tof, tracer); } else { tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); diff --git a/examples/compute_tof_from_files.cpp b/examples/compute_tof_from_files.cpp index 723df378..a6580c2f 100644 --- a/examples/compute_tof_from_files.cpp +++ b/examples/compute_tof_from_files.cpp @@ -148,9 +148,6 @@ main(int argc, char** argv) } else { use_multidim_upwind = param.getDefault("use_multidim_upwind", false); } - if (use_dg && compute_tracer) { - THROW("DG for tracer not yet implemented."); - } // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -178,7 +175,11 @@ main(int argc, char** argv) std::vector tof; std::vector tracer; if (use_dg) { - dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof); + if (compute_tracer) { + dg_solver->solveTofTracer(&flux[0], &porevol[0], &src[0], tracerheads, tof, tracer); + } else { + dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof); + } } else { Opm::TofReorder tofsolver(grid, use_multidim_upwind); if (compute_tracer) { From 61829dd238dda744b89c9b5618bf9ac37270666f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 24 Apr 2013 11:27:04 +0200 Subject: [PATCH 15/21] Added limiter for tracer. --- opm/core/tof/TofDiscGalReorder.cpp | 48 +++++++++++++++++++++++++++++- opm/core/tof/TofDiscGalReorder.hpp | 4 +++ 2 files changed, 51 insertions(+), 1 deletion(-) diff --git a/opm/core/tof/TofDiscGalReorder.cpp b/opm/core/tof/TofDiscGalReorder.cpp index e2e4839a..95ef468b 100644 --- a/opm/core/tof/TofDiscGalReorder.cpp +++ b/opm/core/tof/TofDiscGalReorder.cpp @@ -500,7 +500,11 @@ namespace Opm std::cout << std::endl; #endif applyLimiter(cell, tof_coeff_); - // We do not (yet) apply a limiter to the tracer solution. + if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) { + for (int tr = 0; tr < num_tracers_; ++tr) { + applyTracerLimiter(cell, tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis); + } + } } } @@ -721,5 +725,47 @@ namespace Opm + void TofDiscGalReorder::applyTracerLimiter(const int cell, double* local_coeff) + { + // Evaluate the solution in all corners of all faces. Extract max and min. + const int dim = grid_.dimensions; + const int num_basis = basis_func_->numBasisFunc(); + double min_cornerval = 1e100; + double max_cornerval = -1e100; + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) { + const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode]; + basis_func_->eval(cell, nc, &basis_[0]); + const double tracer_corner = std::inner_product(basis_.begin(), basis_.end(), + local_coeff, 0.0); + min_cornerval = std::min(min_cornerval, tracer_corner); + max_cornerval = std::max(min_cornerval, tracer_corner); + } + } + const double average = basis_func_->functionAverage(local_coeff); + if (average < 0.0 || average > 1.0) { + // Adjust average. Flatten gradient. + std::fill(local_coeff, local_coeff + num_basis, 0.0); + if (average > 1.0) { + basis_func_->addConstant(1.0, local_coeff); + } + } else { + // Possibly adjust gradient. + double factor = 1.0; + if (min_cornerval < 0.0) { + factor = average/(average - min_cornerval); + } + if (max_cornerval > 1.0) { + factor = std::min(factor, (1.0 - average)/(max_cornerval - average)); + } + if (factor != 1.0) { + basis_func_->multiplyGradient(factor, local_coeff); + } + } + } + + + } // namespace Opm diff --git a/opm/core/tof/TofDiscGalReorder.hpp b/opm/core/tof/TofDiscGalReorder.hpp index c4e67f2d..b13affac 100644 --- a/opm/core/tof/TofDiscGalReorder.hpp +++ b/opm/core/tof/TofDiscGalReorder.hpp @@ -174,6 +174,10 @@ namespace Opm void applyLimiterAsSimultaneousPostProcess(); double totalFlux(const int cell) const; double minCornerVal(const int cell, const int face) const; + + // Apply a simple (restrict to [0,1]) limiter. + // Intended for tracers. + void applyTracerLimiter(const int cell, double* local_coeff); }; } // namespace Opm From a6e3bbd6109826988c4fd24973f7c91d46c1e491 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Thu, 2 May 2013 19:04:09 +0200 Subject: [PATCH 16/21] Updated calls to ecl_util_fmt_file() - result is by reference --- opm/core/io/eclipse/EclipseGridParser.cpp | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/opm/core/io/eclipse/EclipseGridParser.cpp b/opm/core/io/eclipse/EclipseGridParser.cpp index 8a1905ec..ed4b246c 100644 --- a/opm/core/io/eclipse/EclipseGridParser.cpp +++ b/opm/core/io/eclipse/EclipseGridParser.cpp @@ -966,9 +966,16 @@ ecl_grid_type * EclipseGridParser::newGrid( ) { void EclipseGridParser::saveEGRID( const std::string & filename) const { bool endian_flip = true;//ECL_ENDIAN_FLIP; - bool fmt_file = ecl_util_fmt_file( filename.c_str() ); + bool fmt_file; struct grdecl grdecl = get_grdecl(); - fortio_type * fortio = fortio_open_writer( filename.c_str() , fmt_file , endian_flip ); + fortio_type * fortio; + + if (!ecl_util_fmt_file( filename.c_str() , &fmt_file)) { + cerr << "Could not determine formatted/unformatted status of file:" << filename << " non-standard name?" << endl; + throw exception(); + } + + fortio = fortio_open_writer( filename.c_str() , fmt_file , endian_flip ); { float * mapaxes = NULL; if (grdecl.mapaxes != NULL) { @@ -1023,9 +1030,16 @@ void EclipseGridParser::save_kw( fortio_type * fortio , const std::string & kw , void EclipseGridParser::saveINIT( const std::string & filename , const ecl_grid_type * ecl_grid) { int phases = ECL_OIL_PHASE + ECL_WATER_PHASE; - bool fmt_file = ecl_util_fmt_file( filename.c_str() ); bool endian_flip = true;//ECL_ENDIAN_FLIP; - fortio_type * fortio = fortio_open_writer( filename.c_str() , fmt_file , endian_flip ); + bool fmt_file; + fortio_type * fortio; + + if (!ecl_util_fmt_file( filename.c_str() , &fmt_file)) { + cerr << "Could not determine formatted/unformatted status of file:" << filename << " non-standard name?" << endl; + throw exception(); + } + fortio = fortio_open_writer( filename.c_str() , fmt_file , endian_flip ); + { ecl_kw_type * poro_kw = newEclKW( PORO_KW , ECL_FLOAT_TYPE ); time_t start_date; From 1281d17b60767736271c15f4f032d58994f1363e Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Thu, 2 May 2013 19:06:43 +0200 Subject: [PATCH 17/21] Changed writeECLDate() to write UNFORMATTED restart files. --- opm/core/io/eclipse/writeECLData.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/io/eclipse/writeECLData.cpp b/opm/core/io/eclipse/writeECLData.cpp index 1729fdb6..ce2b258c 100644 --- a/opm/core/io/eclipse/writeECLData.cpp +++ b/opm/core/io/eclipse/writeECLData.cpp @@ -94,7 +94,7 @@ namespace Opm const std::string& base_name) { ecl_file_enum file_type = ECL_UNIFIED_RESTART_FILE; // Alternatively ECL_RESTART_FILE for multiple restart files. - bool fmt_file = true; + bool fmt_file = false; char * filename = ecl_util_alloc_filename(output_dir.c_str() , base_name.c_str() , file_type , fmt_file , current_step ); int phases = ECL_OIL_PHASE + ECL_WATER_PHASE; From 0fb4f14b28a50dd21cae045bb5ea24bb879f53b6 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Mon, 13 May 2013 10:50:22 +0200 Subject: [PATCH 18/21] Provide DTD schema to describe parameter file format Notice that this definition specifies a superset of the actual format as the ParameterGroup tag is context-dependent (which cannot be captured by DTD schemata). --- opm/core/utility/parameters/param.dtd | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 opm/core/utility/parameters/param.dtd diff --git a/opm/core/utility/parameters/param.dtd b/opm/core/utility/parameters/param.dtd new file mode 100644 index 00000000..37fcf8b6 --- /dev/null +++ b/opm/core/utility/parameters/param.dtd @@ -0,0 +1,26 @@ + + + + + + + + + + + + + + + + From 229cf6d134deb0d069884f179a0c843d84f3d87e Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Thu, 18 Apr 2013 13:35:25 +0200 Subject: [PATCH 19/21] Probe for library without using STL If dune-grid has been compiled with _GLIBCXX_DEBUG, linking the test program will fail since probing doesn't get the flags that are specified for the main project, and thus configuration halts. Instead, we use a smaller test program which doesn't need proper STL linkage to compile and link. --- cmake/Modules/Finddune-grid.cmake | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/cmake/Modules/Finddune-grid.cmake b/cmake/Modules/Finddune-grid.cmake index 1ee8a062..ae648fea 100644 --- a/cmake/Modules/Finddune-grid.cmake +++ b/cmake/Modules/Finddune-grid.cmake @@ -34,12 +34,9 @@ find_opm_package ( # test program "#include -#include int main (void) { - std::vector coords; - Dune::OneDGrid grid(coords); + Dune::OneDGrid grid(1, 0., 1.); return grid.lbegin<0>(0) == grid.lend<0>(0); - return 0; } " # config variables From 238bd6b8919a8d4b1c30e56b2870c50c84fc30fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 13 May 2013 16:17:41 +0200 Subject: [PATCH 20/21] Ensure (average of) tracers sum to one for DG1. --- opm/core/tof/TofDiscGalReorder.cpp | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/opm/core/tof/TofDiscGalReorder.cpp b/opm/core/tof/TofDiscGalReorder.cpp index 95ef468b..a5578f7e 100644 --- a/opm/core/tof/TofDiscGalReorder.cpp +++ b/opm/core/tof/TofDiscGalReorder.cpp @@ -506,6 +506,26 @@ namespace Opm } } } + + // Ensure that tracer averages sum to 1. + if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) { + std::vector tr_aver(num_tracers_); + double tr_sum = 0.0; + for (int tr = 0; tr < num_tracers_; ++tr) { + const double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis; + tr_aver[tr] = basis_func_->functionAverage(local_basis); + tr_sum += tr_aver[tr]; + } + if (tr_sum == 0.0) { + std::cout << "Tracer sum is zero in cell " << cell << std::endl; + } else { + for (int tr = 0; tr < num_tracers_; ++tr) { + const double increment = tr_aver[tr]/tr_sum - tr_aver[tr]; + double* local_basis = tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis; + basis_func_->addConstant(increment, local_basis); + } + } + } } From 53b082af282f5018950ed9a0bf721309bc400ca0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 14 May 2013 10:19:12 +0200 Subject: [PATCH 21/21] Always interpret 'ref_pressure' as a double Otherwise, when specifying (e.g.) ref_pressure=1.0e-5 (1 Pascal in bars), the value gets reinterpreted as ref_pressure=1 which is one bar. --- opm/core/simulator/initState_impl.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index d8e53bd6..b6492bb9 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -337,7 +337,7 @@ namespace Opm } } state.setFirstSat(left_cells, props, State::MaxSat); - const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa; + const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; std::fill(state.pressure().begin(), state.pressure().end(), init_p); } else if (segregation_testcase) { // Warn against error-prone usage. @@ -351,7 +351,7 @@ namespace Opm const double woc = param.get("water_oil_contact"); initWaterOilContact(grid, props, woc, WaterAbove, state); // Initialise pressure to hydrostatic state. - const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa; + const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; 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")) { @@ -366,7 +366,7 @@ namespace Opm const double woc = param.get("water_oil_contact"); initWaterOilContact(grid, props, woc, WaterBelow, state); // Initialise pressure to hydrostatic state. - const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa; + const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; initHydrostaticPressure(grid, props.density(), woc, gravity, woc, ref_p, state); } else if (param.has("init_saturation")) { // Initialise water saturation to init_saturation parameter. @@ -376,7 +376,7 @@ namespace Opm state.saturation()[2*cell + 1] = 1.0 - init_saturation; } // Initialise pressure to hydrostatic state. - const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa; + const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; const double rho = props.density()[0]*init_saturation + props.density()[1]*(1.0 - init_saturation); const double dens[2] = { rho, rho }; const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1]; @@ -384,7 +384,7 @@ namespace Opm } else { // 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 ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; const double rho = props.density()[1]; const double dens[2] = { rho, rho }; const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1]; @@ -432,7 +432,7 @@ namespace Opm } } state.setFirstSat(left_cells, props, State::MaxSat); - const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa; + const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; std::fill(state.pressure().begin(), state.pressure().end(), init_p); } else if (param.has("water_oil_contact")) { // Warn against error-prone usage. @@ -446,12 +446,12 @@ namespace Opm const double woc = param.get("water_oil_contact"); initWaterOilContact(grid, props, woc, WaterBelow, state); // Initialise pressure to hydrostatic state. - const double ref_p = param.getDefault("ref_pressure", 100)*unit::barsa; + const double ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; initHydrostaticPressure(grid, props, woc, gravity, woc, ref_p, state); } else { // 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 ref_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; const double ref_z = grid.cell_centroids[0 + grid.dimensions - 1]; const double woc = -1e100; initHydrostaticPressure(grid, props, woc, gravity, ref_z, ref_p, state);