Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include <array> 7 : #include <boost/functional/hash.hpp> 8 : #include <iterator> 9 : #include <limits> 10 : #include <memory> 11 : #include <pup.h> 12 : #include <type_traits> 13 : #include <unordered_map> 14 : #include <utility> 15 : 16 : #include "DataStructures/DataVector.hpp" 17 : #include "DataStructures/Tags.hpp" 18 : #include "DataStructures/Tensor/Tensor.hpp" 19 : #include "Domain/Structure/Direction.hpp" 20 : #include "Domain/Structure/DirectionMap.hpp" 21 : #include "Domain/Structure/Element.hpp" 22 : #include "Domain/Structure/ElementId.hpp" 23 : #include "Domain/Structure/OrientationMap.hpp" 24 : #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodHelpers.hpp" 25 : #include "Evolution/DiscontinuousGalerkin/Limiters/MinmodType.hpp" 26 : #include "NumericalAlgorithms/LinearOperators/MeanValue.hpp" 27 : #include "NumericalAlgorithms/Spectral/Mesh.hpp" 28 : #include "Utilities/ErrorHandling/Assert.hpp" 29 : #include "Utilities/Gsl.hpp" 30 : #include "Utilities/Literals.hpp" 31 : #include "Utilities/Numeric.hpp" 32 : #include "Utilities/TMPL.hpp" 33 : 34 : namespace Limiters::Minmod_detail { 35 : 36 : // This function combines the evaluation of the troubled-cell indicator with the 37 : // computation of the post-limiter reduced slopes. The returned bool indicates 38 : // whether the slopes are to be reduced. The slopes themselves are returned by 39 : // pointer. 40 : // 41 : // Note: This function is only made available in this header file to facilitate 42 : // testing. 43 : template <size_t VolumeDim> 44 : bool minmod_limited_slopes( 45 : gsl::not_null<DataVector*> u_lin_buffer, 46 : gsl::not_null<BufferWrapper<VolumeDim>*> buffer, 47 : gsl::not_null<double*> u_mean, 48 : gsl::not_null<std::array<double, VolumeDim>*> u_limited_slopes, 49 : Limiters::MinmodType minmod_type, double tvb_constant, const DataVector& u, 50 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element, 51 : const std::array<double, VolumeDim>& element_size, 52 : const DirectionMap<VolumeDim, double>& effective_neighbor_means, 53 : const DirectionMap<VolumeDim, double>& effective_neighbor_sizes); 54 : 55 : // Implements the minmod limiter for one Tensor<DataVector> at a time. 56 : template <size_t VolumeDim, typename Tag, typename PackagedData> 57 : bool minmod_impl( 58 : const gsl::not_null<DataVector*> u_lin_buffer, 59 : const gsl::not_null<BufferWrapper<VolumeDim>*> buffer, 60 : const gsl::not_null<typename Tag::type*> tensor, 61 : const Limiters::MinmodType minmod_type, const double tvb_constant, 62 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element, 63 : const tnsr::I<DataVector, VolumeDim, Frame::ElementLogical>& logical_coords, 64 : const std::array<double, VolumeDim>& element_size, 65 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData, 66 : boost::hash<DirectionalId<VolumeDim>>>& 67 : neighbor_data) { 68 : // True if the mesh is linear-order in every direction 69 : const bool mesh_is_linear = (mesh.extents() == Index<VolumeDim>(2)); 70 : const bool minmod_type_is_linear = 71 : (minmod_type != Limiters::MinmodType::LambdaPiN); 72 : const bool using_linear_limiter_on_non_linear_mesh = 73 : minmod_type_is_linear and not mesh_is_linear; 74 : 75 : // In each direction, average the size of all different neighbors in that 76 : // direction. Note that only the component of neighor_size that is normal 77 : // to the face is needed (and, therefore, computed). Note that this average 78 : // does not depend on the solution on the neighboring elements, so could be 79 : // precomputed outside of `limit_one_tensor`. Changing the code to 80 : // precompute the average may or may not be a measurable optimization. 81 : const auto effective_neighbor_sizes = 82 : compute_effective_neighbor_sizes(element, neighbor_data); 83 : 84 : bool some_component_was_limited = false; 85 : for (size_t i = 0; i < tensor->size(); ++i) { 86 : // In each direction, average the mean of the i'th tensor component over 87 : // all different neighbors in that direction. This produces one effective 88 : // neighbor per direction. 89 : const auto effective_neighbor_means = 90 : compute_effective_neighbor_means<Tag>(i, element, neighbor_data); 91 : 92 : DataVector& u = (*tensor)[i]; 93 : double u_mean = std::numeric_limits<double>::signaling_NaN(); 94 : std::array<double, VolumeDim> u_limited_slopes{}; 95 : const bool reduce_slopes = minmod_limited_slopes( 96 : u_lin_buffer, buffer, make_not_null(&u_mean), 97 : make_not_null(&u_limited_slopes), minmod_type, tvb_constant, u, mesh, 98 : element, element_size, effective_neighbor_means, 99 : effective_neighbor_sizes); 100 : 101 : if (reduce_slopes or using_linear_limiter_on_non_linear_mesh) { 102 : u = u_mean; 103 : for (size_t d = 0; d < VolumeDim; ++d) { 104 : u += logical_coords.get(d) * gsl::at(u_limited_slopes, d); 105 : } 106 : some_component_was_limited = true; 107 : } 108 : } // end for loop over tensor components 109 : 110 : return some_component_was_limited; 111 : } 112 : 113 : } // namespace Limiters::Minmod_detail