/* 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 . */ #include "config.h" #include #include #include #include #include #include #include #include namespace Opm { /// Construct solver. /// \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) : grid_(grid), darcyflux_(0), porevolume_(0), source_(0), tof_(0), tracer_(0), num_tracers_(0), gauss_seidel_tol_(1e-3), use_multidim_upwind_(use_multidim_upwind) { } /// Solve for time-of-flight. /// \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[out] tof Array of time-of-flight values. void TofReorder::solveTof(const double* darcyflux, const double* porevolume, const double* source, std::vector& tof) { 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); } #endif 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); } num_tracers_ = 0; 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; } } /// 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 Array of time-of-flight values (1 per cell). /// \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, const SparseTable& tracerheads, std::vector& tof, std::vector& tracer) { 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); } #endif tof.resize(grid_.number_of_cells); std::fill(tof.begin(), tof.end(), 0.0); tof_ = &tof[0]; // Find the tracer heads (injectors). num_tracers_ = tracerheads.size(); tracer.resize(grid_.number_of_cells*num_tracers_); std::fill(tracer.begin(), tracer.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]; tracer[cell*num_tracers_ + tr] = 1.0; tracerhead_by_cell_[cell] = tr; } } 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; } } void TofReorder::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 // 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) { 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) { upwind_term += flux*tof_[other]; 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 { downwind_flux += flux; } } // Compute tof. tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux; // Compute tracers (if any). // Do not change tracer solution in source cells. 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; } } } void TofReorder::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; // Compute cell flux if (cell == grid_.face_cells[2*f]) { flux = darcyflux_[f]; } else { flux =-darcyflux_[f]; } // Add flux to upwind_term or downwind_term_[face|cell_factor]. if (flux < 0.0) { upwind_term += flux*face_tof_[f]; } else if (flux > 0.0) { 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 TofReorder::solveMultiCell(const int num_cells, const int* cells) { ++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. 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 = tof_[cell]; solveSingleCell(cell); max_delta = std::max(max_delta, std::fabs(tof_[cell] - tof_before)); } // std::cout << "Max delta = " << max_delta << std::endl; } max_iter_multicell_ = std::max(max_iter_multicell_, num_iter); } // 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 TofReorder::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" // 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. // 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]; 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); } // 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); } } } // Indentify adjacent faces with inflows, compute omega_star, omega, // add up contributions. const int num_adj = adj_faces_.size(); // 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; 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; face_term += omega*face_tof_[f]; cell_term_factor += (1.0 - omega); } face_term /= double(num_adj); cell_term_factor /= double(num_adj); } } // namespace Opm