12#include <deal.II/base/function_parser.h>
13#include <deal.II/numerics/data_out.h>
14#include <deal.II/numerics/vector_tools.h>
22 template <
typename Description,
int dim,
typename Number>
24 const MPI_Comm &mpi_communicator,
27 const std::string &subsection )
28 : ParameterAcceptor(subsection)
29 , mpi_communicator_(mpi_communicator)
30 , hyperbolic_system_(&hyperbolic_system)
31 , offline_data_(&offline_data)
34 add_parameter(
"schlieren beta",
36 "Beta factor used in the exponential scale for the schlieren "
37 "and vorticity plots");
39 recompute_bounds_ =
true;
41 "schlieren recompute bounds",
43 "Recompute bounds for every output cycle. If set to false, bounds once "
44 "at the beginning and reused thereafter.");
46 static_assert(HyperbolicSystemView::component_names.size() > 0,
47 "Need at least one scalar quantitity");
48 schlieren_quantities_.push_back(HyperbolicSystemView::component_names[0]);
51 "schlieren quantities",
52 schlieren_quantities_,
53 "List of conserved quantities used for the schlieren postprocessor.");
55 if constexpr (dim > 1) {
57 "vorticity quantities",
58 vorticity_quantities_,
59 "List of conserved quantities used for the vorticity postprocessor.");
64 template <
typename Description,
int dim,
typename Number>
68 std::cout <<
"Postprocessor<dim, Number>::prepare()" << std::endl;
72 component_names_.clear();
73 schlieren_indices_.clear();
74 vorticity_indices_.clear();
76 const auto populate = [&](
const auto &strings,
79 const auto &cons = HyperbolicSystemView::component_names;
80 const auto &prim = HyperbolicSystemView::primitive_component_names;
81 for (
const auto &entry : strings) {
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);
95 dealii::ExcMessage(
"Invalid component name »" + entry +
"«"));
98 populate(schlieren_quantities_, schlieren_indices_,
"schlieren_");
99 populate(vorticity_quantities_, vorticity_indices_,
"vorticity_");
101 const auto &partitioner = offline_data_->scalar_partitioner();
103 quantities_.resize(component_names_.size());
104 for (
auto &it : quantities_)
105 it.reinit(partitioner);
109 template <
typename Description,
int dim,
typename Number>
114 std::cout <<
"Postprocessor<dim, Number>::compute()" << std::endl;
117 using VA = dealii::VectorizedArray<Number>;
119 const auto &affine_constraints = offline_data_->affine_constraints();
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();
126 const unsigned int n_internal = offline_data_->n_locally_internal();
127 const unsigned int n_owned = offline_data_->n_locally_owned();
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());
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>;
150 std::vector<grad_type<T>> local_schlieren_values(n_schlieren);
151 std::vector<curl_type<T>> local_vorticity_values(n_vorticities);
154 for (
unsigned int i = left; i < right; i += stride_size) {
156 for (
auto &it : local_schlieren_values)
158 for (
auto &it : local_vorticity_values)
162 const unsigned int row_length = sparsity_simd.row_length(i);
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) {
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);
174 const auto c_ij = cij_matrix.template get_tensor<T>(i, col_idx);
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]);
183 for (
const auto &[is_primitive, index] : vorticity_indices_) {
185 for (
unsigned int d = 0; d < dim; ++d)
186 q_j[d] = (is_primitive ? prim_j[index + d] : U_j[index + d]);
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);
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);
208 for (
auto &it : local_schlieren_values) {
210 grad_q_i -= 1. * (grad_q_i * normal) * normal;
216 for (
auto &it : local_schlieren_values) {
221 for (
auto &it : local_vorticity_values) {
224 if constexpr (dim == 2) {
226 }
else if constexpr (dim == 3) {
228 curl_q_i = (curl_q_i * normal) * normal;
237 const auto m_i = load_value<T>(lumped_mass_matrix, i);
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);
246 for (
const auto &vorticity : local_vorticity_values) {
248 (dim == 2 ? vorticity[0] / m_i : vorticity.norm() / m_i);
249 store_value<T>(quantities_[k++], value_i, i);
255 loop(Number(), n_internal, n_owned);
257 loop(VA(), 0, n_internal);
267 if (recompute_bounds_)
270 if (bounds_.size() != n_quantities) {
274 std::make_pair(Number(0.), std::numeric_limits<Number>::max()));
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));
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());
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 =
301 std::exp(-beta_ * (std::abs(q) - q_min) / (q_max - q_min + eps));
302 q = std::copysign(magnitude, q);
311 for (
auto &it : quantities_) {
312 affine_constraints.distribute(it);
313 it.update_ghost_values();
typename Description::HyperbolicSystem HyperbolicSystem
dealii::Tensor< 1, dim==2 ? 1 :dim, T > curl_type
void compute(const vector_type &U) const
dealii::Tensor< 1, dim, T > grad_type
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
#define RYUJIN_PARALLEL_REGION_END
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)
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)