From 604be4587133cd240d561621263b4033aefc5717 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 21 Jan 2013 14:55:27 +0100 Subject: [PATCH] Add methods totalFlux() and minCornerVal(). Also started refactoring applyMinUpwindLimiter() using the added methods. --- .../TransportModelTracerTofDiscGal.cpp | 77 ++++++++++++------- .../TransportModelTracerTofDiscGal.hpp | 6 +- 2 files changed, 54 insertions(+), 29 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index fdbec4409..24846cf58 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -444,34 +444,13 @@ namespace Opm // minimum. // 2. The TOF shall not be below zero in any point. - // Find total upstream/downstream fluxes. - double upstream_flux = 0.0; - double downstream_flux = 0.0; - for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { - const int face = grid_.cell_faces[hface]; - double flux = 0.0; - if (cell == grid_.face_cells[2*face]) { - flux = darcyflux_[face]; - } else { - flux = -darcyflux_[face]; - } - if (flux < 0.0) { - upstream_flux += flux; - } else { - downstream_flux += flux; - } - } - // In the presence of sources, significant fluxes may be missing from the computed fluxes, - // setting the total flux to the (positive) maximum avoids this: since source is either - // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. - const double total_flux = std::max(-upstream_flux, downstream_flux); - - // Find minimum tof on upstream faces and for this cell. + // Find minimum tof on upstream faces/cells and for this cell. const int dim = grid_.dimensions; const int num_basis = basis_func_->numBasisFunc(); double min_upstream_tof = 1e100; double min_here_tof = 1e100; int num_upstream_faces = 0; + const double total_flux = totalFlux(cell); for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { const int face = grid_.cell_faces[hface]; double flux = 0.0; @@ -490,12 +469,9 @@ namespace Opm bool interior = (upstream_cell >= 0); // Evaluate the solution in all corners. + min_here_tof = std::min(min_here_tof, minCornerVal(cell, face)); 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 tof_here = std::inner_product(basis_.begin(), basis_.end(), - tof_coeff_ + num_basis*cell, 0.0); - min_here_tof = std::min(min_here_tof, tof_here); if (upstream) { if (interior) { const double* upstream_coef = tof_coeff_ + num_basis*upstream_cell; @@ -583,4 +559,51 @@ namespace Opm + double TransportModelTracerTofDiscGal::totalFlux(const int cell) const + { + // Find total upstream/downstream fluxes. + double upstream_flux = 0.0; + double downstream_flux = 0.0; + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + } else { + flux = -darcyflux_[face]; + } + if (flux < 0.0) { + upstream_flux += flux; + } else { + downstream_flux += flux; + } + } + // In the presence of sources, significant fluxes may be missing from the computed fluxes, + // setting the total flux to the (positive) maximum avoids this: since source is either + // inflow or outflow, not both, either upstream_flux or downstream_flux must be correct. + return std::max(-upstream_flux, downstream_flux); + } + + + + + double TransportModelTracerTofDiscGal::minCornerVal(const int cell, const int face) const + { + // Evaluate the solution in all corners. + const int dim = grid_.dimensions; + const int num_basis = basis_func_->numBasisFunc(); + double min_cornerval = 1e100; + 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 tof_corner = std::inner_product(basis_.begin(), basis_.end(), + tof_coeff_ + num_basis*cell, 0.0); + min_cornerval = std::min(min_cornerval, tof_corner); + } + return min_cornerval; + } + + + + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 9ce8cec22..cf55ca8cb 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -119,8 +119,8 @@ namespace Opm std::vector orig_jac_; // single-cell jacobian (copy) // Below: storage for quantities needed by solveSingleCell(). std::vector coord_; - std::vector basis_; - std::vector basis_nb_; + mutable std::vector basis_; + mutable std::vector basis_nb_; std::vector grad_basis_; std::vector velocity_; @@ -133,6 +133,8 @@ namespace Opm void applyMinUpwindLimiter(const int cell, const bool face_min, double* tof); void applyLimiterAsPostProcess(); void applyLimiterAsSimultaneousPostProcess(); + double totalFlux(const int cell) const; + double minCornerVal(const int cell, const int face) const; }; } // namespace Opm