From 832072a5cbd2f502a8f5736d33de2a051a4758d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 8 Jan 2013 11:04:43 +0100 Subject: [PATCH] Improvements to DG1 flux limiter. Skeleton in place for increased flexibility in methods and usage. (So far behaviour choices are hardcoded, though.) Added relative flux thresholding to existing limiter to avoid flux noise strongly affecting solution. For example no-flow boundaries could be treated as inflow boundaries and make the minimum upwind face limiter meaningless. --- .../TransportModelTracerTofDiscGal.cpp | 104 ++++++++++++++++-- .../TransportModelTracerTofDiscGal.hpp | 14 ++- 2 files changed, 109 insertions(+), 9 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 2ec842f64..de73f47b6 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -133,6 +133,9 @@ namespace Opm : grid_(grid), use_cvi_(use_cvi), use_limiter_(use_limiter), + limiter_relative_flux_threshold_(1e-3), + limiter_method_(MinUpwindFace), + limiter_usage_(DuringComputations), coord_(grid.dimensions), velocity_(grid.dimensions) { @@ -194,6 +197,19 @@ namespace Opm grad_basis_.resize(num_basis*grid_.dimensions); velocity_interpolation_->setupFluxes(darcyflux); 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_); + } } @@ -370,8 +386,8 @@ namespace Opm std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); // Apply limiter. - if (degree_ > 0 && use_limiter_) { - useLimiter(cell); + if (degree_ > 0 && use_limiter_ && limiter_usage_ == DuringComputations) { + applyLimiter(cell, tof_coeff_); } } @@ -389,7 +405,22 @@ namespace Opm - void TransportModelTracerTofDiscGal::useLimiter(const int cell) + void TransportModelTracerTofDiscGal::applyLimiter(const int cell, double* tof) + { + switch (limiter_method_) { + case MinUpwindFace: + applyMinUpwindFaceLimiter(cell, tof); + break; + case MinUpwindAverage: + default: + THROW("Limiter type not implemented: " << limiter_method_); + } + } + + + + + void TransportModelTracerTofDiscGal::applyMinUpwindFaceLimiter(const int cell, double* tof) { if (degree_ != 1) { THROW("This limiter only makes sense for our DG1 implementation."); @@ -399,15 +430,38 @@ namespace Opm // 1. Let M be the minimum TOF value on the upstream faces, // evaluated in the upstream cells. Then the value at all // points in this cell shall be at least M. + // Upstream faces whose flux does not exceed the relative + // flux threshold are not considered for this 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. const int dim = grid_.dimensions; const int num_basis = DGBasis::numBasisFunc(dim, degree_); - double min_upstream_tof = 1e100; double min_here_tof = 1e100; int num_upstream_faces = 0; - // Find minimum tof on upstream 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; @@ -419,7 +473,7 @@ namespace Opm flux = -darcyflux_[face]; upstream_cell = grid_.face_cells[2*face]; } - const bool upstream = (flux < 0.0); + const bool upstream = (flux < -total_flux*limiter_relative_flux_threshold_); if (upstream) { ++num_upstream_faces; } @@ -458,7 +512,7 @@ namespace Opm // Handle by setting a flat solution. std::cout << "Trouble in cell " << cell << std::endl; limiter = 0.0; - tof_coeff_[num_basis*cell] = min_upstream_tof; + tof[num_basis*cell] = min_upstream_tof; } ASSERT(limiter >= 0.0); @@ -466,7 +520,7 @@ namespace Opm if (limiter < 1.0) { std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl; for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) { - tof_coeff_[i] *= limiter; + tof[i] *= limiter; } } else { std::cout << "Not applying limiter in cell " << cell << "!" << std::endl; @@ -476,4 +530,38 @@ namespace Opm + void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess() + { + // Apply the limiter sequentially to all cells. + // This means that a cell's limiting behaviour may be affected by + // any limiting applied to its upstream cells. + const std::vector& seq = TransportModelInterface::sequence(); + const int nc = seq.size(); + ASSERT(nc == grid_.number_of_cells); + for (int i = 0; i < nc; ++i) { + const int cell = seq[i]; + applyLimiter(cell, tof_coeff_); + } + } + + + + + void TransportModelTracerTofDiscGal::applyLimiterAsSimultaneousPostProcess() + { + // Apply the limiter simultaneously to all cells. + // This means that each cell is limited independently from all other cells, + // we write the resulting dofs to a new array instead of writing to tof_coeff_. + // Afterwards we copy the results back to tof_coeff_. + const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_); + std::vector tof_coeffs_new(tof_coeff_, tof_coeff_ + num_basis*grid_.number_of_cells); + for (int c = 0; c < grid_.number_of_cells; ++c) { + applyLimiter(c, &tof_coeffs_new[0]); + } + std::copy(tof_coeffs_new.begin(), tof_coeffs_new.end(), tof_coeff_); + } + + + + } // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 3188b3716..24ad712de 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -88,6 +88,11 @@ namespace Opm boost::shared_ptr velocity_interpolation_; bool use_cvi_; bool use_limiter_; + double limiter_relative_flux_threshold_; + enum LimiterMethod { MinUpwindFace, MinUpwindAverage }; + LimiterMethod limiter_method_; + enum LimiterUsage { DuringComputations, AsPostProcess, AsSimultaneousPostProcess }; + LimiterUsage limiter_usage_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell const double* source_; // one volumetric source term per cell @@ -105,7 +110,14 @@ namespace Opm std::vector velocity_; // Private methods - void useLimiter(const int cell); + + // Apply some limiter, writing to array tof + // (will read data from tof_coeff_, it is ok to call + // with tof_coeff as tof argument. + void applyLimiter(const int cell, double* tof); + void applyMinUpwindFaceLimiter(const int cell, double* tof); + void applyLimiterAsPostProcess(); + void applyLimiterAsSimultaneousPostProcess(); }; } // namespace Opm