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 <cstddef>
9 : #include <limits>
10 : #include <string>
11 : #include <type_traits>
12 : #include <unordered_map>
13 : #include <utility>
14 :
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "DataStructures/Variables.hpp"
17 : #include "Domain/Structure/Direction.hpp"
18 : #include "Domain/Structure/Element.hpp"
19 : #include "Domain/Structure/ElementId.hpp"
20 : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoGridHelpers.hpp"
21 : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoHelpers.hpp"
22 : #include "Evolution/DiscontinuousGalerkin/Limiters/WenoOscillationIndicator.hpp"
23 : #include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp"
24 : #include "NumericalAlgorithms/LinearOperators/MeanValue.hpp"
25 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
26 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
27 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
28 : #include "Utilities/Gsl.hpp"
29 :
30 : namespace Limiters::Weno_detail {
31 :
32 : // Compute the Simple WENO solution for one tensor component
33 : //
34 : // This interface is intended for use in limiters that check the troubled-cell
35 : // indicator independently for each tensor component. These limiters generally
36 : // need to limit only a subset of the tensor components.
37 : //
38 : // When calling `simple_weno_impl`,
39 : // - `interpolator_buffer` may be empty
40 : // - `modified_neighbor_solution_buffer` should contain one DataVector for each
41 : // neighboring element (i.e. for each entry in `neighbor_data`)
42 : template <typename Tag, size_t VolumeDim, typename PackagedData>
43 : void simple_weno_impl(
44 : const gsl::not_null<std::unordered_map<
45 : DirectionalId<VolumeDim>, intrp::RegularGrid<VolumeDim>,
46 : boost::hash<DirectionalId<VolumeDim>>>*>
47 : interpolator_buffer,
48 : const gsl::not_null<
49 : std::unordered_map<DirectionalId<VolumeDim>, DataVector,
50 : boost::hash<DirectionalId<VolumeDim>>>*>
51 : modified_neighbor_solution_buffer,
52 : const gsl::not_null<typename Tag::type*> tensor,
53 : const double neighbor_linear_weight, const size_t tensor_storage_index,
54 : const Mesh<VolumeDim>& mesh, const Element<VolumeDim>& element,
55 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
56 : boost::hash<DirectionalId<VolumeDim>>>&
57 : neighbor_data) {
58 : // Check that basis is LGL or LG
59 : // Note that the SimpleWeno implementation should generalize well to other
60 : // bases beyond LGL or LG, once the oscillation_indicator function has been
61 : // generalized.
62 : ASSERT(mesh.basis() == make_array<VolumeDim>(Spectral::Basis::Legendre),
63 : "Unsupported basis: " << mesh);
64 : ASSERT(mesh.quadrature() ==
65 : make_array<VolumeDim>(Spectral::Quadrature::GaussLobatto) or
66 : mesh.quadrature() ==
67 : make_array<VolumeDim>(Spectral::Quadrature::Gauss),
68 : "Unsupported quadrature: " << mesh);
69 :
70 : ASSERT(
71 : modified_neighbor_solution_buffer->size() == neighbor_data.size(),
72 : "modified_neighbor_solution_buffer->size() = "
73 : << modified_neighbor_solution_buffer->size()
74 : << "\nneighbor_data.size() = " << neighbor_data.size()
75 : << "\nmodified_neighbor_solution_buffer was incorrectly initialized "
76 : "before calling simple_weno_impl.");
77 :
78 : // Compute the modified neighbor solutions.
79 : // First extrapolate neighbor data onto local grid points, then shift the
80 : // extrapolated data so its mean matches the local mean.
81 : DataVector& component_to_limit = (*tensor)[tensor_storage_index];
82 : const double local_mean = mean_value(component_to_limit, mesh);
83 : for (const auto& neighbor_and_data : neighbor_data) {
84 : const auto& neighbor = neighbor_and_data.first;
85 : const auto& data = neighbor_and_data.second;
86 :
87 : if (interpolator_buffer->find(neighbor) == interpolator_buffer->end()) {
88 : // No interpolator found => create one
89 : const auto& direction = neighbor.direction();
90 : const auto& source_mesh = data.mesh;
91 : const auto target_1d_logical_coords =
92 : Weno_detail::local_grid_points_in_neighbor_logical_coords(
93 : mesh, source_mesh, element, direction);
94 : interpolator_buffer->insert(std::make_pair(
95 : neighbor, intrp::RegularGrid<VolumeDim>(source_mesh, mesh,
96 : target_1d_logical_coords)));
97 : }
98 :
99 : // Avoid allocations by working directly in the preallocated buffer
100 : DataVector& buffer = modified_neighbor_solution_buffer->at(neighbor);
101 :
102 : interpolator_buffer->at(neighbor).interpolate(
103 : make_not_null(&buffer),
104 : get<Tag>(data.volume_data)[tensor_storage_index]);
105 : const double neighbor_mean = mean_value(buffer, mesh);
106 : buffer += (local_mean - neighbor_mean);
107 : }
108 :
109 : // Sum local and modified neighbor polynomials for the WENO reconstruction
110 : Weno_detail::reconstruct_from_weighted_sum(
111 : make_not_null(&component_to_limit), neighbor_linear_weight,
112 : Weno_detail::DerivativeWeight::PowTwoEll, mesh,
113 : *modified_neighbor_solution_buffer);
114 : }
115 :
116 : // Compute the Simple WENO solution for one tensor
117 : //
118 : // This interface is intended for use in limiters that check the troubled-cell
119 : // indicator for the whole cell, and apply the limiter to all fields.
120 : //
121 : // When calling `simple_weno_impl`,
122 : // - `interpolator_buffer` may be empty
123 : // - `modified_neighbor_solution_buffer` should contain one DataVector for each
124 : // neighboring element (i.e. for each entry in `neighbor_data`)
125 : template <typename Tag, size_t VolumeDim, typename PackagedData>
126 : void simple_weno_impl(
127 : const gsl::not_null<std::unordered_map<
128 : DirectionalId<VolumeDim>, intrp::RegularGrid<VolumeDim>,
129 : boost::hash<DirectionalId<VolumeDim>>>*>
130 : interpolator_buffer,
131 : const gsl::not_null<
132 : std::unordered_map<DirectionalId<VolumeDim>, DataVector,
133 : boost::hash<DirectionalId<VolumeDim>>>*>
134 : modified_neighbor_solution_buffer,
135 : const gsl::not_null<typename Tag::type*> tensor,
136 : const double neighbor_linear_weight, const Mesh<VolumeDim>& mesh,
137 : const Element<VolumeDim>& element,
138 : const std::unordered_map<DirectionalId<VolumeDim>, PackagedData,
139 : boost::hash<DirectionalId<VolumeDim>>>&
140 : neighbor_data) {
141 : for (size_t tensor_storage_index = 0; tensor_storage_index < tensor->size();
142 : ++tensor_storage_index) {
143 : simple_weno_impl<Tag>(interpolator_buffer,
144 : modified_neighbor_solution_buffer, tensor,
145 : neighbor_linear_weight, tensor_storage_index, mesh,
146 : element, neighbor_data);
147 : }
148 : }
149 :
150 : } // namespace Limiters::Weno_detail
|
