diff --git a/opm/core/tof/TofDiscGalReorder.cpp b/opm/core/tof/TofDiscGalReorder.cpp index 2c2331eb..5bfa5294 100644 --- a/opm/core/tof/TofDiscGalReorder.cpp +++ b/opm/core/tof/TofDiscGalReorder.cpp @@ -59,6 +59,8 @@ namespace Opm basis_func_.reset(new DGBasisBoundedTotalDegree(grid_, dg_degree)); } + tracers_ensure_unity_ = param.getDefault("tracers_ensure_unity", true); + use_cvi_ = param.getDefault("use_cvi", use_cvi_); use_limiter_ = param.getDefault("use_limiter", use_limiter_); if (use_limiter_) { @@ -273,13 +275,66 @@ namespace Opm // right-hand-sides, first for tof and then (optionally) for // all tracers. - 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); + // Add cell contributions to res_ and jac_. + cellContribs(cell); + + // Add face contributions to res_ and jac_. + faceContribs(cell); + + // Solve linear equation. + solveLinearSystem(cell); + + // The solution ends up in rhs_, so we must copy it. + 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) { + applyLimiter(cell, tof_coeff_); + 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); + } + } + } + + // Ensure that tracer averages sum to 1. + if (num_tracers_ && tracers_ensure_unity_ && 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); + } + } + } + } + + + + + void TofDiscGalReorder::cellContribs(const int cell) + { + const int num_basis = basis_func_->numBasisFunc(); + const int dim = grid_.dimensions; + // Compute cell residual contribution. { const int deg_needed = basis_func_->degree(); @@ -296,7 +351,70 @@ namespace Opm } } - // Compute upstream residual contribution. + // Compute cell jacobian contribution. We use Fortran ordering + // for jac_, i.e. rows cycling fastest. + { + // Even with ECVI velocity interpolation, degree of precision 1 + // is sufficient for optimal convergence order for DG1 when we + // use linear (total degree 1) basis functions. + // With bi(tri)-linear basis functions, it still seems sufficient + // for convergence order 2, but the solution looks much better and + // has significantly lower error with degree of precision 2. + // For now, we err on the side of caution, and use 2*degree, even + // though this is wasteful for the pure linear basis functions. + // const int deg_needed = 2*basis_func_->degree() - 1; + const int deg_needed = 2*basis_func_->degree(); + CellQuadrature quad(grid_, cell, deg_needed); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // b_i (v \cdot \grad b_j) + quad.quadPtCoord(quad_pt, &coord_[0]); + basis_func_->eval(cell, &coord_[0], &basis_[0]); + basis_func_->evalGrad(cell, &coord_[0], &grad_basis_[0]); + velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + for (int dd = 0; dd < dim; ++dd) { + jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd]; + } + } + } + } + } + + // Compute downstream jacobian contribution from sink terms. + // Contribution from inflow sources would be + // similar to the contribution from upstream faces, but + // it is zero since we let all external inflow be associated + // with a zero tof. + if (source_[cell] < 0.0) { + // A sink. + const double flux = -source_[cell]; // Sign convention for flux: outflux > 0. + const double flux_density = flux / grid_.cell_volumes[cell]; + // Do quadrature over the cell to compute + // \int_{K} b_i flux b_j dx + CellQuadrature quad(grid_, cell, 2*basis_func_->degree()); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + quad.quadPtCoord(quad_pt, &coord_[0]); + basis_func_->eval(cell, &coord_[0], &basis_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j]; + } + } + } + } + } + + + + + void TofDiscGalReorder::faceContribs(const int cell) + { + const int num_basis = basis_func_->numBasisFunc(); + + // Compute upstream residual contribution from faces. 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; @@ -352,37 +470,6 @@ namespace Opm } } - // Compute cell jacobian contribution. We use Fortran ordering - // for jac_, i.e. rows cycling fastest. - { - // Even with ECVI velocity interpolation, degree of precision 1 - // is sufficient for optimal convergence order for DG1 when we - // use linear (total degree 1) basis functions. - // With bi(tri)-linear basis functions, it still seems sufficient - // for convergence order 2, but the solution looks much better and - // has significantly lower error with degree of precision 2. - // For now, we err on the side of caution, and use 2*degree, even - // though this is wasteful for the pure linear basis functions. - // const int deg_needed = 2*basis_func_->degree() - 1; - const int deg_needed = 2*basis_func_->degree(); - CellQuadrature quad(grid_, cell, deg_needed); - for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { - // b_i (v \cdot \grad b_j) - quad.quadPtCoord(quad_pt, &coord_[0]); - basis_func_->eval(cell, &coord_[0], &basis_[0]); - basis_func_->evalGrad(cell, &coord_[0], &grad_basis_[0]); - velocity_interpolation_->interpolate(cell, &coord_[0], &velocity_[0]); - const double w = quad.quadPtWeight(quad_pt); - for (int j = 0; j < num_basis; ++j) { - for (int i = 0; i < num_basis; ++i) { - for (int dd = 0; dd < dim; ++dd) { - jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd]; - } - } - } - } - } - // Compute downstream jacobian contribution from faces. for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { const int face = grid_.cell_faces[hface]; @@ -412,33 +499,17 @@ namespace Opm } } } + } - // Compute downstream jacobian contribution from sink terms. - // Contribution from inflow sources would be - // similar to the contribution from upstream faces, but - // it is zero since we let all external inflow be associated - // with a zero tof. - if (source_[cell] < 0.0) { - // A sink. - const double flux = -source_[cell]; // Sign convention for flux: outflux > 0. - const double flux_density = flux / grid_.cell_volumes[cell]; - // Do quadrature over the cell to compute - // \int_{K} b_i flux b_j dx - CellQuadrature quad(grid_, cell, 2*basis_func_->degree()); - for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { - quad.quadPtCoord(quad_pt, &coord_[0]); - basis_func_->eval(cell, &coord_[0], &basis_[0]); - const double w = quad.quadPtWeight(quad_pt); - for (int j = 0; j < num_basis; ++j) { - for (int i = 0; i < num_basis; ++i) { - jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j]; - } - } - } - } - // Solve linear equation. - MAT_SIZE_T n = num_basis; + + // This function assumes that jac_ and rhs_ contain the + // linear system to be solved. They are stored in orig_jac_ + // and orig_rhs_, then the system is solved via LAPACK, + // overwriting the input data (jac_ and rhs_). + void TofDiscGalReorder::solveLinearSystem(const int cell) + { + MAT_SIZE_T n = basis_func_->numBasisFunc(); int num_tracer_to_compute = num_tracers_; if (num_tracers_) { if (tracerhead_by_cell_[cell] != NoTracerHead) { @@ -446,9 +517,9 @@ namespace Opm } } 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; + MAT_SIZE_T lda = n; + std::vector piv(n); + MAT_SIZE_T ldb = n; MAT_SIZE_T info = 0; orig_jac_ = jac_; orig_rhs_ = rhs_; @@ -469,65 +540,6 @@ namespace Opm } OPM_THROW(std::runtime_error, "Lapack error: " << info << " encountered in cell " << cell); } - - // The solution ends up in rhs_, so we must copy it. - 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) { -#ifdef EXTRA_VERBOSE - std::cout << "Cell: " << cell << " "; - std::cout << "v = "; - for (int dd = 0; dd < dim; ++dd) { - std::cout << velocity_[dd] << ' '; - } - std::cout << " grad tau = "; - for (int dd = 0; dd < dim; ++dd) { - std::cout << tof_coeff_[num_basis*cell + dd + 1] << ' '; - } - const double prod = std::inner_product(velocity_.begin(), velocity_.end(), - tof_coeff_ + num_basis*cell + 1, 0.0); - const double vv = std::inner_product(velocity_.begin(), velocity_.end(), - velocity_.begin(), 0.0); - const double gg = std::inner_product(tof_coeff_ + num_basis*cell + 1, - tof_coeff_ + num_basis*cell + num_basis, - tof_coeff_ + num_basis*cell + 1, 0.0); - std::cout << " prod = " << std::inner_product(velocity_.begin(), velocity_.end(), - tof_coeff_ + num_basis*cell + 1, 0.0); - std::cout << " normalized = " << prod/std::sqrt(vv*gg); - std::cout << " angle = " << std::acos(prod/std::sqrt(vv*gg))*360.0/(2.0*M_PI); - std::cout << std::endl; -#endif - applyLimiter(cell, tof_coeff_); - 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); - } - } - } - - // 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); - } - } - } } diff --git a/opm/core/tof/TofDiscGalReorder.hpp b/opm/core/tof/TofDiscGalReorder.hpp index 2e8fa317..880e864d 100644 --- a/opm/core/tof/TofDiscGalReorder.hpp +++ b/opm/core/tof/TofDiscGalReorder.hpp @@ -121,6 +121,10 @@ namespace Opm virtual void solveSingleCell(const int cell); virtual void solveMultiCell(const int num_cells, const int* cells); + void cellContribs(const int cell); + void faceContribs(const int cell); + void solveLinearSystem(const int cell); + private: // Disable copying and assignment. TofDiscGalReorder(const TofDiscGalReorder&); @@ -146,6 +150,7 @@ namespace Opm int num_tracers_; enum { NoTracerHead = -1 }; std::vector tracerhead_by_cell_; + bool tracers_ensure_unity_; // Used by solveSingleCell(). std::vector rhs_; // single-cell right-hand-sides std::vector jac_; // single-cell jacobian diff --git a/opm/core/tof/TofReorder.cpp b/opm/core/tof/TofReorder.cpp index 04a169e5..54953e5b 100644 --- a/opm/core/tof/TofReorder.cpp +++ b/opm/core/tof/TofReorder.cpp @@ -42,8 +42,6 @@ namespace Opm porevolume_(0), source_(0), tof_(0), - tracer_(0), - num_tracers_(0), gauss_seidel_tol_(1e-3), use_multidim_upwind_(use_multidim_upwind) { @@ -71,7 +69,8 @@ namespace Opm // 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) { - OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src); + // OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src); + OPM_MESSAGE("Warning: sources do not sum to zero: " << cum_src); } #endif tof.resize(grid_.number_of_cells); @@ -80,17 +79,11 @@ namespace Opm if (use_multidim_upwind_) { face_tof_.resize(grid_.number_of_faces); std::fill(face_tof_.begin(), face_tof_.end(), 0.0); + face_part_tof_.resize(grid_.face_nodepos[grid_.number_of_faces]); + std::fill(face_part_tof_.begin(), face_part_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; - } + compute_tracer_ = false; + executeSolve(); } @@ -117,43 +110,72 @@ namespace Opm darcyflux_ = darcyflux; porevolume_ = porevolume; source_ = source; + const int num_cells = grid_.number_of_cells; #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) { + const double cum_src = std::accumulate(source, source + num_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + num_cells)*1e-2) { OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src); } #endif - tof.resize(grid_.number_of_cells); + tof.resize(num_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); + if (use_multidim_upwind_) { + face_tof_.resize(grid_.number_of_faces); + std::fill(face_tof_.begin(), face_tof_.end(), 0.0); + face_part_tof_.resize(grid_.face_nodepos[grid_.number_of_faces]); + std::fill(face_part_tof_.begin(), face_part_tof_.end(), 0.0); } - for (int tr = 0; tr < num_tracers_; ++tr) { + + // Execute solve for tof + compute_tracer_ = false; + executeSolve(); + + // Find the tracer heads (injectors). + const int num_tracers = tracerheads.size(); + tracer.resize(num_cells*num_tracers); + std::fill(tracer.begin(), tracer.end(), 0.0); + if (num_tracers > 0) { + tracerhead_by_cell_.clear(); + tracerhead_by_cell_.resize(num_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; + tracer[num_cells * tr + cell] = 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); - OPM_THROW(std::runtime_error, "Multidimensional upwind not yet implemented for tracer."); + // Execute solve for tracers. + std::vector fake_pv(num_cells, 0.0); + porevolume_ = fake_pv.data(); + for (int tr = 0; tr < num_tracers; ++tr) { + tof_ = tracer.data() + tr * num_cells; + compute_tracer_ = true; + executeSolve(); } + + // Write output tracer data (transposing the computed data). + std::vector computed = tracer; + for (int cell = 0; cell < num_cells; ++cell) { + for (int tr = 0; tr < num_tracers; ++tr) { + tracer[num_tracers * cell + tr] = computed[num_cells * tr + cell]; + } + } + } + + + + + void TofReorder::executeSolve() + { num_multicell_ = 0; max_size_multicell_ = 0; max_iter_multicell_ = 0; - reorderAndTransport(grid_, darcyflux); + reorderAndTransport(grid_, darcyflux_); if (num_multicell_ > 0) { std::cout << num_multicell_ << " multicell blocks with max size " << max_size_multicell_ << " cells in upto " @@ -163,7 +185,6 @@ namespace Opm - void TofReorder::solveSingleCell(const int cell) { if (use_multidim_upwind_) { @@ -176,10 +197,9 @@ 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; - } + if (compute_tracer_ && tracerhead_by_cell_[cell] != NoTracerHead) { + // This is a tracer head cell, already has solution. + return; } double upwind_term = 0.0; double downwind_flux = std::max(-source_[cell], 0.0); @@ -202,11 +222,6 @@ namespace Opm // 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; @@ -215,14 +230,6 @@ namespace Opm // 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; - } - } } @@ -260,7 +267,11 @@ namespace Opm } // Compute tof for cell. - tof_[cell] = (porevolume_[cell] - upwind_term - downwind_term_face)/downwind_term_cell_factor; + if (compute_tracer_ && tracerhead_by_cell_[cell] != NoTracerHead) { + // Do nothing to the value in this cell, since we are at a tracer head. + } else { + 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) { @@ -270,6 +281,20 @@ namespace Opm double fterm, cterm_factor; multidimUpwindTerms(f, cell, fterm, cterm_factor); face_tof_[f] = fterm + cterm_factor*tof_[cell]; + + // Combine locally computed (for each adjacent vertex) terms, with uniform weighting. + const int* face_nodes_beg = grid_.face_nodes + grid_.face_nodepos[f]; + const int* face_nodes_end = grid_.face_nodes + grid_.face_nodepos[f + 1]; + const int num_terms = face_nodes_end - face_nodes_beg; + assert(num_terms == 2 || grid_.dimensions != 2); + for (const int* fn_iter = face_nodes_beg; fn_iter < face_nodes_end; ++fn_iter) { + double loc_face_term = 0.0; + double loc_cell_term_factor = 0.0; + const int node_pos = fn_iter - grid_.face_nodes; + localMultidimUpwindTerms(f, cell, node_pos, + loc_face_term, loc_cell_term_factor); + face_part_tof_[node_pos] = loc_face_term + loc_cell_term_factor * tof_[cell]; + } } } } @@ -303,7 +328,10 @@ namespace Opm - // Assumes that face_tof_[f] is known for all upstream faces f of upwind_cell. + // Assumes that face_part_tof_[node_pos] is known for all inflow + // faces to 'upwind_cell' sharing vertices with 'face'. The index + // 'node_pos' is the same as the one used for the grid face-node + // connectivity. // 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). @@ -314,71 +342,112 @@ namespace Opm double& face_term, double& cell_term_factor) const { - // Implements multidim upwind according to + // Implements multidim upwind inspired by // "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. + // a good idea. We have here attempted some generalization, by treating each face-part + // (association of a face and a vertex) as possibly influencing all downwind face-parts + // of the neighbouring cell that share the same vertex. + // The current implementation aims to reproduce 2d results for extruded 3d grids. - // Identify the adjacent faces of the upwind cell. + // Combine locally computed (for each adjacent vertex) terms, with uniform weighting. 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]); + const int num_terms = face_nodes_end - face_nodes_beg; + assert(num_terms == 2 || grid_.dimensions != 2); 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); + for (const int* fn_iter = face_nodes_beg; fn_iter < face_nodes_end; ++fn_iter) { + double loc_face_term = 0.0; + double loc_cell_term_factor = 0.0; + localMultidimUpwindTerms(face, upwind_cell, fn_iter - grid_.face_nodes, + loc_face_term, loc_cell_term_factor); + face_term += loc_face_term; + cell_term_factor += loc_cell_term_factor; } - face_term /= double(num_adj); - cell_term_factor /= double(num_adj); + face_term /= double(num_terms); + cell_term_factor /= double(num_terms); + } + + namespace { + double weightFunc(const double w) + { + // SPU + // return 0.0; + // TMU + return w > 0.0 ? std::min(w, 1.0) : 0.0; + // SMU + // return w > 0.0 ? w/(1.0 + w) : 0.0; + } + } + + + + + void TofReorder::localMultidimUpwindTerms(const int face, + const int upwind_cell, + const int node_pos, + double& face_term, + double& cell_term_factor) const + { + // Loop over all faces adjacent to the given cell and the + // vertex in position node_pos. + // If that part's influx is positive, we store it, and also its associated + // node position. + std::vector influx; + std::vector node_pos_influx; + influx.reserve(5); + node_pos_influx.reserve(5); + const int node = grid_.face_nodes[node_pos]; + 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) { + // Find out if the face 'f' is adjacent to vertex 'node'. + 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]; + const int* pos = std::find(f_nodes_beg, f_nodes_end, node); + const int node_pos2 = pos - grid_.face_nodes; + const bool is_adj = (pos != f_nodes_end); + if (is_adj) { + const int num_parts = f_nodes_end - f_nodes_beg; + const double influx_sign = (grid_.face_cells[2*f] == upwind_cell) ? -1.0 : 1.0; + const double part_influx = influx_sign * darcyflux_[f] / double(num_parts); + if (part_influx > 0.0) { + influx.push_back(part_influx); + node_pos_influx.push_back(node_pos2); + } + } + } + } + + // Now we may compute the weighting of the upwind terms. + const int num_parts = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + const double outflux_sign = (grid_.face_cells[2*face] == upwind_cell) ? 1.0 : -1.0; + const double part_outflux = outflux_sign * darcyflux_[face] / double(num_parts); + const double sum_influx = std::accumulate(influx.begin(), influx.end(), 0.0); + const double w_factor = weightFunc(sum_influx / part_outflux); + const int num_influx = influx.size(); + std::vector w(num_influx); + face_term = 0.0; + for (int ii = 0; ii < num_influx; ++ii) { + w[ii] = (influx[ii] / sum_influx) * w_factor; + face_term += w[ii] * face_part_tof_[node_pos_influx[ii]]; + } + const double sum_w = std::accumulate(w.begin(), w.end(), 0.0); + cell_term_factor = 1.0 - sum_w; + const double tol = 1e-5; + if (cell_term_factor < -tol && cell_term_factor > 1.0 + tol) { + OPM_THROW(std::logic_error, "cell_term_factor outside [0,1]: " << cell_term_factor); + } + cell_term_factor = std::min(std::max(cell_term_factor, 0.0), 1.0); + assert(cell_term_factor >= 0.0); + assert(cell_term_factor <= 1.0); + } + } // namespace Opm diff --git a/opm/core/tof/TofReorder.hpp b/opm/core/tof/TofReorder.hpp index 4abd640a..8dc88b53 100644 --- a/opm/core/tof/TofReorder.hpp +++ b/opm/core/tof/TofReorder.hpp @@ -80,6 +80,7 @@ namespace Opm std::vector& tracer); private: + void executeSolve(); virtual void solveSingleCell(const int cell); void solveSingleCellMultidimUpwind(const int cell); void assembleSingleCell(const int cell, @@ -90,6 +91,8 @@ namespace Opm void multidimUpwindTerms(const int face, const int upwind_cell, double& face_term, double& cell_term_factor) const; + void localMultidimUpwindTerms(const int face, const int upwind_cell, const int node_pos, + double& face_term, double& cell_term_factor) const; private: const UnstructuredGrid& grid_; @@ -97,8 +100,7 @@ namespace Opm const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell double* tof_; - double* tracer_; - int num_tracers_; + bool compute_tracer_; enum { NoTracerHead = -1 }; std::vector tracerhead_by_cell_; // For solveMultiCell(): @@ -109,7 +111,7 @@ namespace Opm // 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. + std::vector face_part_tof_; // For multidim upwind face tofs. }; } // namespace Opm