12#include <deal.II/base/function_parser.h>
13#include <deal.II/numerics/data_out.h>
14#include <deal.II/numerics/vector_tools.h>
18 template <
typename Description,
int dim,
typename Number>
20 const MPI_Comm &mpi_communicator,
23 const std::string &subsection )
24 : ParameterAcceptor(subsection)
25 , mpi_communicator_(mpi_communicator)
26 , hyperbolic_system_(&hyperbolic_system)
27 , offline_data_(&offline_data)
30 add_parameter(
"schlieren beta",
32 "Beta factor used in the exponential scale for the schlieren "
33 "and vorticity plots");
35 recompute_bounds_ =
true;
37 "schlieren recompute bounds",
39 "Recompute bounds for every output cycle. If set to false, bounds once "
40 "at the beginning and reused thereafter.");
42 static_assert(View::component_names.size() > 0,
43 "Need at least one scalar quantitity");
44 schlieren_quantities_.push_back(View::component_names[0]);
47 "schlieren quantities",
48 schlieren_quantities_,
49 "List of conserved quantities used for the schlieren postprocessor.");
51 if constexpr (dim > 1) {
53 "vorticity quantities",
54 vorticity_quantities_,
55 "List of conserved quantities used for the vorticity postprocessor.");
60 template <
typename Description,
int dim,
typename Number>
64 std::cout <<
"Postprocessor<dim, Number>::prepare()" << std::endl;
68 component_names_.clear();
69 schlieren_indices_.clear();
70 vorticity_indices_.clear();
72 const auto populate = [&](
const auto &strings,
75 const auto &cons = View::component_names;
76 const auto &prim = View::primitive_component_names;
77 for (
const auto &entry : strings) {
79 for (
const auto &[is_primitive, names] :
80 {std::make_pair(
false, cons), std::make_pair(
true, prim)}) {
81 const auto pos = std::find(std::begin(names), std::end(names), entry);
82 if (!found && pos != std::end(names)) {
83 const auto index = std::distance(std::begin(names), pos);
84 indices.push_back(std::make_pair(is_primitive, index));
85 component_names_.push_back(pre + entry);
91 dealii::ExcMessage(
"Invalid component name »" + entry +
"«"));
94 populate(schlieren_quantities_, schlieren_indices_,
"schlieren_");
95 populate(vorticity_quantities_, vorticity_indices_,
"vorticity_");
97 const auto &partitioner = offline_data_->scalar_partitioner();
99 quantities_.resize(component_names_.size());
100 for (
auto &it : quantities_)
101 it.reinit(partitioner);
105 template <
typename Description,
int dim,
typename Number>
110 std::cout <<
"Postprocessor<dim, Number>::compute()" << std::endl;
113 const auto &U = std::get<0>(state_vector);
115 using VA = dealii::VectorizedArray<Number>;
117 const auto &affine_constraints = offline_data_->affine_constraints();
119 const auto &sparsity_simd = offline_data_->sparsity_pattern_simd();
120 const auto &lumped_mass_matrix = offline_data_->lumped_mass_matrix();
121 const auto &cij_matrix = offline_data_->cij_matrix();
122 const auto &boundary_map = offline_data_->boundary_map();
124 const unsigned int n_internal = offline_data_->n_locally_internal();
125 const unsigned int n_owned = offline_data_->n_locally_owned();
127 const unsigned int n_schlieren = schlieren_indices_.size();
128 Assert(n_schlieren == schlieren_quantities_.size(),
129 dealii::ExcInternalError());
130 const unsigned int n_vorticities = vorticity_indices_.size();
131 Assert(n_vorticities == vorticity_quantities_.size(),
132 dealii::ExcInternalError());
133 const unsigned int n_quantities = n_schlieren + n_vorticities;
134 Assert(n_quantities == quantities_.size(), dealii::ExcInternalError());
135 Assert(n_quantities == component_names_.size(), dealii::ExcInternalError());
144 auto loop = [&](
auto sentinel,
unsigned int left,
unsigned int right) {
145 using T =
decltype(sentinel);
146 unsigned int stride_size = get_stride_size<T>;
148 std::vector<grad_type<T>> local_schlieren_values(n_schlieren);
149 std::vector<curl_type<T>> local_vorticity_values(n_vorticities);
152 for (
unsigned int i = left; i < right; i += stride_size) {
154 for (
auto &it : local_schlieren_values)
156 for (
auto &it : local_vorticity_values)
160 const unsigned int row_length = sparsity_simd.row_length(i);
164 const unsigned int *js = sparsity_simd.columns(i);
165 for (
unsigned int col_idx = 0; col_idx < row_length;
166 ++col_idx, js += stride_size) {
168 const auto U_j = U.template get_tensor<T>(js);
169 const auto view = hyperbolic_system_->template view<dim, T>();
170 const auto prim_j = view.to_primitive_state(U_j);
172 const auto c_ij = cij_matrix.template get_tensor<T>(i, col_idx);
175 for (
const auto &[is_primitive, index] : schlieren_indices_) {
176 local_schlieren_values[k++] -=
177 c_ij * (is_primitive ? prim_j[index] : U_j[index]);
181 for (
const auto &[is_primitive, index] : vorticity_indices_) {
183 for (
unsigned int d = 0; d < dim; ++d)
184 q_j[d] = (is_primitive ? prim_j[index + d] : U_j[index + d]);
186 if constexpr (dim == 2) {
187 local_vorticity_values[k++][0] -= cross_product_2d(c_ij) * q_j;
188 }
else if constexpr (dim == 3) {
189 local_vorticity_values[k++] -= cross_product_3d(c_ij, q_j);
197 for (
unsigned int k = 0; k < stride_size; ++k) {
198 const auto range = boundary_map.equal_range(i + k);
199 for (
auto it = range.first; it != range.second; ++it) {
200 const auto normal = std::get<0>(it->second);
201 const auto id = std::get<3>(it->second);
206 for (
auto &it : local_schlieren_values) {
208 grad_q_i -= 1. * (grad_q_i * normal) * normal;
214 for (
auto &it : local_schlieren_values) {
219 for (
auto &it : local_vorticity_values) {
222 if constexpr (dim == 2) {
224 }
else if constexpr (dim == 3) {
226 curl_q_i = (curl_q_i * normal) * normal;
235 const auto m_i = get_entry<T>(lumped_mass_matrix, i);
239 for (
const auto &schlieren : local_schlieren_values) {
240 const auto value_i = schlieren.norm() / m_i;
241 write_entry<T>(quantities_[k++], value_i, i);
244 for (
const auto &vorticity : local_vorticity_values) {
246 (dim == 2 ? vorticity[0] / m_i : vorticity.norm() / m_i);
247 write_entry<T>(quantities_[k++], value_i, i);
253 loop(Number(), n_internal, n_owned);
255 loop(VA(), 0, n_internal);
265 if (recompute_bounds_)
268 if (bounds_.size() != n_quantities) {
272 std::make_pair(Number(0.), std::numeric_limits<Number>::max()));
274 for (
unsigned int d = 0; d < n_quantities; ++d) {
275 auto &[q_max, q_min] = bounds_[d];
276 for (
unsigned int i = 0; i < n_owned; ++i) {
277 const auto q = quantities_[d].local_element(i);
278 q_max = std::max(q_max, std::abs(q));
279 q_min = std::min(q_min, std::abs(q));
281 q_max = dealii::Utilities::MPI::max(q_max, mpi_communicator_);
282 q_min = dealii::Utilities::MPI::min(q_min, mpi_communicator_);
283 Assert(q_max >= q_min, dealii::ExcInternalError());
292 for (
unsigned int d = 0; d < n_quantities; ++d) {
293 auto &[q_max, q_min] = bounds_[d];
294 for (
unsigned int i = 0; i < n_owned; ++i) {
295 auto &q = quantities_[d].local_element(i);
296 constexpr auto eps = std::numeric_limits<Number>::epsilon();
297 const auto magnitude =
299 std::exp(-beta_ * (std::abs(q) - q_min) / (q_max - q_min + eps));
300 q = std::copysign(magnitude, q);
309 for (
auto &it : quantities_) {
310 affine_constraints.distribute(it);
311 it.update_ghost_values();
typename Description::HyperbolicSystem HyperbolicSystem
dealii::Tensor< 1, dim==2 ? 1 :dim, T > curl_type
dealii::Tensor< 1, dim, T > grad_type
void compute(const StateVector &state_vector) const
typename View::StateVector StateVector
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)