deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-76.h
Go to the documentation of this file.
1true);
282
283 for (unsigned int face = range.first; face < range.second; ++face)
284 {
285 phi_m.reinit(face);
286 phi_m.gather_evaluate(src, face_evaluation_flags);
287 phi_p.reinit(face);
288 phi_p.gather_evaluate(src, face_evaluation_flags);
289
290 // some operations on the face quadrature points
291
292 phi_m.integrate_scatter(face_evaluation_flags, dst);
293 phi_p.integrate_scatter(face_evaluation_flags, dst);
294 }
295 },
296 [&](const auto &data, auto &dst, const auto &src, const auto range) {
297 // operation performed boundary faces
298
300
301 for (unsigned int face = range.first; face < range.second; ++face)
302 {
303 phi_m.reinit(face);
304 phi_m.gather_evaluate(src, face_evaluation_flags);
305
306 // some operations on the face quadrature points
307
308 phi_m.integrate_scatter(face_evaluation_flags, dst);
309 }
310 },
311 dst,
312 src);
313@endcode
314
315in the following way:
316
317@code
319 [&](const auto &data, auto &dst, const auto &src, const auto range) {
323
324 for (unsigned int cell = range.first; cell < range.second; ++cell)
325 {
326 phi.reinit(cell);
327 phi.gather_evaluate(src, cell_evaluation_flags);
328
329 // some operations on the cell quadrature points
330
331 phi.integrate_scatter(cell_evaluation_flags, dst);
332
333 // loop over all faces of cell
334 for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
335 {
336 if (data.get_faces_by_cells_boundary_id(cell, face)[0] ==
338 {
339 // internal face
340 phi_m.reinit(cell, face);
341 phi_m.gather_evaluate(src, face_evaluation_flags);
342 phi_p.reinit(cell, face);
343 phi_p.gather_evaluate(src, face_evaluation_flags);
344
345 // some operations on the face quadrature points
346
347 phi_m.integrate_scatter(face_evaluation_flags, dst);
348 }
349 else
350 {
351 // boundary face
352 phi_m.reinit(cell, face);
353 phi_m.gather_evaluate(src, face_evaluation_flags);
354
355 // some operations on the face quadrature points
356
357 phi_m.integrate_scatter(face_evaluation_flags, dst);
358 }
359 }
360 }
361 },
362 dst,
363 src);
364@endcode
365
367the cell number and the local face number. The given example only
372
375
376@code
377typename MatrixFree<dim, Number>::AdditionalData additional_data;
378
379// set flags as usual (not shown)
380
381additional_data.hold_all_faces_to_owned_cells = true;
382additional_data.mapping_update_flags_faces_by_cells =
383 additional_data.mapping_update_flags_inner_faces |
384 additional_data.mapping_update_flags_boundary_faces;
385
386data.reinit(mapping, dof_handler, constraint, quadrature, additional_data);
387@endcode
388
390for all faces of the cells.
391
392Currently, cell-centric loops in deal.II only work for uniformly refined meshes
393and if no constraints are applied (which is the standard case DG is normally
394used).
395
396
397<a name="step_76-ProvidinglambdastoMatrixFreeloops"></a><h3>Providing lambdas to MatrixFree loops</h3>
398
399
402a version where a class and a pointer to one of its methods are used and a
403variant where lambdas are utilized.
404
407
408@code
409void
410local_apply_cell(const MatrixFree<dim, Number> & data,
411 VectorType & dst,
412 const VectorType & src,
413 const std::pair<unsigned int, unsigned int> &range) const
414{
416 for (unsigned int cell = range.first; cell < range.second; ++cell)
417 {
418 phi.reinit(cell);
419 phi.gather_evaluate(src, cell_evaluation_flags);
420
421 // some operations on the quadrature points
422
423 phi.integrate_scatter(cell_evaluation_flags, dst);
424 }
425}
426@endcode
427
428@code
429matrix_free.cell_loop(&Operator::local_apply_cell, this, dst, src);
430@endcode
431
432However, it is also possible to pass an anonymous function via a lambda function
434
435@code
436matrix_free.template cell_loop<VectorType, VectorType>(
437 [&](const auto &data, auto &dst, const auto &src, const auto range) {
439 for (unsigned int cell = range.first; cell < range.second; ++cell)
440 {
441 phi.reinit(cell);
442 phi.gather_evaluate(src, cell_evaluation_flags);
443
444 // some operations on the quadrature points
445
446 phi.integrate_scatter(cell_evaluation_flags, dst);
447 }
448 },
449 dst,
450 src);
451@endcode
452
453<a name="step_76-VectorizedArrayType"></a><h3>VectorizedArrayType</h3>
454
455
456The class VectorizedArray<Number> is a key component to achieve the high
458It is a wrapper class around a short vector of @f$n@f$ entries of type Number and
460(SIMD) concepts by intrinsic functions. The length of the vector can be
463
464In the default case (<code>VectorizedArray<Number></code>), the vector length is
465set at compile time of the library to
468specified as <code>VectorizedArray<Number, size></code>, where <code>size</code> explicitly
470set. A full list of supported vector lengths is presented in the following table:
471
472<table align="center" class="doxtable">
473 <tr>
474 <th>double</th>
475 <th>float</th>
476 <th>ISA</th>
477 </tr>
478 <tr>
479 <td><code>VectorizedArray<double, 1></code></td>
480 <td><code>VectorizedArray<float, 1></code></td>
481 <td>(auto-vectorization)</td>
482 </tr>
483 <tr>
484 <td><code>VectorizedArray<double, 2></code></td>
485 <td><code>VectorizedArray<float, 4></code></td>
486 <td>SSE2/AltiVec</td>
487 </tr>
488 <tr>
489 <td><code>VectorizedArray<double, 4></code></td>
490 <td><code>VectorizedArray<float, 8></code></td>
491 <td>AVX/AVX2</td>
492 </tr>
493 <tr>
494 <td><code>VectorizedArray<double, 8></code></td>
495 <td><code>VectorizedArray<float, 16></code></td>
496 <td>AVX-512</td>
497 </tr>
498</table>
499
500This allows users to select the vector length/ISA and, as a consequence, the
501number of cells to be processed at once in matrix-free operator evaluations,
503degrees (and dimensions).
504
507cells, one can concentrate on a single cell.
508
510a matching interface. Specifically, this prepares deal.II for the <code>std::simd</code>
513
514
515<table align="center" class="doxtable">
516 <tr>
517 <th>VectorizedArray (deal.II)</th>
518 <th>std::simd (C++23)</th>
519 </tr>
520 <tr>
521 <td><code>VectorizedArray<Number></code></td>
522 <td><code>std::experimental::native_simd<Number></code></td>
523 </tr>
524 <tr>
525 <td><code>VectorizedArray<Number, size></code></td>
526 <td><code>std::experimental::fixed_size_simd<Number, size></code></td>
527 </tr>
528</table>
529 *
530 *
531 * <a name="step_76-CommProg"></a>
532 * <h1> The commented program</h1>
533 *
534 *
535 * <a name="step_76-Parametersandutilityfunctions"></a>
537 *
538
539 *
540 * The same includes as in @ref step_67 "step-67":
541 *
542 * @code
543 *   #include <deal.II/base/conditional_ostream.h>
544 *   #include <deal.II/base/function.h>
545 *   #include <deal.II/base/time_stepping.h>
546 *   #include <deal.II/base/timer.h>
547 *   #include <deal.II/base/utilities.h>
548 *   #include <deal.II/base/vectorization.h>
549 *  
550 *   #include <deal.II/distributed/tria.h>
551 *  
552 *   #include <deal.II/dofs/dof_handler.h>
553 *  
554 *   #include <deal.II/fe/fe_dgq.h>
555 *   #include <deal.II/fe/fe_system.h>
556 *  
557 *   #include <deal.II/grid/grid_generator.h>
558 *   #include <deal.II/grid/tria.h>
559 *   #include <deal.II/grid/tria_accessor.h>
560 *   #include <deal.II/grid/tria_iterator.h>
561 *  
562 *   #include <deal.II/lac/affine_constraints.h>
563 *   #include <deal.II/lac/la_parallel_vector.h>
564 *  
565 *   #include <deal.II/matrix_free/fe_evaluation.h>
566 *   #include <deal.II/matrix_free/matrix_free.h>
567 *   #include <deal.II/matrix_free/operators.h>
568 *  
569 *   #include <deal.II/numerics/data_out.h>
570 *  
571 *   #include <fstream>
572 *   #include <iomanip>
573 *   #include <iostream>
574 *  
575 * @endcode
576 *
577 * A new include for categorizing of cells according to their boundary IDs:
578 *
579 * @code
580 *   #include <deal.II/matrix_free/tools.h>
581 *  
582 *  
583 *  
585 *   {
586 *   using namespace dealii;
587 *  
588 * @endcode
589 *
590 * The same input parameters as in @ref step_67 "step-67":
591 *
592 * @code
593 *   constexpr unsigned int testcase = 1;
594 *   constexpr unsigned int dimension = 2;
595 *   constexpr unsigned int n_global_refinements = 2;
596 *   constexpr unsigned int fe_degree = 5;
597 *   constexpr unsigned int n_q_points_1d = fe_degree + 2;
598 *  
599 * @endcode
600 *
601 * This parameter specifies the size of the shared-memory group. Currently,
604 * having access to the same shared-memory domain are grouped together.
605 *
606 * @code
607 *   constexpr unsigned int group_size = numbers::invalid_unsigned_int;
608 *  
609 *   using Number = double;
610 *  
611 * @endcode
612 *
614 * default case, VectorizedArray<Number> is used, i.e., the highest
616 * the maximum number of vector lanes is used. However, one might reduce
617 * the number of filled lanes, e.g., by writing
618 * <code>using VectorizedArrayType = VectorizedArray<Number, 4></code> to only
619 * process 4 cells.
620 *
621 * @code
622 *   using VectorizedArrayType = VectorizedArray<Number>;
623 *  
624 * @endcode
625 *
626 * The following parameters have not changed:
627 *
628 * @code
629 *   constexpr double gamma = 1.4;
630 *   constexpr double final_time = testcase == 0 ? 10 : 2.0;
631 *   constexpr double output_tick = testcase == 0 ? 1 : 0.05;
632 *  
633 *   const double courant_number = 0.15 / std::pow(fe_degree, 1.5);
634 *  
635 * @endcode
636 *
637 * Specify max number of time steps useful for performance studies.
638 *
639 * @code
640 *   constexpr unsigned int max_time_steps = numbers::invalid_unsigned_int;
641 *  
642 * @endcode
643 *
645 * with the purpose to minimize global vector access:
646 *
647 * @code
649 *   {
654 *   };
656 *  
657 *  
658 *  
660 *   {
661 *   public:
663 *   {
665 *   switch (scheme)
666 *   {
667 *   case stage_3_order_3:
669 *   break;
670 *   case stage_5_order_4:
672 *   break;
673 *   case stage_7_order_4:
675 *   break;
676 *   case stage_9_order_5:
678 *   break;
679 *  
680 *   default:
681 *   AssertThrow(false, ExcNotImplemented());
682 *   }
686 *   std::vector<double> ci; // not used
687 *   rk_integrator.get_coefficients(ai, bi, ci);
688 *   }
689 *  
690 *   unsigned int n_stages() const
691 *   {
692 *   return bi.size();
693 *   }
694 *  
695 *   template <typename VectorType, typename Operator>
697 *   const double current_time,
698 *   const double time_step,
699 *   VectorType &solution,
700 *   VectorType &vec_ri,
701 *   VectorType &vec_ki) const
702 *   {
703 *   AssertDimension(ai.size() + 1, bi.size());
704 *  
705 *   vec_ki.swap(solution);
706 *  
707 *   double sum_previous_bi = 0;
708 *   for (unsigned int stage = 0; stage < bi.size(); ++stage)
709 *   {
710 *   const double c_i = stage == 0 ? 0 : sum_previous_bi + ai[stage - 1];
711 *  
712 *   pde_operator.perform_stage(stage,
713 *   current_time + c_i * time_step,
714 *   bi[stage] * time_step,
715 *   (stage == bi.size() - 1 ?
716 *   0 :
717 *   ai[stage] * time_step),
718 *   (stage % 2 == 0 ? vec_ki : vec_ri),
719 *   (stage % 2 == 0 ? vec_ri : vec_ki),
720 *   solution);
721 *  
722 *   if (stage > 0)
723 *   sum_previous_bi += bi[stage - 1];
724 *   }
725 *   }
726 *  
727 *   private:
728 *   std::vector<double> bi;
729 *   std::vector<double> ai;
730 *   };
731 *  
732 *  
733 * @endcode
734 *
735 * Euler-specific utility functions from @ref step_67 "step-67":
736 *
737 * @code
739 *   {
742 *   };
744 *  
745 *  
746 *  
747 *   template <int dim>
748 *   class ExactSolution : public Function<dim>
749 *   {
750 *   public:
751 *   ExactSolution(const double time)
752 *   : Function<dim>(dim + 2, time)
753 *   {}
754 *  
755 *   virtual double value(const Point<dim> &p,
756 *   const unsigned int component = 0) const override;
757 *   };
758 *  
759 *  
760 *  
761 *   template <int dim>
763 *   const unsigned int component) const
764 *   {
765 *   const double t = this->get_time();
766 *  
767 *   switch (testcase)
768 *   {
769 *   case 0:
770 *   {
771 *   Assert(dim == 2, ExcNotImplemented());
772 *   const double beta = 5;
773 *  
775 *   x0[0] = 5.;
776 *   const double radius_sqr =
777 *   (x - x0).norm_square() - 2. * (x[0] - x0[0]) * t + t * t;
778 *   const double factor =
779 *   beta / (numbers::PI * 2) * std::exp(1. - radius_sqr);
780 *   const double density_log = std::log2(
781 *   std::abs(1. - (gamma - 1.) / gamma * 0.25 * factor * factor));
782 *   const double density = std::exp2(density_log * (1. / (gamma - 1.)));
783 *   const double u = 1. - factor * (x[1] - x0[1]);
784 *   const double v = factor * (x[0] - t - x0[0]);
785 *  
786 *   if (component == 0)
787 *   return density;
788 *   else if (component == 1)
789 *   return density * u;
790 *   else if (component == 2)
791 *   return density * v;
792 *   else
793 *   {
794 *   const double pressure =
795 *   std::exp2(density_log * (gamma / (gamma - 1.)));
796 *   return pressure / (gamma - 1.) +
797 *   0.5 * (density * u * u + density * v * v);
798 *   }
799 *   }
800 *  
801 *   case 1:
802 *   {
803 *   if (component == 0)
804 *   return 1.;
805 *   else if (component == 1)
806 *   return 0.4;
807 *   else if (component == dim + 1)
808 *   return 3.097857142857143;
809 *   else
810 *   return 0.;
811 *   }
812 *  
813 *   default:
815 *   return 0.;
816 *   }
817 *   }
818 *  
819 *  
820 *  
821 *   template <int dim, typename Number>
825 *   {
826 *   const Number inverse_density = Number(1.) / conserved_variables[0];
827 *  
829 *   for (unsigned int d = 0; d < dim; ++d)
831 *  
832 *   return velocity;
833 *   }
834 *  
835 *   template <int dim, typename Number>
837 *   Number
839 *   {
842 *  
843 *   Number rho_u_dot_u = conserved_variables[1] * velocity[0];
844 *   for (unsigned int d = 1; d < dim; ++d)
846 *  
847 *   return (gamma - 1.) * (conserved_variables[dim + 1] - 0.5 * rho_u_dot_u);
848 *   }
849 *  
850 *   template <int dim, typename Number>
854 *   {
858 *  
860 *   for (unsigned int d = 0; d < dim; ++d)
861 *   {
862 *   flux[0][d] = conserved_variables[1 + d];
863 *   for (unsigned int e = 0; e < dim; ++e)
864 *   flux[e + 1][d] = conserved_variables[e + 1] * velocity[d];
865 *   flux[d + 1][d] += pressure;
866 *   flux[dim + 1][d] =
867 *   velocity[d] * (conserved_variables[dim + 1] + pressure);
868 *   }
869 *  
870 *   return flux;
871 *   }
872 *  
873 *   template <int n_components, int dim, typename Number>
876 *   operator*(const Tensor<1, n_components, Tensor<1, dim, Number>> &matrix,
877 *   const Tensor<1, dim, Number> &vector)
878 *   {
880 *   for (unsigned int d = 0; d < n_components; ++d)
881 *   result[d] = matrix[d] * vector;
882 *   return result;
883 *   }
884 *  
885 *   template <int dim, typename Number>
890 *   const Tensor<1, dim, Number> &normal)
891 *   {
892 *   const auto velocity_m = euler_velocity<dim>(u_m);
893 *   const auto velocity_p = euler_velocity<dim>(u_p);
894 *  
895 *   const auto pressure_m = euler_pressure<dim>(u_m);
896 *   const auto pressure_p = euler_pressure<dim>(u_p);
897 *  
898 *   const auto flux_m = euler_flux<dim>(u_m);
899 *   const auto flux_p = euler_flux<dim>(u_p);
900 *  
901 *   switch (numerical_flux_type)
902 *   {
904 *   {
905 *   const auto lambda =
907 *   gamma * pressure_p * (1. / u_p[0]),
909 *   gamma * pressure_m * (1. / u_m[0])));
910 *  
911 *   return 0.5 * (flux_m * normal + flux_p * normal) +
912 *   0.5 * lambda * (u_m - u_p);
913 *   }
914 *  
916 *   {
917 *   const auto avg_velocity_normal =
918 *   0.5 * ((velocity_m + velocity_p) * normal);
919 *   const auto avg_c = std::sqrt(std::abs(
920 *   0.5 * gamma *
921 *   (pressure_p * (1. / u_p[0]) + pressure_m * (1. / u_m[0]))));
922 *   const Number s_pos =
923 *   std::max(Number(), avg_velocity_normal + avg_c);
924 *   const Number s_neg =
925 *   std::min(Number(), avg_velocity_normal - avg_c);
926 *   const Number inverse_s = Number(1.) / (s_pos - s_neg);
927 *  
928 *   return inverse_s *
929 *   ((s_pos * (flux_m * normal) - s_neg * (flux_p * normal)) -
930 *   s_pos * s_neg * (u_m - u_p));
931 *   }
932 *  
933 *   default:
934 *   {
936 *   return {};
937 *   }
938 *   }
939 *   }
940 *  
941 *  
942 *  
943 * @endcode
944 *
945 * General-purpose utility functions from @ref step_67 "step-67":
946 *
947 * @code
948 *   template <int dim, typename VectorizedArrayType>
949 *   VectorizedArrayType
950 *   evaluate_function(const Function<dim> &function,
952 *   const unsigned int component)
953 *   {
954 *   VectorizedArrayType result;
955 *   for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
956 *   {
957 *   Point<dim> p;
958 *   for (unsigned int d = 0; d < dim; ++d)
959 *   p[d] = p_vectorized[d][v];
960 *   result[v] = function.value(p, component);
961 *   }
962 *   return result;
963 *   }
964 *  
965 *  
966 *   template <int dim, typename VectorizedArrayType, int n_components = dim + 2>
968 *   evaluate_function(const Function<dim> &function,
970 *   {
971 *   AssertDimension(function.n_components, n_components);
973 *   for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
974 *   {
975 *   Point<dim> p;
976 *   for (unsigned int d = 0; d < dim; ++d)
977 *   p[d] = p_vectorized[d][v];
978 *   for (unsigned int d = 0; d < n_components; ++d)
979 *   result[d][v] = function.value(p, d);
980 *   }
981 *   return result;
982 *   }
983 *  
984 *  
985 * @endcode
986 *
987 *
988 * <a name="step_76-EuleroperatorusingacellcentricloopandMPI30sharedmemory"></a>
989 * <h3>Euler operator using a cell-centric loop and MPI-3.0 shared memory</h3>
990 *
991
992 *
993 * Euler operator from @ref step_67 "step-67" with some changes as detailed below:
994 *
995 * @code
996 *   template <int dim, int degree, int n_points_1d>
997 *   class EulerOperator
998 *   {
999 *   public:
1000 *   static constexpr unsigned int n_quadrature_points_1d = n_points_1d;
1001 *  
1003 *  
1004 *   ~EulerOperator();
1005 *  
1006 *   void reinit(const Mapping<dim> &mapping,
1007 *   const DoFHandler<dim> &dof_handler);
1008 *  
1009 *   void set_inflow_boundary(const types::boundary_id boundary_id,
1010 *   std::unique_ptr<Function<dim>> inflow_function);
1011 *  
1013 *   const types::boundary_id boundary_id,
1014 *   std::unique_ptr<Function<dim>> outflow_energy);
1015 *  
1016 *   void set_wall_boundary(const types::boundary_id boundary_id);
1017 *  
1018 *   void set_body_force(std::unique_ptr<Function<dim>> body_force);
1019 *  
1020 *   void
1021 *   perform_stage(const unsigned int stage,
1022 *   const Number cur_time,
1023 *   const Number bi,
1024 *   const Number ai,
1028 *  
1029 *   void project(const Function<dim> &function,
1031 *  
1032 *   std::array<double, 3> compute_errors(
1033 *   const Function<dim> &function,
1034 *   const LinearAlgebra::distributed::Vector<Number> &solution) const;
1035 *  
1037 *   const LinearAlgebra::distributed::Vector<Number> &solution) const;
1038 *  
1039 *   void
1041 *  
1042 *   private:
1043 * @endcode
1044 *
1047 * shared-memory capabilities:
1048 *
1049 * @code
1051 *  
1052 *   MatrixFree<dim, Number, VectorizedArrayType> data;
1053 *  
1054 *   TimerOutput &timer;
1055 *  
1056 *   std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1058 *   std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1060 *   std::set<types::boundary_id> wall_boundaries;
1061 *   std::unique_ptr<Function<dim>> body_force;
1062 *   };
1063 *  
1064 *  
1065 *  
1066 * @endcode
1067 *
1068 * New constructor, which creates a sub-communicator. The user can specify
1069 * the size of the sub-communicator via the global parameter group_size. If
1070 * the size is set to -1, all MPI processes of a
1071 * shared-memory domain are combined to a group. The specified size is
1073 * and, therefore, setting the <code>size</code> to <code>-1</code> is a
1075 * disable the MPI-3.0 shared-memory features of MatrixFree and rely
1077 * <code>MPI_Irecv</code>.
1078 *
1079 * @code
1080 *   template <int dim, int degree, int n_points_1d>
1081 *   EulerOperator<dim, degree, n_points_1d>::EulerOperator(TimerOutput &timer)
1082 *   : timer(timer)
1083 *   {
1085 *   if (group_size == 1)
1086 *   {
1087 *   this->subcommunicator = MPI_COMM_SELF;
1088 *   }
1090 *   {
1092 *  
1095 *   rank,
1098 *   }
1099 *   else
1100 *   {
1102 *   }
1103 *   #else
1104 *   (void)subcommunicator;
1105 *   (void)group_size;
1106 *   this->subcommunicator = MPI_COMM_SELF;
1107 *   #endif
1108 *   }
1109 *  
1110 *  
1111 * @endcode
1112 *
1113 * New destructor responsible for freeing of the sub-communicator.
1114 *
1115 * @code
1116 *   template <int dim, int degree, int n_points_1d>
1118 *   {
1120 *   if (this->subcommunicator != MPI_COMM_SELF)
1122 *   #endif
1123 *   }
1124 *  
1125 *  
1126 * @endcode
1127 *
1128 * Modified reinit() function to set up the internal data structures in
1130 * the MPI-3.0 shared-memory capabilities are used:
1131 *
1132 * @code
1133 *   template <int dim, int degree, int n_points_1d>
1134 *   void EulerOperator<dim, degree, n_points_1d>::reinit(
1135 *   const Mapping<dim> &mapping,
1136 *   const DoFHandler<dim> &dof_handler)
1137 *   {
1138 *   const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler};
1139 *   const AffineConstraints<double> dummy;
1140 *   const std::vector<const AffineConstraints<double> *> constraints = {&dummy};
1141 *   const std::vector<Quadrature<1>> quadratures = {QGauss<1>(n_q_points_1d),
1142 *   QGauss<1>(fe_degree + 1)};
1143 *  
1145 *   additional_data;
1146 *   additional_data.mapping_update_flags =
1149 *   additional_data.mapping_update_flags_inner_faces =
1152 *   additional_data.mapping_update_flags_boundary_faces =
1155 *   additional_data.tasks_parallel_scheme =
1157 *  
1158 * @endcode
1159 *
1160 * Categorize cells so that all lanes have the same boundary IDs for each
1165 *
1166 * @code
1167 *   MatrixFreeTools::categorize_by_boundary_ids(dof_handler.get_triangulation(),
1168 *   additional_data);
1169 *  
1170 * @endcode
1171 *
1172 * Enable MPI-3.0 shared-memory capabilities within MatrixFree by providing
1173 * the sub-communicator:
1174 *
1175 * @code
1176 *   additional_data.communicator_sm = subcommunicator;
1177 *  
1178 *   data.reinit(
1179 *   mapping, dof_handlers, constraints, quadratures, additional_data);
1180 *   }
1181 *  
1182 *  
1183 * @endcode
1184 *
1187 * compared to @ref step_67 "step-67".
1188 *
1189
1190 *
1191 * In contrast to @ref step_67 "step-67", we are not executing the advection step
1192 * (using MatrixFree::loop()) and the inverse mass-matrix step
1193 * (using MatrixFree::cell_loop()) in sequence, but evaluate everything in
1194 * one go inside of MatrixFree::loop_cell_centric(). This function expects
1195 * a single function that is executed on each locally-owned (macro) cell as
1196 * parameter so that we need to loop over all faces of that cell and perform
1197 * needed integration steps on our own.
1198 *
1199
1200 *
1202 * functions from @ref step_67 "step-67" so that comments related the evaluation of the weak
1203 * form are skipped here:
1204 * - <code>EulerDG::EulerOperator::local_apply_cell</code>
1208 *
1209 * @code
1210 *   template <int dim, int degree, int n_points_1d>
1211 *   void EulerOperator<dim, degree, n_points_1d>::perform_stage(
1212 *   const unsigned int stage,
1213 *   const Number current_time,
1214 *   const Number bi,
1215 *   const Number ai,
1216 *   const LinearAlgebra::distributed::Vector<Number> &current_ri,
1217 *   LinearAlgebra::distributed::Vector<Number> &vec_ki,
1218 *   LinearAlgebra::distributed::Vector<Number> &solution) const
1219 *   {
1220 *   for (auto &i : inflow_boundaries)
1221 *   i.second->set_time(current_time);
1222 *   for (auto &i : subsonic_outflow_boundaries)
1223 *   i.second->set_time(current_time);
1224 *  
1225 * @endcode
1226 *
1228 * providing a lambda containing the effects of the cell, face and
1229 * boundary-face integrals:
1230 *
1231 * @code
1232 *   data.template loop_cell_centric<LinearAlgebra::distributed::Vector<Number>,
1233 *   LinearAlgebra::distributed::Vector<Number>>(
1234 *   [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
1235 *   using FECellIntegral = FEEvaluation<dim,
1236 *   degree,
1237 *   n_points_1d,
1238 *   dim + 2,
1239 *   Number,
1240 *   VectorizedArrayType>;
1241 *   using FEFaceIntegral = FEFaceEvaluation<dim,
1242 *   degree,
1243 *   n_points_1d,
1244 *   dim + 2,
1245 *   Number,
1246 *   VectorizedArrayType>;
1247 *  
1250 *   FEFaceIntegral phi_m(data, true);
1251 *   FEFaceIntegral phi_p(data, false);
1252 *  
1255 *   dynamic_cast<Functions::ConstantFunction<dim> *>(body_force.get());
1256 *  
1257 *   if (constant_function)
1261 *  
1264 *   dim,
1265 *   n_points_1d,
1266 *   n_points_1d,
1267 *   VectorizedArrayType,
1268 *   Number>
1269 *   eval({},
1270 *   data.get_shape_info().data[0].shape_gradients_collocation_eo,
1271 *   {});
1272 *  
1274 *   phi.n_components);
1275 *  
1276 * @endcode
1277 *
1278 * Loop over all cell batches:
1279 *
1280 * @code
1281 *   for (unsigned int cell = cell_range.first; cell < cell_range.second;
1282 *   ++cell)
1283 *   {
1284 *   phi.reinit(cell);
1285 *  
1286 *   if (ai != Number())
1287 *   phi_temp.reinit(cell);
1288 *  
1289 * @endcode
1290 *
1291 * Read values from global vector and compute the values at the
1292 * quadrature points:
1293 *
1294 * @code
1295 *   if (ai != Number() && stage == 0)
1296 *   {
1297 *   phi.read_dof_values(src);
1298 *  
1299 *   for (unsigned int i = 0;
1300 *   i < phi.static_dofs_per_component * (dim + 2);
1301 *   ++i)
1302 *   phi_temp.begin_dof_values()[i] = phi.begin_dof_values()[i];
1303 *  
1304 *   phi.evaluate(EvaluationFlags::values);
1305 *   }
1306 *   else
1307 *   {
1308 *   phi.gather_evaluate(src, EvaluationFlags::values);
1309 *   }
1310 *  
1311 * @endcode
1312 *
1313 * Buffer the computed values at the quadrature points, since
1315 * step, however, are needed later on for the face integrals:
1316 *
1317 * @code
1318 *   for (unsigned int i = 0; i < phi.static_n_q_points * (dim + 2); ++i)
1319 *   buffer[i] = phi.begin_values()[i];
1320 *  
1321 * @endcode
1322 *
1323 * Apply the cell integral at the cell quadrature points. See also
1324 * the function <code>EulerOperator::local_apply_cell()</code> from
1325 * @ref step_67 "step-67":
1326 *
1327 * @code
1328 *   for (const unsigned int q : phi.quadrature_point_indices())
1329 *   {
1330 *   const auto w_q = phi.get_value(q);
1331 *   phi.submit_gradient(euler_flux<dim>(w_q), q);
1332 *   if (body_force.get() != nullptr)
1333 *   {
1338 *   *body_force, phi.quadrature_point(q));
1339 *  
1341 *   for (unsigned int d = 0; d < dim; ++d)
1342 *   forcing[d + 1] = w_q[0] * force[d];
1343 *   for (unsigned int d = 0; d < dim; ++d)
1344 *   forcing[dim + 1] += force[d] * w_q[d + 1];
1345 *  
1346 *   phi.submit_value(forcing, q);
1347 *   }
1348 *   }
1349 *  
1350 * @endcode
1351 *
1352 * Test with the gradient of the test functions in the quadrature
1353 * points. We skip the interpolation back to the support points
1354 * of the element, since we first collect all contributions in the
1355 * cell quadrature points and only perform the interpolation back
1356 * as the final step.
1357 *
1358 * @code
1359 *   {
1360 *   auto *values_ptr = phi.begin_values();
1361 *   auto *gradient_ptr = phi.begin_gradients();
1362 *  
1363 *   for (unsigned int c = 0; c < dim + 2; ++c)
1364 *   {
1365 *   if (dim >= 1 && body_force.get() == nullptr)
1367 *   values_ptr);
1368 *   else if (dim >= 1)
1370 *   values_ptr);
1371 *   if (dim >= 2)
1373 *   1,
1374 *   values_ptr);
1375 *   if (dim >= 3)
1377 *   2,
1378 *   values_ptr);
1379 *  
1380 *   values_ptr += phi.static_n_q_points;
1381 *   gradient_ptr += phi.static_n_q_points * dim;
1382 *   }
1383 *   }
1384 *  
1385 * @endcode
1386 *
1387 * Loop over all faces of the current cell:
1388 *
1389 * @code
1390 *   for (unsigned int face = 0;
1392 *   ++face)
1393 *   {
1394 * @endcode
1395 *
1396 * Determine the boundary ID of the current face. Since we have
1397 * set up MatrixFree in a way that all filled lanes have
1398 * guaranteed the same boundary ID, we can select the
1399 * boundary ID of the first lane.
1400 *
1401 * @code
1402 *   const auto boundary_ids =
1403 *   data.get_faces_by_cells_boundary_id(cell, face);
1404 *  
1405 *   Assert(std::equal(boundary_ids.begin(),
1406 *   boundary_ids.begin() +
1407 *   data.n_active_entries_per_cell_batch(cell),
1408 *   boundary_ids.begin()),
1409 *   ExcMessage("Boundary IDs of lanes differ."));
1410 *  
1411 *   const auto boundary_id = boundary_ids[0];
1412 *  
1413 *   phi_m.reinit(cell, face);
1414 *  
1415 * @endcode
1416 *
1417 * Interpolate the values from the cell quadrature points to the
1418 * quadrature points of the current face via a simple 1d
1419 * interpolation:
1420 *
1421 * @code
1423 *   n_points_1d - 1,
1424 *   VectorizedArrayType>::
1426 *   dim + 2,
1428 *   data.get_shape_info(),
1429 *   buffer.data(),
1430 *   phi_m.begin_values(),
1431 *   face);
1432 *  
1433 * @endcode
1434 *
1435 * Check if the face is an internal or a boundary face and
1436 * select a different code path based on this information:
1437 *
1438 * @code
1439 *   if (boundary_id == numbers::internal_face_boundary_id)
1440 *   {
1441 * @endcode
1442 *
1443 * Process and internal face. The following lines of code
1444 * are a copy of the function
1445 * <code>EulerDG::EulerOperator::local_apply_face</code>
1446 * from @ref step_67 "step-67":
1447 *
1448 * @code
1449 *   phi_p.reinit(cell, face);
1450 *   phi_p.gather_evaluate(src, EvaluationFlags::values);
1451 *  
1452 *   for (const unsigned int q :
1453 *   phi_m.quadrature_point_indices())
1454 *   {
1456 *   euler_numerical_flux<dim>(phi_m.get_value(q),
1457 *   phi_p.get_value(q),
1458 *   phi_m.normal_vector(q));
1459 *   phi_m.submit_value(-numerical_flux, q);
1460 *   }
1461 *   }
1462 *   else
1463 *   {
1464 * @endcode
1465 *
1466 * Process a boundary face. These following lines of code
1467 * are a copy of the function
1468 * <code>EulerDG::EulerOperator::local_apply_boundary_face</code>
1469 * from @ref step_67 "step-67":
1470 *
1471 * @code
1472 *   for (const unsigned int q :
1473 *   phi_m.quadrature_point_indices())
1474 *   {
1475 *   const auto w_m = phi_m.get_value(q);
1476 *   const auto normal = phi_m.normal_vector(q);
1477 *  
1478 *   auto rho_u_dot_n = w_m[1] * normal[0];
1479 *   for (unsigned int d = 1; d < dim; ++d)
1480 *   rho_u_dot_n += w_m[1 + d] * normal[d];
1481 *  
1482 *   bool at_outflow = false;
1483 *  
1485 *  
1486 *   if (wall_boundaries.find(boundary_id) !=
1487 *   wall_boundaries.end())
1488 *   {
1489 *   w_p[0] = w_m[0];
1490 *   for (unsigned int d = 0; d < dim; ++d)
1491 *   w_p[d + 1] =
1492 *   w_m[d + 1] - 2. * rho_u_dot_n * normal[d];
1493 *   w_p[dim + 1] = w_m[dim + 1];
1494 *   }
1495 *   else if (inflow_boundaries.find(boundary_id) !=
1496 *   inflow_boundaries.end())
1498 *   *inflow_boundaries.find(boundary_id)->second,
1499 *   phi_m.quadrature_point(q));
1500 *   else if (subsonic_outflow_boundaries.find(
1501 *   boundary_id) !=
1503 *   {
1504 *   w_p = w_m;
1505 *   w_p[dim + 1] =
1507 *   .find(boundary_id)
1508 *   ->second,
1509 *   phi_m.quadrature_point(q),
1510 *   dim + 1);
1511 *   at_outflow = true;
1512 *   }
1513 *   else
1514 *   AssertThrow(false,
1515 *   ExcMessage(
1516 *   "Unknown boundary id, did "
1517 *   "you set a boundary condition for "
1518 *   "this part of the domain boundary?"));
1519 *  
1520 *   auto flux = euler_numerical_flux<dim>(w_m, w_p, normal);
1521 *  
1522 *   if (at_outflow)
1523 *   for (unsigned int v = 0;
1524 *   v < VectorizedArrayType::size();
1525 *   ++v)
1526 *   {
1527 *   if (rho_u_dot_n[v] < -1e-12)
1528 *   for (unsigned int d = 0; d < dim; ++d)
1529 *   flux[d + 1][v] = 0.;
1530 *   }
1531 *  
1532 *   phi_m.submit_value(-flux, q);
1533 *   }
1534 *   }
1535 *  
1536 * @endcode
1537 *
1538 * Evaluate local integrals related to cell by quadrature and
1539 * add into cell contribution via a simple 1d interpolation:
1540 *
1541 * @code
1543 *   n_points_1d - 1,
1544 *   VectorizedArrayType>::
1546 *   dim + 2,
1548 *   data.get_shape_info(),
1549 *   phi_m.begin_values(),
1550 *   phi.begin_values(),
1551 *   face);
1552 *   }
1553 *  
1554 * @endcode
1555 *
1556 * Apply inverse mass matrix in the cell quadrature points. See
1557 * also the function
1558 * <code>EulerDG::EulerOperator::local_apply_inverse_mass_matrix()</code>
1559 * from @ref step_67 "step-67":
1560 *
1561 * @code
1562 *   for (unsigned int q = 0; q < phi.static_n_q_points; ++q)
1563 *   {
1564 *   const auto factor = VectorizedArrayType(1.0) / phi.JxW(q);
1565 *   for (unsigned int c = 0; c < dim + 2; ++c)
1566 *   phi.begin_values()[c * phi.static_n_q_points + q] =
1567 *   phi.begin_values()[c * phi.static_n_q_points + q] * factor;
1568 *   }
1569 *  
1570 * @endcode
1571 *
1573 * Gauss-Lobatto space:
1574 *
1575 * @code
1579 *   dim,
1580 *   degree + 1,
1581 *   n_points_1d>::do_backward(dim + 2,
1582 *   data.get_shape_info()
1583 *   .data[0]
1584 *   .inverse_shape_values_eo,
1585 *   false,
1586 *   phi.begin_values(),
1587 *   phi.begin_dof_values());
1588 *  
1589 * @endcode
1590 *
1591 * Perform Runge-Kutta update and write results back to global
1592 * vectors:
1593 *
1594 * @code
1595 *   if (ai == Number())
1596 *   {
1597 *   for (unsigned int q = 0; q < phi.static_dofs_per_cell; ++q)
1598 *   phi.begin_dof_values()[q] = bi * phi.begin_dof_values()[q];
1599 *   phi.distribute_local_to_global(solution);
1600 *   }
1601 *   else
1602 *   {
1603 *   if (stage != 0)
1604 *   phi_temp.read_dof_values(solution);
1605 *  
1606 *   for (unsigned int q = 0; q < phi.static_dofs_per_cell; ++q)
1607 *   {
1608 *   const auto K_i = phi.begin_dof_values()[q];
1609 *  
1610 *   phi.begin_dof_values()[q] =
1611 *   phi_temp.begin_dof_values()[q] + (ai * K_i);
1612 *  
1613 *   phi_temp.begin_dof_values()[q] += bi * K_i;
1614 *   }
1615 *   phi.set_dof_values(dst);
1616 *   phi_temp.set_dof_values(solution);
1617 *   }
1618 *   }
1619 *   },
1620 *   vec_ki,
1621 *   current_ri,
1622 *   true,
1624 *   }
1625 *  
1626 *  
1627 *  
1628 * @endcode
1629 *
1630 * From here, the code of @ref step_67 "step-67" has not changed.
1631 *
1632 * @code
1633 *   template <int dim, int degree, int n_points_1d>
1636 *   {
1637 *   data.initialize_dof_vector(vector);
1638 *   }
1639 *  
1640 *  
1641 *  
1642 *   template <int dim, int degree, int n_points_1d>
1644 *   const types::boundary_id boundary_id,
1645 *   std::unique_ptr<Function<dim>> inflow_function)
1646 *   {
1647 *   AssertThrow(subsonic_outflow_boundaries.find(boundary_id) ==
1649 *   wall_boundaries.find(boundary_id) == wall_boundaries.end(),
1650 *   ExcMessage("You already set the boundary with id " +
1651 *   std::to_string(static_cast<int>(boundary_id)) +
1652 *   " to another type of boundary before now setting " +
1653 *   "it as inflow"));
1654 *   AssertThrow(inflow_function->n_components == dim + 2,
1655 *   ExcMessage("Expected function with dim+2 components"));
1656 *  
1658 *   }
1659 *  
1660 *  
1661 *  
1662 *   template <int dim, int degree, int n_points_1d>
1664 *   const types::boundary_id boundary_id,
1665 *   std::unique_ptr<Function<dim>> outflow_function)
1666 *   {
1667 *   AssertThrow(inflow_boundaries.find(boundary_id) ==
1668 *   inflow_boundaries.end() &&
1669 *   wall_boundaries.find(boundary_id) == wall_boundaries.end(),
1670 *   ExcMessage("You already set the boundary with id " +
1671 *   std::to_string(static_cast<int>(boundary_id)) +
1672 *   " to another type of boundary before now setting " +
1673 *   "it as subsonic outflow"));
1674 *   AssertThrow(outflow_function->n_components == dim + 2,
1675 *   ExcMessage("Expected function with dim+2 components"));
1676 *  
1678 *   }
1679 *  
1680 *  
1681 *  
1682 *   template <int dim, int degree, int n_points_1d>
1684 *   const types::boundary_id boundary_id)
1685 *   {
1686 *   AssertThrow(inflow_boundaries.find(boundary_id) ==
1687 *   inflow_boundaries.end() &&
1688 *   subsonic_outflow_boundaries.find(boundary_id) ==
1690 *   ExcMessage("You already set the boundary with id " +
1691 *   std::to_string(static_cast<int>(boundary_id)) +
1692 *   " to another type of boundary before now setting " +
1693 *   "it as wall boundary"));
1694 *  
1695 *   wall_boundaries.insert(boundary_id);
1696 *   }
1697 *  
1698 *  
1699 *  
1700 *   template <int dim, int degree, int n_points_1d>
1702 *   std::unique_ptr<Function<dim>> body_force)
1703 *   {
1704 *   AssertDimension(body_force->n_components, dim);
1705 *  
1706 *   this->body_force = std::move(body_force);
1707 *   }
1708 *  
1709 *  
1710 *  
1711 *   template <int dim, int degree, int n_points_1d>
1713 *   const Function<dim> &function,
1715 *   {
1717 *   phi(data, 0, 1);
1719 *   degree,
1720 *   dim + 2,
1721 *   Number,
1722 *   VectorizedArrayType>
1723 *   inverse(phi);
1724 *   solution.zero_out_ghost_values();
1725 *   for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
1726 *   {
1727 *   phi.reinit(cell);
1728 *   for (const unsigned int q : phi.quadrature_point_indices())
1729 *   phi.submit_dof_value(evaluate_function(function,
1730 *   phi.quadrature_point(q)),
1731 *   q);
1732 *   inverse.transform_from_q_points_to_basis(dim + 2,
1733 *   phi.begin_dof_values(),
1734 *   phi.begin_dof_values());
1735 *   phi.set_dof_values(solution);
1736 *   }
1737 *   }
1738 *  
1739 *  
1740 *  
1741 *   template <int dim, int degree, int n_points_1d>
1743 *   const Function<dim> &function,
1744 *   const LinearAlgebra::distributed::Vector<Number> &solution) const
1745 *   {
1746 *   TimerOutput::Scope t(timer, "compute errors");
1747 *   double errors_squared[3] = {};
1749 *   phi(data, 0, 0);
1750 *  
1751 *   for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
1752 *   {
1753 *   phi.reinit(cell);
1754 *   phi.gather_evaluate(solution, EvaluationFlags::values);
1755 *   VectorizedArrayType local_errors_squared[3] = {};
1756 *   for (const unsigned int q : phi.quadrature_point_indices())
1757 *   {
1758 *   const auto error =
1759 *   evaluate_function(function, phi.quadrature_point(q)) -
1760 *   phi.get_value(q);
1761 *   const auto JxW = phi.JxW(q);
1762 *  
1763 *   local_errors_squared[0] += error[0] * error[0] * JxW;
1764 *   for (unsigned int d = 0; d < dim; ++d)
1765 *   local_errors_squared[1] += (error[d + 1] * error[d + 1]) * JxW;
1766 *   local_errors_squared[2] += (error[dim + 1] * error[dim + 1]) * JxW;
1767 *   }
1768 *   for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
1769 *   ++v)
1770 *   for (unsigned int d = 0; d < 3; ++d)
1772 *   }
1773 *  
1775 *  
1776 *   std::array<double, 3> errors;
1777 *   for (unsigned int d = 0; d < 3; ++d)
1779 *  
1780 *   return errors;
1781 *   }
1782 *  
1783 *  
1784 *  
1785 *   template <int dim, int degree, int n_points_1d>
1787 *   const LinearAlgebra::distributed::Vector<Number> &solution) const
1788 *   {
1789 *   TimerOutput::Scope t(timer, "compute transport speed");
1790 *   Number max_transport = 0;
1792 *   phi(data, 0, 1);
1793 *  
1794 *   for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
1795 *   {
1796 *   phi.reinit(cell);
1797 *   phi.gather_evaluate(solution, EvaluationFlags::values);
1798 *   VectorizedArrayType local_max = 0.;
1799 *   for (const unsigned int q : phi.quadrature_point_indices())
1800 *   {
1801 *   const auto solution = phi.get_value(q);
1802 *   const auto velocity = euler_velocity<dim>(solution);
1803 *   const auto pressure = euler_pressure<dim>(solution);
1804 *  
1805 *   const auto inverse_jacobian = phi.inverse_jacobian(q);
1806 *   const auto convective_speed = inverse_jacobian * velocity;
1807 *   VectorizedArrayType convective_limit = 0.;
1808 *   for (unsigned int d = 0; d < dim; ++d)
1811 *  
1812 *   const auto speed_of_sound =
1813 *   std::sqrt(gamma * pressure * (1. / solution[0]));
1814 *  
1816 *   for (unsigned int d = 0; d < dim; ++d)
1817 *   eigenvector[d] = 1.;
1818 *   for (unsigned int i = 0; i < 5; ++i)
1819 *   {
1820 *   eigenvector = transpose(inverse_jacobian) *
1821 *   (inverse_jacobian * eigenvector);
1822 *   VectorizedArrayType eigenvector_norm = 0.;
1823 *   for (unsigned int d = 0; d < dim; ++d)
1827 *   }
1828 *   const auto jac_times_ev = inverse_jacobian * eigenvector;
1829 *   const auto max_eigenvalue = std::sqrt(
1831 *   local_max =
1833 *   max_eigenvalue * speed_of_sound + convective_limit);
1834 *   }
1835 *  
1836 *   for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
1837 *   ++v)
1839 *   }
1840 *  
1842 *  
1843 *   return max_transport;
1844 *   }
1845 *  
1846 *  
1847 *  
1848 *   template <int dim>
1849 *   class EulerProblem
1850 *   {
1851 *   public:
1852 *   EulerProblem();
1853 *  
1854 *   void run();
1855 *  
1856 *   private:
1857 *   void make_grid_and_dofs();
1858 *  
1859 *   void output_results(const unsigned int result_number);
1860 *  
1862 *  
1864 *  
1867 *   #else
1869 *   #endif
1870 *  
1871 *   const FESystem<dim> fe;
1872 *   const MappingQ<dim> mapping;
1873 *   DoFHandler<dim> dof_handler;
1874 *  
1875 *   TimerOutput timer;
1876 *  
1878 *  
1879 *   double time, time_step;
1880 *  
1881 *   class Postprocessor : public DataPostprocessor<dim>
1882 *   {
1883 *   public:
1884 *   Postprocessor();
1885 *  
1886 *   virtual void evaluate_vector_field(
1888 *   std::vector<Vector<double>> &computed_quantities) const override;
1889 *  
1890 *   virtual std::vector<std::string> get_names() const override;
1891 *  
1892 *   virtual std::vector<
1894 *   get_data_component_interpretation() const override;
1895 *  
1896 *   virtual UpdateFlags get_needed_update_flags() const override;
1897 *  
1898 *   private:
1899 *   const bool do_schlieren_plot;
1900 *   };
1901 *   };
1902 *  
1903 *  
1904 *  
1905 *   template <int dim>
1907 *   : do_schlieren_plot(dim == 2)
1908 *   {}
1909 *  
1910 *  
1911 *  
1912 *   template <int dim>
1915 *   std::vector<Vector<double>> &computed_quantities) const
1916 *   {
1917 *   const unsigned int n_evaluation_points = inputs.solution_values.size();
1918 *  
1919 *   if (do_schlieren_plot == true)
1920 *   Assert(inputs.solution_gradients.size() == n_evaluation_points,
1921 *   ExcInternalError());
1922 *  
1924 *   ExcInternalError());
1925 *   Assert(inputs.solution_values[0].size() == dim + 2, ExcInternalError());
1927 *   dim + 2 + (do_schlieren_plot == true ? 1 : 0),
1928 *   ExcInternalError());
1929 *  
1930 *   for (unsigned int p = 0; p < n_evaluation_points; ++p)
1931 *   {
1932 *   Tensor<1, dim + 2> solution;
1933 *   for (unsigned int d = 0; d < dim + 2; ++d)
1934 *   solution[d] = inputs.solution_values[p](d);
1935 *  
1936 *   const double density = solution[0];
1937 *   const Tensor<1, dim> velocity = euler_velocity<dim>(solution);
1938 *   const double pressure = euler_pressure<dim>(solution);
1939 *  
1940 *   for (unsigned int d = 0; d < dim; ++d)
1942 *   computed_quantities[p](dim) = pressure;
1943 *   computed_quantities[p](dim + 1) = std::sqrt(gamma * pressure / density);
1944 *  
1945 *   if (do_schlieren_plot == true)
1946 *   computed_quantities[p](dim + 2) =
1947 *   inputs.solution_gradients[p][0] * inputs.solution_gradients[p][0];
1948 *   }
1949 *   }
1950 *  
1951 *  
1952 *  
1953 *   template <int dim>
1954 *   std::vector<std::string> EulerProblem<dim>::Postprocessor::get_names() const
1955 *   {
1956 *   std::vector<std::string> names;
1957 *   for (unsigned int d = 0; d < dim; ++d)
1958 *   names.emplace_back("velocity");
1959 *   names.emplace_back("pressure");
1960 *   names.emplace_back("speed_of_sound");
1961 *  
1962 *   if (do_schlieren_plot == true)
1963 *   names.emplace_back("schlieren_plot");
1964 *  
1965 *   return names;
1966 *   }
1967 *  
1968 *  
1969 *  
1970 *   template <int dim>
1971 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1973 *   {
1974 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1976 *   for (unsigned int d = 0; d < dim; ++d)
1977 *   interpretation.push_back(
1981 *  
1982 *   if (do_schlieren_plot == true)
1983 *   interpretation.push_back(
1985 *  
1986 *   return interpretation;
1987 *   }
1988 *  
1989 *  
1990 *  
1991 *   template <int dim>
1993 *   {
1994 *   if (do_schlieren_plot == true)
1996 *   else
1997 *   return update_values;
1998 *   }
1999 *  
2000 *  
2001 *  
2002 *   template <int dim>
2007 *   #endif
2008 *   , fe(FE_DGQ<dim>(fe_degree), dim + 2)
2009 *   , mapping(fe_degree)
2010 *   , dof_handler(triangulation)
2012 *   , euler_operator(timer)
2013 *   , time(0)
2014 *   , time_step(0)
2015 *   {}
2016 *  
2017 *  
2018 *  
2019 *   template <int dim>
2021 *   {
2022 *   switch (testcase)
2023 *   {
2024 *   case 0:
2025 *   {
2027 *   for (unsigned int d = 1; d < dim; ++d)
2028 *   lower_left[d] = -5;
2029 *  
2031 *   upper_right[0] = 10;
2032 *   for (unsigned int d = 1; d < dim; ++d)
2033 *   upper_right[d] = 5;
2034 *  
2036 *   lower_left,
2037 *   upper_right);
2038 *   triangulation.refine_global(2);
2039 *  
2040 *   euler_operator.set_inflow_boundary(
2041 *   0, std::make_unique<ExactSolution<dim>>(0));
2042 *  
2043 *   break;
2044 *   }
2045 *  
2046 *   case 1:
2047 *   {
2049 *   triangulation, 0.03, 1, 0, true);
2050 *  
2051 *   euler_operator.set_inflow_boundary(
2052 *   0, std::make_unique<ExactSolution<dim>>(0));
2053 *   euler_operator.set_subsonic_outflow_boundary(
2054 *   1, std::make_unique<ExactSolution<dim>>(0));
2055 *  
2056 *   euler_operator.set_wall_boundary(2);
2057 *   euler_operator.set_wall_boundary(3);
2058 *  
2059 *   if (dim == 3)
2060 *   euler_operator.set_body_force(
2061 *   std::make_unique<Functions::ConstantFunction<dim>>(
2062 *   std::vector<double>({0., 0., -0.2})));
2063 *  
2064 *   break;
2065 *   }
2066 *  
2067 *   default:
2069 *   }
2070 *  
2071 *   triangulation.refine_global(n_global_refinements);
2072 *  
2073 *   dof_handler.distribute_dofs(fe);
2074 *  
2075 *   euler_operator.reinit(mapping, dof_handler);
2076 *   euler_operator.initialize_vector(solution);
2077 *  
2078 *   std::locale s = pcout.get_stream().getloc();
2079 *   pcout.get_stream().imbue(std::locale(""));
2080 *   pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
2081 *   << " ( = " << (dim + 2) << " [vars] x "
2082 *   << triangulation.n_global_active_cells() << " [cells] x "
2083 *   << Utilities::pow(fe_degree + 1, dim) << " [dofs/cell/var] )"
2084 *   << std::endl;
2085 *   pcout.get_stream().imbue(s);
2086 *   }
2087 *  
2088 *  
2089 *  
2090 *   template <int dim>
2091 *   void EulerProblem<dim>::output_results(const unsigned int result_number)
2092 *   {
2093 *   const std::array<double, 3> errors =
2094 *   euler_operator.compute_errors(ExactSolution<dim>(time), solution);
2095 *   const std::string quantity_name = testcase == 0 ? "error" : "norm";
2096 *  
2097 *   pcout << "Time:" << std::setw(8) << std::setprecision(3) << time
2098 *   << ", dt: " << std::setw(8) << std::setprecision(2) << time_step
2099 *   << ", " << quantity_name << " rho: " << std::setprecision(4)
2100 *   << std::setw(10) << errors[0] << ", rho * u: " << std::setprecision(4)
2101 *   << std::setw(10) << errors[1] << ", energy:" << std::setprecision(4)
2102 *   << std::setw(10) << errors[2] << std::endl;
2103 *  
2104 *   {
2105 *   TimerOutput::Scope t(timer, "output");
2106 *  
2107 *   Postprocessor postprocessor;
2108 *   DataOut<dim> data_out;
2109 *  
2111 *   flags.write_higher_order_cells = true;
2112 *   data_out.set_flags(flags);
2113 *  
2114 *   data_out.attach_dof_handler(dof_handler);
2115 *   {
2116 *   std::vector<std::string> names;
2117 *   names.emplace_back("density");
2118 *   for (unsigned int d = 0; d < dim; ++d)
2119 *   names.emplace_back("momentum");
2120 *   names.emplace_back("energy");
2121 *  
2122 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
2124 *   interpretation.push_back(
2126 *   for (unsigned int d = 0; d < dim; ++d)
2127 *   interpretation.push_back(
2129 *   interpretation.push_back(
2131 *  
2132 *   data_out.add_data_vector(dof_handler, solution, names, interpretation);
2133 *   }
2134 *   data_out.add_data_vector(solution, postprocessor);
2135 *  
2137 *   if (testcase == 0 && dim == 2)
2138 *   {
2139 *   reference.reinit(solution);
2140 *   euler_operator.project(ExactSolution<dim>(time), reference);
2141 *   reference.sadd(-1., 1, solution);
2142 *   std::vector<std::string> names;
2143 *   names.emplace_back("error_density");
2144 *   for (unsigned int d = 0; d < dim; ++d)
2145 *   names.emplace_back("error_momentum");
2146 *   names.emplace_back("error_energy");
2147 *  
2148 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
2150 *   interpretation.push_back(
2152 *   for (unsigned int d = 0; d < dim; ++d)
2153 *   interpretation.push_back(
2155 *   interpretation.push_back(
2157 *  
2158 *   data_out.add_data_vector(dof_handler,
2159 *   reference,
2160 *   names,
2162 *   }
2163 *  
2164 *   Vector<double> mpi_owner(triangulation.n_active_cells());
2166 *   data_out.add_data_vector(mpi_owner, "owner");
2167 *  
2168 *   data_out.build_patches(mapping,
2169 *   fe.degree,
2171 *  
2172 *   const std::string filename =
2173 *   "solution_" + Utilities::int_to_string(result_number, 3) + ".vtu";
2174 *   data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD);
2175 *   }
2176 *   }
2177 *  
2178 *  
2179 *  
2180 *   template <int dim>
2182 *   {
2183 *   {
2184 *   const unsigned int n_vect_number = VectorizedArrayType::size();
2185 *   const unsigned int n_vect_bits = 8 * sizeof(Number) * n_vect_number;
2186 *  
2187 *   pcout << "Running with "
2189 *   << " MPI processes" << std::endl;
2190 *   pcout << "Vectorization over " << n_vect_number << ' '
2191 *   << (std::is_same_v<Number, double> ? "doubles" : "floats") << " = "
2192 *   << n_vect_bits << " bits ("
2194 *   << std::endl;
2195 *   }
2196 *  
2198 *  
2200 *  
2203 *   rk_register_1.reinit(solution);
2204 *   rk_register_2.reinit(solution);
2205 *  
2206 *   euler_operator.project(ExactSolution<dim>(time), solution);
2207 *  
2208 *  
2209 *   double min_vertex_distance = std::numeric_limits<double>::max();
2210 *   for (const auto &cell : triangulation.active_cell_iterators())
2211 *   if (cell->is_locally_owned())
2213 *   std::min(min_vertex_distance, cell->minimum_vertex_distance());
2216 *  
2217 *   time_step = courant_number * integrator.n_stages() /
2218 *   euler_operator.compute_cell_transport_speed(solution);
2219 *   pcout << "Time step size: " << time_step
2220 *   << ", minimal h: " << min_vertex_distance
2221 *   << ", initial transport scaling: "
2222 *   << 1. / euler_operator.compute_cell_transport_speed(solution)
2223 *   << std::endl
2224 *   << std::endl;
2225 *  
2226 *   output_results(0);
2227 *  
2228 *   unsigned int timestep_number = 0;
2229 *  
2230 *   while (time < final_time - 1e-12 && timestep_number < max_time_steps)
2231 *   {
2233 *   if (timestep_number % 5 == 0)
2234 *   time_step =
2235 *   courant_number * integrator.n_stages() /
2237 *   euler_operator.compute_cell_transport_speed(solution), 3);
2238 *  
2239 *   {
2240 *   TimerOutput::Scope t(timer, "rk time stepping total");
2241 *   integrator.perform_time_step(euler_operator,
2242 *   time,
2243 *   time_step,
2244 *   solution,
2247 *   }
2248 *  
2249 *   time += time_step;
2250 *  
2251 *   if (static_cast<int>(time / output_tick) !=
2252 *   static_cast<int>((time - time_step) / output_tick) ||
2253 *   time >= final_time - 1e-12)
2255 *   static_cast<unsigned int>(std::round(time / output_tick)));
2256 *   }
2257 *  
2258 *   timer.print_wall_time_statistics(MPI_COMM_WORLD);
2259 *   pcout << std::endl;
2260 *   }
2261 *  
2262 *   } // namespace Euler_DG
2263 *  
2264 *  
2265 *   int main(int argc, char **argv)
2266 *   {
2267 *   using namespace Euler_DG;
2268 *   using namespace dealii;
2269 *  
2271 *  
2272 *   try
2273 *   {
2275 *   euler_problem.run();
2276 *   }
2277 *   catch (std::exception &exc)
2278 *   {
2279 *   std::cerr << std::endl
2280 *   << std::endl
2281 *   << "----------------------------------------------------"
2282 *   << std::endl;
2283 *   std::cerr << "Exception on processing: " << std::endl
2284 *   << exc.what() << std::endl
2285 *   << "Aborting!" << std::endl
2286 *   << "----------------------------------------------------"
2287 *   << std::endl;
2288 *  
2289 *   return 1;
2290 *   }
2291 *   catch (...)
2292 *   {
2293 *   std::cerr << std::endl
2294 *   << std::endl
2295 *   << "----------------------------------------------------"
2296 *   << std::endl;
2297 *   std::cerr << "Unknown exception!" << std::endl
2298 *   << "Aborting!" << std::endl
2299 *   << "----------------------------------------------------"
2300 *   << std::endl;
2301 *   return 1;
2302 *   }
2303 *  
2304 *   return 0;
2305 *   }
2306 * @endcode
2307<a name="step_76-Results"></a><h1>Results</h1>
2308
2309
2310Running the program with the default settings on a machine with 40 processes
2311produces the following output:
2312
2313@code
2316Number of degrees of freedom: 27.648.000 ( = 5 [vars] x 25.600 [cells] x 216 [dofs/cell/var] )
2317Time step size: 0.000295952, minimal h: 0.0075, initial transport scaling: 0.00441179
2318Time: 0, dt: 0.0003, norm rho: 5.385e-16, rho * u: 1.916e-16, energy: 1.547e-15
2319+--------------------------------------+------------------+------------+------------------+
2320| Total wallclock time elapsed | 17.52s 10 | 17.52s | 17.52s 11 |
2321| | | |
2322| Section | no. calls | min time rank | avg time | max time rank |
2323+--------------------------------------+------------------+------------+------------------+
2324| compute errors | 1 | 0.009594s 16 | 0.009705s | 0.009819s 8 |
2325| compute transport speed | 22 | 0.1366s 0 | 0.1367s | 0.1368s 18 |
2326| output | 1 | 1.233s 0 | 1.233s | 1.233s 32 |
2327| rk time stepping total | 100 | 8.746s 35 | 8.746s | 8.746s 0 |
2328| rk_stage - integrals L_h | 500 | 8.742s 36 | 8.742s | 8.743s 2 |
2329+--------------------------------------+------------------+------------+------------------+
2330@endcode
2331
2332and the following visual output:
2333
2334<table align="center" class="doxtable" style="width:85%">
2335 <tr>
2336 <td>
2337 <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_010.png" alt="" width="100%">
2338 </td>
2339 <td>
2340 <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_025.png" alt="" width="100%">
2341 </td>
2342 </tr>
2343 <tr>
2344 <td>
2345 <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_050.png" alt="" width="100%">
2346 </td>
2347 <td>
2348 <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_100.png" alt="" width="100%">
2349 </td>
2350 </tr>
2351</table>
2352
2353As a reference, the results of @ref step_67 "step-67" using FCL are:
2354
2355@code
2358Number of degrees of freedom: 27.648.000 ( = 5 [vars] x 25.600 [cells] x 216 [dofs/cell/var] )
2359Time step size: 0.000295952, minimal h: 0.0075, initial transport scaling: 0.00441179
2360Time: 0, dt: 0.0003, norm rho: 5.385e-16, rho * u: 1.916e-16, energy: 1.547e-15
2361+-------------------------------------------+------------------+------------+------------------+
2362| Total wallclock time elapsed | 13.33s 0 | 13.34s | 13.35s 34 |
2363| | | |
2364| Section | no. calls | min time rank | avg time | max time rank |
2365+-------------------------------------------+------------------+------------+------------------+
2366| compute errors | 1 | 0.007977s 10 | 0.008053s | 0.008161s 30 |
2367| compute transport speed | 22 | 0.1228s 34 | 0.2227s | 0.3845s 0 |
2368| output | 1 | 1.255s 3 | 1.257s | 1.259s 27 |
2369| rk time stepping total | 100 | 11.15s 0 | 11.32s | 11.42s 34 |
2370| rk_stage - integrals L_h | 500 | 8.719s 10 | 8.932s | 9.196s 0 |
2371| rk_stage - inv mass + vec upd | 500 | 1.944s 0 | 2.377s | 2.55s 10 |
2372+-------------------------------------------+------------------+------------+------------------+
2373@endcode
2374
2377
2378<a name="step_76-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
2379
2380
2382<a href="https://github.com/hyperdeal/hyperdeal/blob/a9e67b4e625ff1dde2fed93ad91cdfacfaa3acdf/include/hyper.deal/operators/advection/advection_operation.h#L219-L569">advection operator based on cell-centric loops</a>
2385
2386<a name="step_76-ExtensiontothecompressibleNavierStokesequations"></a><h4>Extension to the compressible Navier-Stokes equations</h4>
2387
2388
2392the additional cost of elliptic terms, e.g. via an interior penalty method, that
2394the @ref step_59 "step-59" tutorial program. The reasoning behind this switch is that in the
2395case of FE_DGQ all values of neighboring cells (i.e., @f$k+1@f$ layers) are needed,
2402
2403<a name="step_76-BlockGaussSeidellikepreconditioners"></a><h4>Block Gauss-Seidel-like preconditioners</h4>
2404
2405
2406Cell-centric loops could be used to create block Gauss-Seidel preconditioners
2411
2412@code
2413// vector monitor if cells have been updated or not
2414Vector<Number> visit_flags(data.n_cell_batches () + data.n_ghost_cell_batches ());
2415
2416// element centric loop with a modified kernel
2418 [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
2419
2420 for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2421 {
2422 // cell integral as usual (not shown)
2423
2424 for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face)
2425 {
2426 const auto boundary_id = data.get_faces_by_cells_boundary_id(cell, face)[0];
2427
2428 if (boundary_id == numbers::internal_face_boundary_id)
2429 {
2430 phi_p.reinit(cell, face);
2431
2432 const auto flags = phi_p.read_cell_data(visit_flags);
2434 std::min(flags.begin(),
2435 flags().begin() + data.n_active_entries_per_cell_batch(cell) == 1;
2436
2438 phi_p.gather_evaluate(dst, EvaluationFlags::values);
2439 else
2440 phi_p.gather_evaluate(src, EvaluationFlags::values);
2441
2442 // continue as usual (not shown)
2443 }
2444 else
2445 {
2446 // boundary integral as usual (not shown)
2447 }
2448 }
2449
2450 // continue as above and apply your favorite algorithm to invert
2451 // the cell-local operator (not shown)
2452
2453 // make cells as updated
2454 phi.set_cell_data(visit_flags, VectorizedArrayType(1.0));
2455 }
2456 },
2457 dst,
2458 src,
2459 true,
2461@endcode
2462
2463For this purpose, one can exploit the cell-data vector capabilities of
2465
2466Please note that in the given example we process <code>VectorizedArrayType::size()</code>
2467number of blocks, since each lane corresponds to one block. We consider blocks
2468as updated if all blocks processed by a vector register have been updated. In
2471efficiency of the preconditioner. A reduction of cells processed in parallel
2472by explicitly reducing the number of lanes used by <code>VectorizedArrayType</code>
2475"possibility for extension": vectorization within an element.
2476 *
2477 *
2478<a name="step_76-PlainProg"></a>
2479<h1> The plain program</h1>
2480@include "step-76.cc"
2481*/
void submit_value(const value_type val_in, const unsigned int q_point)
Abstract base class for mapping classes.
Definition mapping.h:320
void loop_cell_centric(void(CLASS::*cell_operation)(const MatrixFree &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &) const, const CLASS *owning_class, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
friend class Tensor
Definition tensor.h:865
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
@ wall_times
Definition timer.h:651
std::enable_if_t< std::is_floating_point_v< T > &&std::is_floating_point_v< U > &&!std::is_same_v< T, U >, typename ProductType< std::complex< T >, std::complex< U > >::type > operator*(const std::complex< T > &left, const std::complex< U > &right)
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
unsigned int level
Definition grid_out.cc:4632
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
UpdateFlags
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
const Event initial
Definition event.cc:68
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
@ matrix
Contents is actually a matrix.
@ general
No special properties.
constexpr char A
constexpr types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition divergence.h:471
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
@ LOW_STORAGE_RK_STAGE9_ORDER5
@ LOW_STORAGE_RK_STAGE3_ORDER3
@ LOW_STORAGE_RK_STAGE7_ORDER4
@ LOW_STORAGE_RK_STAGE5_ORDER4
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:97
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:112
T reduce(const T &local_value, const MPI_Comm comm, const std::function< T(const T &, const T &)> &combiner, const unsigned int root_process=0)
std::string get_time()
std::string get_current_vectorization_level()
Definition utilities.cc:942
Number truncate_to_n_digits(const Number number, const unsigned int n_digits)
Definition utilities.cc:582
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:474
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >()), const bool project_to_boundary_first=false)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
constexpr double PI
Definition numbers.h:241
constexpr unsigned int invalid_unsigned_int
Definition types.h:232
constexpr types::boundary_id internal_face_boundary_id
Definition types.h:323
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Function &function, const unsigned int grainsize)
Definition parallel.h:204
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition types.h:32
unsigned int boundary_id
Definition types.h:153
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation