ryujin 2.1.1 revision 9dcb748690310d6a540ebb8b066d1a0834fc7604
postprocessor.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: MIT
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include "postprocessor.h"
9
10#include <simd.h>
11
12#include <deal.II/base/function_parser.h>
13#include <deal.II/numerics/data_out.h>
14#include <deal.II/numerics/vector_tools.h>
15
16#include <atomic>
17#include <chrono>
18#include <fstream>
19
20namespace ryujin
21{
22 template <typename Description, int dim, typename Number>
24 const MPI_Comm &mpi_communicator,
25 const HyperbolicSystem &hyperbolic_system,
26 const OfflineData<dim, Number> &offline_data,
27 const std::string &subsection /*= "Postprocessor"*/)
28 : ParameterAcceptor(subsection)
29 , mpi_communicator_(mpi_communicator)
30 , hyperbolic_system_(&hyperbolic_system)
31 , offline_data_(&offline_data)
32 {
33 beta_ = 10.;
34 add_parameter("schlieren beta",
35 beta_,
36 "Beta factor used in the exponential scale for the schlieren "
37 "and vorticity plots");
38
39 recompute_bounds_ = true;
40 add_parameter(
41 "schlieren recompute bounds",
42 recompute_bounds_,
43 "Recompute bounds for every output cycle. If set to false, bounds once "
44 "at the beginning and reused thereafter.");
45
46 static_assert(HyperbolicSystemView::component_names.size() > 0,
47 "Need at least one scalar quantitity");
48 schlieren_quantities_.push_back(HyperbolicSystemView::component_names[0]);
49
50 add_parameter(
51 "schlieren quantities",
52 schlieren_quantities_,
53 "List of conserved quantities used for the schlieren postprocessor.");
54
55 if constexpr (dim > 1) {
56 add_parameter(
57 "vorticity quantities",
58 vorticity_quantities_,
59 "List of conserved quantities used for the vorticity postprocessor.");
60 }
61 }
62
63
64 template <typename Description, int dim, typename Number>
66 {
67#ifdef DEBUG_OUTPUT
68 std::cout << "Postprocessor<dim, Number>::prepare()" << std::endl;
69#endif
70
71 bounds_.clear();
72 component_names_.clear();
73 schlieren_indices_.clear();
74 vorticity_indices_.clear();
75
76 const auto populate = [&](const auto &strings,
77 auto &indices,
78 const auto &pre) {
79 const auto &cons = HyperbolicSystemView::component_names;
80 const auto &prim = HyperbolicSystemView::primitive_component_names;
81 for (const auto &entry : strings) {
82 bool found = false;
83 for (const auto &[is_primitive, names] :
84 {std::make_pair(false, cons), std::make_pair(true, prim)}) {
85 const auto pos = std::find(std::begin(names), std::end(names), entry);
86 if (!found && pos != std::end(names)) {
87 const auto index = std::distance(std::begin(names), pos);
88 indices.push_back(std::make_pair(is_primitive, index));
89 component_names_.push_back(pre + entry);
90 found = true;
91 }
92 }
93 AssertThrow(
94 found,
95 dealii::ExcMessage("Invalid component name »" + entry + "«"));
96 }
97 };
98 populate(schlieren_quantities_, schlieren_indices_, "schlieren_");
99 populate(vorticity_quantities_, vorticity_indices_, "vorticity_");
100
101 const auto &partitioner = offline_data_->scalar_partitioner();
103 quantities_.resize(component_names_.size());
104 for (auto &it : quantities_)
105 it.reinit(partitioner);
106 }
107
108
109 template <typename Description, int dim, typename Number>
110 void
112 {
113#ifdef DEBUG_OUTPUT
114 std::cout << "Postprocessor<dim, Number>::compute()" << std::endl;
115#endif
116
117 using VA = dealii::VectorizedArray<Number>;
118
119 const auto &affine_constraints = offline_data_->affine_constraints();
120
121 const auto &sparsity_simd = offline_data_->sparsity_pattern_simd();
122 const auto &lumped_mass_matrix = offline_data_->lumped_mass_matrix();
123 const auto &cij_matrix = offline_data_->cij_matrix();
124 const auto &boundary_map = offline_data_->boundary_map();
125
126 const unsigned int n_internal = offline_data_->n_locally_internal();
127 const unsigned int n_owned = offline_data_->n_locally_owned();
128
129 const unsigned int n_schlieren = schlieren_indices_.size();
130 Assert(n_schlieren == schlieren_quantities_.size(),
131 dealii::ExcInternalError());
132 const unsigned int n_vorticities = vorticity_indices_.size();
133 Assert(n_vorticities == vorticity_quantities_.size(),
134 dealii::ExcInternalError());
135 const unsigned int n_quantities = n_schlieren + n_vorticities;
136 Assert(n_quantities == quantities_.size(), dealii::ExcInternalError());
137 Assert(n_quantities == component_names_.size(), dealii::ExcInternalError());
138
139 /*
140 * Step 1: Compute quantities:
141 */
143 {
145
146 auto loop = [&](auto sentinel, unsigned int left, unsigned int right) {
147 using T = decltype(sentinel);
148 unsigned int stride_size = get_stride_size<T>;
149
150 std::vector<grad_type<T>> local_schlieren_values(n_schlieren);
151 std::vector<curl_type<T>> local_vorticity_values(n_vorticities);
152
154 for (unsigned int i = left; i < right; i += stride_size) {
155
156 for (auto &it : local_schlieren_values)
157 it = grad_type<T>();
158 for (auto &it : local_vorticity_values)
159 it = curl_type<T>();
160
161 /* Skip constrained degrees of freedom: */
162 const unsigned int row_length = sparsity_simd.row_length(i);
163 if (row_length == 1)
164 continue;
165
166 const unsigned int *js = sparsity_simd.columns(i);
167 for (unsigned int col_idx = 0; col_idx < row_length;
168 ++col_idx, js += stride_size) {
169
170 const auto U_j = U.template get_tensor<T>(js);
171 const auto view = hyperbolic_system_->template view<dim, T>();
172 const auto prim_j = view.to_primitive_state(U_j);
173
174 const auto c_ij = cij_matrix.template get_tensor<T>(i, col_idx);
175
176 unsigned int k = 0;
177 for (const auto &[is_primitive, index] : schlieren_indices_) {
178 local_schlieren_values[k++] -=
179 c_ij * (is_primitive ? prim_j[index] : U_j[index]);
180 }
181
182 k = 0;
183 for (const auto &[is_primitive, index] : vorticity_indices_) {
184 grad_type<T> q_j;
185 for (unsigned int d = 0; d < dim; ++d)
186 q_j[d] = (is_primitive ? prim_j[index + d] : U_j[index + d]);
187
188 if constexpr (dim == 2) {
189 local_vorticity_values[k++][0] -= cross_product_2d(c_ij) * q_j;
190 } else if constexpr (dim == 3) {
191 local_vorticity_values[k++] -= cross_product_3d(c_ij, q_j);
192 }
193 }
194 }
195
196 /* Fix up boundaries: */
197
198 /* Serialize over the stride: */
199 for (unsigned int k = 0; k < stride_size; ++k) {
200 const auto range = boundary_map.equal_range(i + k);
201 for (auto it = range.first; it != range.second; ++it) {
202 const auto normal = std::get<0>(it->second);
203 const auto id = std::get<3>(it->second);
204
205 if (id == Boundary::slip || id == Boundary::no_slip) {
206 /* Remove normal component of schlieren values at slip
207 * and no-slip boundaries: */
208 for (auto &it : local_schlieren_values) {
209 auto grad_q_i = serialize_tensor(it, k);
210 grad_q_i -= 1. * (grad_q_i * normal) * normal;
211 assign_serial_tensor(it, grad_q_i, k);
212 }
213
214 } else {
215 /* Set schlieren values to zero everywhere else: */
216 for (auto &it : local_schlieren_values) {
218 }
219 }
220
221 for (auto &it : local_vorticity_values) {
222 /* Retain only the normal component of the curl on the
223 * boundary: */
224 if constexpr (dim == 2) {
226 } else if constexpr (dim == 3) {
227 auto curl_q_i = serialize_tensor(it, k);
228 curl_q_i = (curl_q_i * normal) * normal;
229 assign_serial_tensor(it, curl_q_i, k);
230 }
231 }
232 }
233 }
234
235 /* Populate quantities: */
236
237 const auto m_i = load_value<T>(lumped_mass_matrix, i);
238
239 unsigned int k = 0;
240
241 for (const auto &schlieren : local_schlieren_values) {
242 const auto value_i = schlieren.norm() / m_i;
243 store_value<T>(quantities_[k++], value_i, i);
244 }
245
246 for (const auto &vorticity : local_vorticity_values) {
247 auto value_i =
248 (dim == 2 ? vorticity[0] / m_i : vorticity.norm() / m_i);
249 store_value<T>(quantities_[k++], value_i, i);
250 }
251 } /* i */
252 };
253
254 /* Parallel non-vectorized loop: */
255 loop(Number(), n_internal, n_owned);
256 /* Parallel vectorized SIMD loop: */
257 loop(VA(), 0, n_internal);
258
260 }
261
262 /*
263 * Step 2: Compute bounds and synchronize over MPI ranks:
264 */
265
266 /* Force recomputation of bounds: */
267 if (recompute_bounds_)
268 bounds_.clear();
269
270 if (bounds_.size() != n_quantities) {
271 bounds_.clear();
272 bounds_.resize(
273 n_quantities,
274 std::make_pair(Number(0.), std::numeric_limits<Number>::max()));
275
276 for (unsigned int d = 0; d < n_quantities; ++d) {
277 auto &[q_max, q_min] = bounds_[d];
278 for (unsigned int i = 0; i < n_owned; ++i) {
279 const auto q = quantities_[d].local_element(i);
280 q_max = std::max(q_max, std::abs(q));
281 q_min = std::min(q_min, std::abs(q));
282 }
283 q_max = dealii::Utilities::MPI::max(q_max, mpi_communicator_);
284 q_min = dealii::Utilities::MPI::min(q_min, mpi_communicator_);
285 Assert(q_max >= q_min, dealii::ExcInternalError());
286 }
287 }
288
289 /*
290 * Step 3: Normalize quantities on exponential scale:
291 */
292
293 {
294 for (unsigned int d = 0; d < n_quantities; ++d) {
295 auto &[q_max, q_min] = bounds_[d];
296 for (unsigned int i = 0; i < n_owned; ++i) {
297 auto &q = quantities_[d].local_element(i);
298 constexpr auto eps = std::numeric_limits<Number>::epsilon();
299 const auto magnitude =
300 Number(1.) -
301 std::exp(-beta_ * (std::abs(q) - q_min) / (q_max - q_min + eps));
302 q = std::copysign(magnitude, q);
303 }
304 }
305 }
306
307 /*
308 * Step 4: Fix up constraints and distribute:
309 */
310
311 for (auto &it : quantities_) {
312 affine_constraints.distribute(it);
313 it.update_ghost_values();
314 }
315 }
316
317} // namespace ryujin
typename Description::HyperbolicSystem HyperbolicSystem
Definition: postprocessor.h:44
void compute(const vector_type &U) const
Postprocessor(const MPI_Comm &mpi_communicator, const HyperbolicSystem &hyperbolic_system, const OfflineData< dim, Number > &offline_data, const std::string &subsection="/Postprocessor")
#define RYUJIN_PARALLEL_REGION_BEGIN
Definition: openmp.h:54
#define RYUJIN_OMP_FOR
Definition: openmp.h:70
#define RYUJIN_PARALLEL_REGION_END
Definition: openmp.h:63
DEAL_II_ALWAYS_INLINE void assign_serial_tensor(dealii::Tensor< rank, dim, dealii::VectorizedArray< Number, width > > &result, const dealii::Tensor< rank, dim, Number > &serial, const unsigned int k)
Definition: simd.h:453
DEAL_II_ALWAYS_INLINE dealii::Tensor< rank, dim, Number > serialize_tensor(const dealii::Tensor< rank, dim, dealii::VectorizedArray< Number, width > > &vectorized, const unsigned int k)
Definition: simd.h:411