Reference documentation for deal.II version 9.5.0
\(\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
NavierStokes_TRBDF2_DG.h
Go to the documentation of this file.
1
183 *  
184 * @endcode
185 *
186 * We declare class that describes the boundary conditions and initial one for velocity:
187 *
188
189 *
190 *
191 * @code
192 *   template<int dim>
193 *   class Velocity: public Function<dim> {
194 *   public:
195 *   Velocity(const double initial_time = 0.0);
196 *  
197 *   virtual double value(const Point<dim>& p,
198 *   const unsigned int component = 0) const override;
199 *  
200 *   virtual void vector_value(const Point<dim>& p,
201 *   Vector<double>& values) const override;
202 *   };
203 *  
204 *  
205 *   template<int dim>
206 *   Velocity<dim>::Velocity(const double initial_time): Function<dim>(dim, initial_time) {}
207 *  
208 *  
209 *   template<int dim>
210 *   double Velocity<dim>::value(const Point<dim>& p, const unsigned int component) const {
211 *   AssertIndexRange(component, 3);
212 *   if(component == 0) {
213 *   const double Um = 1.5;
214 *   const double H = 4.1;
215 *  
216 *   return 4.0*Um*p(1)*(H - p(1))/(H*H);
217 *   }
218 *   else
219 *   return 0.0;
220 *   }
221 *  
222 *  
223 *   template<int dim>
224 *   void Velocity<dim>::vector_value(const Point<dim>& p, Vector<double>& values) const {
225 *   Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
226 *  
227 *   for(unsigned int i = 0; i < dim; ++i)
228 *   values[i] = value(p, i);
229 *   }
230 *  
231 *  
232 * @endcode
233 *
234 * We do the same for the pressure
235 *
236
237 *
238 *
239 * @code
240 *   template<int dim>
241 *   class Pressure: public Function<dim> {
242 *   public:
243 *   Pressure(const double initial_time = 0.0);
244 *  
245 *   virtual double value(const Point<dim>& p,
246 *   const unsigned int component = 0) const override;
247 *   };
248 *  
249 *  
250 *   template<int dim>
251 *   Pressure<dim>::Pressure(const double initial_time): Function<dim>(1, initial_time) {}
252 *  
253 *  
254 *   template<int dim>
255 *   double Pressure<dim>::value(const Point<dim>& p, const unsigned int component) const {
256 *   (void)component;
257 *   AssertIndexRange(component, 1);
258 *  
259 *   return 22.0 - p(0);
260 *   }
261 *  
262 *   } // namespace EquationData
263 * @endcode
264
265
266<a name="ann-navier_stokes_TRBDF2_DG.cc"></a>
267<h1>Annotated version of navier_stokes_TRBDF2_DG.cc</h1>
268 *
269 *
270 *
271 *
272 * @code
273 *   /* Author: Giuseppe Orlando, 2022. */
274 *  
275 * @endcode
276 *
277 * We start by including all the necessary deal.II header files and some C++
278 * related ones.
279 *
280
281 *
282 *
283 * @code
284 *   #include <deal.II/base/quadrature_lib.h>
285 *   #include <deal.II/base/multithread_info.h>
286 *   #include <deal.II/base/thread_management.h>
287 *   #include <deal.II/base/work_stream.h>
288 *   #include <deal.II/base/parallel.h>
289 *   #include <deal.II/base/utilities.h>
290 *   #include <deal.II/base/conditional_ostream.h>
291 *  
292 *   #include <deal.II/lac/vector.h>
293 *   #include <deal.II/lac/solver_cg.h>
294 *   #include <deal.II/lac/precondition.h>
295 *   #include <deal.II/lac/solver_gmres.h>
296 *   #include <deal.II/lac/affine_constraints.h>
297 *  
298 *   #include <deal.II/grid/tria.h>
299 *   #include <deal.II/grid/grid_generator.h>
300 *   #include <deal.II/grid/grid_tools.h>
301 *   #include <deal.II/grid/grid_refinement.h>
302 *   #include <deal.II/grid/tria_accessor.h>
303 *   #include <deal.II/grid/tria_iterator.h>
304 *   #include <deal.II/distributed/grid_refinement.h>
305 *  
306 *   #include <deal.II/dofs/dof_handler.h>
307 *   #include <deal.II/dofs/dof_accessor.h>
308 *   #include <deal.II/dofs/dof_tools.h>
309 *  
310 *   #include <deal.II/fe/fe_q.h>
311 *   #include <deal.II/fe/fe_dgq.h>
312 *   #include <deal.II/fe/fe_values.h>
313 *   #include <deal.II/fe/fe_tools.h>
314 *   #include <deal.II/fe/fe_system.h>
315 *  
316 *   #include <deal.II/numerics/matrix_tools.h>
317 *   #include <deal.II/numerics/vector_tools.h>
318 *   #include <deal.II/numerics/data_out.h>
319 *  
320 *   #include <fstream>
321 *   #include <cmath>
322 *   #include <iostream>
323 *  
324 *   #include <deal.II/matrix_free/matrix_free.h>
325 *   #include <deal.II/matrix_free/operators.h>
326 *   #include <deal.II/matrix_free/fe_evaluation.h>
327 *   #include <deal.II/fe/component_mask.h>
328 *  
329 *   #include <deal.II/base/timer.h>
330 *   #include <deal.II/distributed/solution_transfer.h>
331 *   #include <deal.II/numerics/error_estimator.h>
332 *  
333 *   #include <deal.II/multigrid/multigrid.h>
334 *   #include <deal.II/multigrid/mg_transfer_matrix_free.h>
335 *   #include <deal.II/multigrid/mg_tools.h>
336 *   #include <deal.II/multigrid/mg_coarse.h>
337 *   #include <deal.II/multigrid/mg_smoother.h>
338 *   #include <deal.II/multigrid/mg_matrix.h>
339 *  
340 *   #include <deal.II/meshworker/mesh_loop.h>
341 *  
342 *   #include "runtime_parameters.h"
343 *   #include "equation_data.h"
344 *  
345 *   namespace MatrixFreeTools {
346 *   using namespace dealii;
347 *  
348 *   template<int dim, typename Number, typename VectorizedArrayType>
350 *   LinearAlgebra::distributed::Vector<Number>& diagonal_global,
351 *   const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType>&,
353 *   const unsigned int&,
354 *   const std::pair<unsigned int, unsigned int>&)>& cell_operation,
355 *   const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType>&,
357 *   const unsigned int&,
358 *   const std::pair<unsigned int, unsigned int>&)>& face_operation,
359 *   const std::function<void(const MatrixFree<dim, Number, VectorizedArrayType>&,
361 *   const unsigned int&,
362 *   const std::pair<unsigned int, unsigned int>&)>& boundary_operation,
363 *   const unsigned int dof_no = 0) {
364 * @endcode
365 *
366 * initialize vector
367 *
368 * @code
369 *   matrix_free.initialize_dof_vector(diagonal_global, dof_no);
370 *  
371 *   const unsigned int dummy = 0;
372 *  
373 *   matrix_free.loop(cell_operation, face_operation, boundary_operation,
374 *   diagonal_global, dummy, false,
377 *   }
378 *   }
379 *  
380 * @endcode
381 *
382 * We include the code in a suitable namespace:
383 *
384
385 *
386 *
387 * @code
388 *   namespace NS_TRBDF2 {
389 *   using namespace dealii;
390 *  
391 * @endcode
392 *
393 * The following class is an auxiliary one for post-processing of the vorticity
394 *
395
396 *
397 *
398 * @code
399 *   template<int dim>
400 *   class PostprocessorVorticity: public DataPostprocessor<dim> {
401 *   public:
402 *   virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector<dim>& inputs,
403 *   std::vector<Vector<double>>& computed_quantities) const override;
404 *  
405 *   virtual std::vector<std::string> get_names() const override;
406 *  
407 *   virtual std::vector<DataComponentInterpretation::DataComponentInterpretation>
408 *   get_data_component_interpretation() const override;
409 *  
410 *   virtual UpdateFlags get_needed_update_flags() const override;
411 *   };
412 *  
413 * @endcode
414 *
415 * This function evaluates the vorticty in both 2D and 3D cases
416 *
417
418 *
419 *
420 * @code
421 *   template <int dim>
422 *   void PostprocessorVorticity<dim>::evaluate_vector_field(const DataPostprocessorInputs::Vector<dim>& inputs,
423 *   std::vector<Vector<double>>& computed_quantities) const {
424 *   const unsigned int n_quadrature_points = inputs.solution_values.size();
425 *  
426 *   /*--- Check the correctness of all data structres ---*/
427 *   Assert(inputs.solution_gradients.size() == n_quadrature_points, ExcInternalError());
428 *   Assert(computed_quantities.size() == n_quadrature_points, ExcInternalError());
429 *  
430 *   Assert(inputs.solution_values[0].size() == dim, ExcInternalError());
431 *  
432 *   if(dim == 2) {
433 *   Assert(computed_quantities[0].size() == 1, ExcInternalError());
434 *   }
435 *   else {
436 *   Assert(computed_quantities[0].size() == dim, ExcInternalError());
437 *   }
438 *  
439 *   /*--- Compute the vorticty ---*/
440 *   if(dim == 2) {
441 *   for(unsigned int q = 0; q < n_quadrature_points; ++q)
442 *   computed_quantities[q](0) = inputs.solution_gradients[q][1][0] - inputs.solution_gradients[q][0][1];
443 *   }
444 *   else {
445 *   for(unsigned int q = 0; q < n_quadrature_points; ++q) {
446 *   computed_quantities[q](0) = inputs.solution_gradients[q][2][1] - inputs.solution_gradients[q][1][2];
447 *   computed_quantities[q](1) = inputs.solution_gradients[q][0][2] - inputs.solution_gradients[q][2][0];
448 *   computed_quantities[q](2) = inputs.solution_gradients[q][1][0] - inputs.solution_gradients[q][0][1];
449 *   }
450 *   }
451 *   }
452 *  
453 * @endcode
454 *
455 * This auxiliary function is required by the base class DataProcessor and simply
456 * sets the name for the output file
457 *
458
459 *
460 *
461 * @code
462 *   template<int dim>
463 *   std::vector<std::string> PostprocessorVorticity<dim>::get_names() const {
464 *   std::vector<std::string> names;
465 *   names.emplace_back("vorticity");
466 *   if(dim == 3) {
467 *   names.emplace_back("vorticity");
468 *   names.emplace_back("vorticity");
469 *   }
470 *  
471 *   return names;
472 *   }
473 *  
474 * @endcode
475 *
476 * This auxiliary function is required by the base class DataProcessor and simply
477 * specifies if the vorticity is a scalar (2D) or a vector (3D)
478 *
479
480 *
481 *
482 * @code
483 *   template<int dim>
484 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
485 *   PostprocessorVorticity<dim>::get_data_component_interpretation() const {
486 *   std::vector<DataComponentInterpretation::DataComponentInterpretation> interpretation;
487 *   if(dim == 2)
488 *   interpretation.push_back(DataComponentInterpretation::component_is_scalar);
489 *   else {
493 *   }
494 *  
495 *   return interpretation;
496 *   }
497 *  
498 * @endcode
499 *
500 * This auxiliary function is required by the base class DataProcessor and simply
501 * sets which variables have to updated (only the gradients)
502 *
503
504 *
505 *
506 * @code
507 *   template<int dim>
508 *   UpdateFlags PostprocessorVorticity<dim>::get_needed_update_flags() const {
509 *   return update_gradients;
510 *   }
511 *  
512 *  
513 * @endcode
514 *
515 * The following structs are auxiliary objects for mesh refinement. ScratchData simply sets
516 * the FEValues object
517 *
518
519 *
520 *
521 * @code
522 *   template <int dim>
523 *   struct ScratchData {
524 *   ScratchData(const FiniteElement<dim>& fe,
525 *   const unsigned int quadrature_degree,
526 *   const UpdateFlags update_flags): fe_values(fe, QGauss<dim>(quadrature_degree), update_flags) {}
527 *  
528 *   ScratchData(const ScratchData<dim>& scratch_data): fe_values(scratch_data.fe_values.get_fe(),
529 *   scratch_data.fe_values.get_quadrature(),
530 *   scratch_data.fe_values.get_update_flags()) {}
531 *   FEValues<dim> fe_values;
532 *   };
533 *  
534 *  
535 * @endcode
536 *
537 * CopyData simply sets the cell index
538 *
539
540 *
541 *
542 * @code
543 *   struct CopyData {
544 *   CopyData() : cell_index(numbers::invalid_unsigned_int), value(0.0) {}
545 *  
546 *   CopyData(const CopyData &) = default;
547 *  
548 *   unsigned int cell_index;
549 *   double value;
550 *   };
551 *  
552 *  
553 * @endcode
554 *
555 *
556 * <a name=""></a>
557 * @sect{ <code>NavierStokesProjectionOperator::NavierStokesProjectionOperator</code> }
558 *
559
560 *
561 * The following class sets effecively the weak formulation of the problems for the different stages
562 * and for both velocity and pressure.
563 * The template parameters are the dimnesion of the problem, the polynomial degree for the pressure,
564 * the polynomial degree for the velocity, the number of quadrature points for integrals for the pressure step,
565 * the number of quadrature points for integrals for the velocity step, the type of vector for storage and the type
566 * of floating point data (in general double or float for preconditioners structures if desired).
567 *
568
569 *
570 *
571 * @code
572 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
573 *   class NavierStokesProjectionOperator: public MatrixFreeOperators::Base<dim, Vec> {
574 *   public:
575 *   using Number = typename Vec::value_type;
576 *  
577 *   NavierStokesProjectionOperator();
578 *  
579 *   NavierStokesProjectionOperator(RunTimeParameters::Data_Storage& data);
580 *  
581 *   void set_dt(const double time_step);
582 *  
583 *   void set_TR_BDF2_stage(const unsigned int stage);
584 *  
585 *   void set_NS_stage(const unsigned int stage);
586 *  
587 *   void set_u_extr(const Vec& src);
588 *  
589 *   void vmult_rhs_velocity(Vec& dst, const std::vector<Vec>& src) const;
590 *  
591 *   void vmult_rhs_pressure(Vec& dst, const std::vector<Vec>& src) const;
592 *  
593 *   void vmult_grad_p_projection(Vec& dst, const Vec& src) const;
594 *  
595 *   virtual void compute_diagonal() override;
596 *  
597 *   protected:
598 *   double Re;
599 *   double dt;
600 *  
601 *   /*--- Parameters of time-marching scheme ---*/
602 *   double gamma;
603 *   double a31;
604 *   double a32;
605 *   double a33;
606 *  
607 *   unsigned int TR_BDF2_stage; /*--- Flag to denote at which stage of the TR-BDF2 are ---*/
608 *   unsigned int NS_stage; /*--- Flag to denote at which stage of NS solution inside each TR-BDF2 stage we are
609 *   (solution of the velocity or of the pressure)---*/
610 *  
611 *   virtual void apply_add(Vec& dst, const Vec& src) const override;
612 *  
613 *   private:
614 *   /*--- Auxiliary variable for the TR stage
615 *   (just to avoid to report a lot of 0.5 and for my personal choice to be coherent with the article) ---*/
616 *   const double a21 = 0.5;
617 *   const double a22 = 0.5;
618 *  
619 *   /*--- Penalty method parameters, theta = 1 means SIP, while C_p and C_u are the penalization coefficients ---*/
620 *   const double theta_v = 1.0;
621 *   const double theta_p = 1.0;
622 *   const double C_p = 1.0*(fe_degree_p + 1)*(fe_degree_p + 1);
623 *   const double C_u = 1.0*(fe_degree_v + 1)*(fe_degree_v + 1);
624 *  
625 *   Vec u_extr; /*--- Auxiliary variable to update the extrapolated velocity ---*/
626 *  
627 *   EquationData::Velocity<dim> vel_boundary_inflow; /*--- Auxiliary variable to impose velocity boundary conditions ---*/
628 *  
629 *   /*--- The following functions basically assemble the linear and bilinear forms. Their syntax is due to
630 *   the base class MatrixFreeOperators::Base ---*/
631 *   void assemble_rhs_cell_term_velocity(const MatrixFree<dim, Number>& data,
632 *   Vec& dst,
633 *   const std::vector<Vec>& src,
634 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
635 *   void assemble_rhs_face_term_velocity(const MatrixFree<dim, Number>& data,
636 *   Vec& dst,
637 *   const std::vector<Vec>& src,
638 *   const std::pair<unsigned int, unsigned int>& face_range) const;
639 *   void assemble_rhs_boundary_term_velocity(const MatrixFree<dim, Number>& data,
640 *   Vec& dst,
641 *   const std::vector<Vec>& src,
642 *   const std::pair<unsigned int, unsigned int>& face_range) const;
643 *  
644 *   void assemble_rhs_cell_term_pressure(const MatrixFree<dim, Number>& data,
645 *   Vec& dst,
646 *   const std::vector<Vec>& src,
647 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
648 *   void assemble_rhs_face_term_pressure(const MatrixFree<dim, Number>& data,
649 *   Vec& dst,
650 *   const std::vector<Vec>& src,
651 *   const std::pair<unsigned int, unsigned int>& face_range) const;
652 *   void assemble_rhs_boundary_term_pressure(const MatrixFree<dim, Number>& data,
653 *   Vec& dst,
654 *   const std::vector<Vec>& src,
655 *   const std::pair<unsigned int, unsigned int>& face_range) const;
656 *  
657 *   void assemble_cell_term_velocity(const MatrixFree<dim, Number>& data,
658 *   Vec& dst,
659 *   const Vec& src,
660 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
661 *   void assemble_face_term_velocity(const MatrixFree<dim, Number>& data,
662 *   Vec& dst,
663 *   const Vec& src,
664 *   const std::pair<unsigned int, unsigned int>& face_range) const;
665 *   void assemble_boundary_term_velocity(const MatrixFree<dim, Number>& data,
666 *   Vec& dst,
667 *   const Vec& src,
668 *   const std::pair<unsigned int, unsigned int>& face_range) const;
669 *  
670 *   void assemble_cell_term_pressure(const MatrixFree<dim, Number>& data,
671 *   Vec& dst,
672 *   const Vec& src,
673 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
674 *   void assemble_face_term_pressure(const MatrixFree<dim, Number>& data,
675 *   Vec& dst,
676 *   const Vec& src,
677 *   const std::pair<unsigned int, unsigned int>& face_range) const;
678 *   void assemble_boundary_term_pressure(const MatrixFree<dim, Number>& data,
679 *   Vec& dst,
680 *   const Vec& src,
681 *   const std::pair<unsigned int, unsigned int>& face_range) const;
682 *  
683 *   void assemble_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
684 *   Vec& dst,
685 *   const Vec& src,
686 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
687 *   void assemble_rhs_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
688 *   Vec& dst,
689 *   const Vec& src,
690 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
691 *  
692 *   void assemble_diagonal_cell_term_velocity(const MatrixFree<dim, Number>& data,
693 *   Vec& dst,
694 *   const unsigned int& src,
695 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
696 *   void assemble_diagonal_face_term_velocity(const MatrixFree<dim, Number>& data,
697 *   Vec& dst,
698 *   const unsigned int& src,
699 *   const std::pair<unsigned int, unsigned int>& face_range) const;
700 *   void assemble_diagonal_boundary_term_velocity(const MatrixFree<dim, Number>& data,
701 *   Vec& dst,
702 *   const unsigned int& src,
703 *   const std::pair<unsigned int, unsigned int>& face_range) const;
704 *  
705 *   void assemble_diagonal_cell_term_pressure(const MatrixFree<dim, Number>& data,
706 *   Vec& dst,
707 *   const unsigned int& src,
708 *   const std::pair<unsigned int, unsigned int>& cell_range) const;
709 *   void assemble_diagonal_face_term_pressure(const MatrixFree<dim, Number>& data,
710 *   Vec& dst,
711 *   const unsigned int& src,
712 *   const std::pair<unsigned int, unsigned int>& face_range) const;
713 *   void assemble_diagonal_boundary_term_pressure(const MatrixFree<dim, Number>& data,
714 *   Vec& dst,
715 *   const unsigned int& src,
716 *   const std::pair<unsigned int, unsigned int>& face_range) const;
717 *   };
718 *  
719 *  
720 * @endcode
721 *
722 * We start with the default constructor. It is important for MultiGrid, so it is fundamental
723 * to properly set the parameters of the time scheme.
724 *
725
726 *
727 *
728 * @code
729 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
730 *   NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
731 *   NavierStokesProjectionOperator():
732 *   MatrixFreeOperators::Base<dim, Vec>(), Re(), dt(), gamma(2.0 - std::sqrt(2.0)), a31((1.0 - gamma)/(2.0*(2.0 - gamma))),
733 *   a32(a31), a33(1.0/(2.0 - gamma)), TR_BDF2_stage(1), NS_stage(1), u_extr() {}
734 *  
735 *  
736 * @endcode
737 *
738 * We focus now on the constructor with runtime parameters storage
739 *
740
741 *
742 *
743 * @code
744 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
745 *   NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
746 *   NavierStokesProjectionOperator(RunTimeParameters::Data_Storage& data):
747 *   MatrixFreeOperators::Base<dim, Vec>(), Re(data.Reynolds), dt(data.dt),
748 *   gamma(2.0 - std::sqrt(2.0)), a31((1.0 - gamma)/(2.0*(2.0 - gamma))),
749 *   a32(a31), a33(1.0/(2.0 - gamma)), TR_BDF2_stage(1), NS_stage(1), u_extr(),
750 *   vel_boundary_inflow(data.initial_time) {}
751 *  
752 *  
753 * @endcode
754 *
755 * Setter of time-step (called by Multigrid and in case a smaller time-step towards the end is needed)
756 *
757
758 *
759 *
760 * @code
761 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
762 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
763 *   set_dt(const double time_step) {
764 *   dt = time_step;
765 *   }
766 *  
767 *  
768 * @endcode
769 *
770 * Setter of TR-BDF2 stage (this can be known only during the effective execution
771 * and so it has to be demanded to the class that really solves the problem)
772 *
773
774 *
775 *
776 * @code
777 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
778 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
779 *   set_TR_BDF2_stage(const unsigned int stage) {
780 *   AssertIndexRange(stage, 3);
781 *   Assert(stage > 0, ExcInternalError());
782 *  
783 *   TR_BDF2_stage = stage;
784 *   }
785 *  
786 *  
787 * @endcode
788 *
789 * Setter of NS stage (this can be known only during the effective execution
790 * and so it has to be demanded to the class that really solves the problem)
791 *
792
793 *
794 *
795 * @code
796 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
797 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
798 *   set_NS_stage(const unsigned int stage) {
799 *   AssertIndexRange(stage, 4);
800 *   Assert(stage > 0, ExcInternalError());
801 *  
802 *   NS_stage = stage;
803 *   }
804 *  
805 *  
806 * @endcode
807 *
808 * Setter of extrapolated velocity for different stages
809 *
810
811 *
812 *
813 * @code
814 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
815 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
816 *   set_u_extr(const Vec& src) {
817 *   u_extr = src;
818 *   u_extr.update_ghost_values();
819 *   }
820 *  
821 *  
822 * @endcode
823 *
824 * We are in a DG-MatrixFree framework, so it is convenient to compute separately cell contribution,
825 * internal faces contributions and boundary faces contributions. We start by
826 * assembling the rhs cell term for the velocity.
827 *
828
829 *
830 *
831 * @code
832 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
833 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
834 *   assemble_rhs_cell_term_velocity(const MatrixFree<dim, Number>& data,
835 *   Vec& dst,
836 *   const std::vector<Vec>& src,
837 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
838 *   if(TR_BDF2_stage == 1) {
839 *   /*--- We first start by declaring the suitable instances to read the old velocity, the
840 *   extrapolated velocity and the old pressure. 'phi' will be used only to submit the result.
841 *   The second argument specifies which dof handler has to be used (in this implementation 0 stands for
842 *   velocity and 1 for pressure). ---*/
844 *   phi_old(data, 0),
845 *   phi_old_extr(data, 0);
847 *  
848 *   /*--- We loop over the cells in the range ---*/
849 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
850 *   /*--- Now we need to assign the current cell to each FEEvaluation object and then to specify which src vector
851 *   it has to read (the proper order is clearly delegated to the user, which has to pay attention in the function
852 *   call to be coherent). ---*/
853 *   phi_old.reinit(cell);
854 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
855 *   /*--- The 'gather_evaluate' function reads data from the vector.
856 *   The second and third parameter specifies if you want to read
857 *   values and/or derivative related quantities ---*/
858 *   phi_old_extr.reinit(cell);
859 *   phi_old_extr.gather_evaluate(src[1], EvaluationFlags::values);
860 *   phi_old_press.reinit(cell);
861 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
862 *   phi.reinit(cell);
863 *  
864 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
865 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
866 *   const auto& u_n = phi_old.get_value(q);
867 *   const auto& grad_u_n = phi_old.get_gradient(q);
868 *   const auto& u_n_gamma_ov_2 = phi_old_extr.get_value(q);
869 *   const auto& tensor_product_u_n = outer_product(u_n, u_n_gamma_ov_2);
870 *   const auto& p_n = phi_old_press.get_value(q);
871 *   auto p_n_times_identity = tensor_product_u_n;
872 *   p_n_times_identity = 0;
873 *   for(unsigned int d = 0; d < dim; ++d)
874 *   p_n_times_identity[d][d] = p_n;
875 *  
876 *   phi.submit_value(1.0/(gamma*dt)*u_n, q); /*--- 'submit_value' contains quantites that we want to test against the
877 *   test function ---*/
878 *   phi.submit_gradient(-a21/Re*grad_u_n + a21*tensor_product_u_n + p_n_times_identity, q);
879 *   /*--- 'submit_gradient' contains quantites that we want to test against the gradient of test function ---*/
880 *   }
881 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
882 *   /*--- 'integrate_scatter' is the responsible of distributing into dst.
883 *   The flag parameter specifies if we are testing against the test function and/or its gradient ---*/
884 *   }
885 *   }
886 *   else {
887 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
889 *   phi_old(data, 0),
890 *   phi_int(data, 0);
892 *  
893 *   /*--- We loop over the cells in the range ---*/
894 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
895 *   phi_old.reinit(cell);
896 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
897 *   phi_int.reinit(cell);
898 *   phi_int.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
899 *   phi_old_press.reinit(cell);
900 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
901 *   phi.reinit(cell);
902 *  
903 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
904 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
905 *   const auto& u_n = phi_old.get_value(q);
906 *   const auto& grad_u_n = phi_old.get_gradient(q);
907 *   const auto& u_n_gamma = phi_int.get_value(q);
908 *   const auto& grad_u_n_gamma = phi_int.get_gradient(q);
909 *   const auto& tensor_product_u_n = outer_product(u_n, u_n);
910 *   const auto& tensor_product_u_n_gamma = outer_product(u_n_gamma, u_n_gamma);
911 *   const auto& p_n = phi_old_press.get_value(q);
912 *   auto p_n_times_identity = tensor_product_u_n;
913 *   p_n_times_identity = 0;
914 *   for(unsigned int d = 0; d < dim; ++d)
915 *   p_n_times_identity[d][d] = p_n;
916 *  
917 *   phi.submit_value(1.0/((1.0 - gamma)*dt)*u_n_gamma, q);
918 *   phi.submit_gradient(a32*tensor_product_u_n_gamma + a31*tensor_product_u_n -
919 *   a32/Re*grad_u_n_gamma - a31/Re*grad_u_n + p_n_times_identity, q);
920 *   }
921 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
922 *   }
923 *   }
924 *   }
925 *  
926 *  
927 * @endcode
928 *
929 * The followinf function assembles rhs face term for the velocity
930 *
931
932 *
933 *
934 * @code
935 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
936 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
937 *   assemble_rhs_face_term_velocity(const MatrixFree<dim, Number>& data,
938 *   Vec& dst,
939 *   const std::vector<Vec>& src,
940 *   const std::pair<unsigned int, unsigned int>& face_range) const {
941 *   if(TR_BDF2_stage == 1) {
942 *   /*--- We first start by declaring the suitable instances to read already available quantities. In this case
943 *   we are at the face between two elements and this is the reason of 'FEFaceEvaluation'. It contains an extra
944 *   input argument, the second one, that specifies if it is from 'interior' or not---*/
946 *   phi_m(data, false, 0),
947 *   phi_old_p(data, true, 0),
948 *   phi_old_m(data, false, 0),
949 *   phi_old_extr_p(data, true, 0),
950 *   phi_old_extr_m(data, false, 0);
951 *   FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number> phi_old_press_p(data, true, 1),
952 *   phi_old_press_m(data, false, 1);
953 *  
954 *   /*--- We loop over the faces in the range ---*/
955 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
956 *   phi_old_p.reinit(face);
957 *   phi_old_p.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
958 *   phi_old_m.reinit(face);
959 *   phi_old_m.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
960 *   phi_old_extr_p.reinit(face);
961 *   phi_old_extr_p.gather_evaluate(src[1], EvaluationFlags::values);
962 *   phi_old_extr_m.reinit(face);
963 *   phi_old_extr_m.gather_evaluate(src[1], EvaluationFlags::values);
964 *   phi_old_press_p.reinit(face);
965 *   phi_old_press_p.gather_evaluate(src[2], EvaluationFlags::values);
966 *   phi_old_press_m.reinit(face);
967 *   phi_old_press_m.gather_evaluate(src[2], EvaluationFlags::values);
968 *   phi_p.reinit(face);
969 *   phi_m.reinit(face);
970 *  
971 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
972 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
973 *   const auto& n_plus = phi_p.get_normal_vector(q); /*--- The normal vector is the same
974 *   for both phi_p and phi_m. If the face is interior,
975 *   it correspond to the outer normal ---*/
976 *  
977 *   const auto& avg_grad_u_old = 0.5*(phi_old_p.get_gradient(q) + phi_old_m.get_gradient(q));
978 *   const auto& avg_tensor_product_u_n = 0.5*(outer_product(phi_old_p.get_value(q), phi_old_extr_p.get_value(q)) +
979 *   outer_product(phi_old_m.get_value(q), phi_old_extr_m.get_value(q)));
980 *   const auto& avg_p_old = 0.5*(phi_old_press_p.get_value(q) + phi_old_press_m.get_value(q));
981 *  
982 *   phi_p.submit_value((a21/Re*avg_grad_u_old - a21*avg_tensor_product_u_n)*n_plus - avg_p_old*n_plus, q);
983 *   phi_m.submit_value(-(a21/Re*avg_grad_u_old - a21*avg_tensor_product_u_n)*n_plus + avg_p_old*n_plus, q);
984 *   }
985 *   phi_p.integrate_scatter(EvaluationFlags::values, dst);
986 *   phi_m.integrate_scatter(EvaluationFlags::values, dst);
987 *   }
988 *   }
989 *   else {
990 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
992 *   phi_m(data, false, 0),
993 *   phi_old_p(data, true, 0),
994 *   phi_old_m(data, false, 0),
995 *   phi_int_p(data, true, 0),
996 *   phi_int_m(data, false, 0);
997 *   FEFaceEvaluation<dim, fe_degree_p, n_q_points_1d_v, 1, Number> phi_old_press_p(data, true, 1),
998 *   phi_old_press_m(data, false, 1);
999 *  
1000 *   /*--- We loop over the faces in the range ---*/
1001 *   for(unsigned int face = face_range.first; face < face_range.second; ++ face) {
1002 *   phi_old_p.reinit(face);
1003 *   phi_old_p.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1004 *   phi_old_m.reinit(face);
1005 *   phi_old_m.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1006 *   phi_int_p.reinit(face);
1007 *   phi_int_p.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
1008 *   phi_int_m.reinit(face);
1009 *   phi_int_m.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
1010 *   phi_old_press_p.reinit(face);
1011 *   phi_old_press_p.gather_evaluate(src[2], EvaluationFlags::values);
1012 *   phi_old_press_m.reinit(face);
1013 *   phi_old_press_m.gather_evaluate(src[2], EvaluationFlags::values);
1014 *   phi_p.reinit(face);
1015 *   phi_m.reinit(face);
1016 *  
1017 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1018 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1019 *   const auto& n_plus = phi_p.get_normal_vector(q);
1020 *  
1021 *   const auto& avg_grad_u_old = 0.5*(phi_old_p.get_gradient(q) + phi_old_m.get_gradient(q));
1022 *   const auto& avg_grad_u_int = 0.5*(phi_int_p.get_gradient(q) + phi_int_m.get_gradient(q));
1023 *   const auto& avg_tensor_product_u_n = 0.5*(outer_product(phi_old_p.get_value(q), phi_old_p.get_value(q)) +
1024 *   outer_product(phi_old_m.get_value(q), phi_old_m.get_value(q)));
1025 *   const auto& avg_tensor_product_u_n_gamma = 0.5*(outer_product(phi_int_p.get_value(q), phi_int_p.get_value(q)) +
1026 *   outer_product(phi_int_m.get_value(q), phi_int_m.get_value(q)));
1027 *   const auto& avg_p_old = 0.5*(phi_old_press_p.get_value(q) + phi_old_press_m.get_value(q));
1028 *  
1029 *   phi_p.submit_value((a31/Re*avg_grad_u_old + a32/Re*avg_grad_u_int -
1030 *   a31*avg_tensor_product_u_n - a32*avg_tensor_product_u_n_gamma)*n_plus - avg_p_old*n_plus, q);
1031 *   phi_m.submit_value(-(a31/Re*avg_grad_u_old + a32/Re*avg_grad_u_int -
1032 *   a31*avg_tensor_product_u_n - a32*avg_tensor_product_u_n_gamma)*n_plus + avg_p_old*n_plus, q);
1033 *   }
1034 *   phi_p.integrate_scatter(EvaluationFlags::values, dst);
1035 *   phi_m.integrate_scatter(EvaluationFlags::values, dst);
1036 *   }
1037 *   }
1038 *   }
1039 *  
1040 *  
1041 * @endcode
1042 *
1043 * The followinf function assembles rhs boundary term for the velocity
1044 *
1045
1046 *
1047 *
1048 * @code
1049 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1050 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1051 *   assemble_rhs_boundary_term_velocity(const MatrixFree<dim, Number>& data,
1052 *   Vec& dst,
1053 *   const std::vector<Vec>& src,
1054 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1055 *   if(TR_BDF2_stage == 1) {
1056 *   /*--- We first start by declaring the suitable instances to read already available quantities. Clearly on the boundary
1057 *   the second argument has to be true. ---*/
1059 *   phi_old(data, true, 0),
1060 *   phi_old_extr(data, true, 0);
1062 *  
1063 *   /*--- We loop over the faces in the range ---*/
1064 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1065 *   phi_old.reinit(face);
1066 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1067 *   phi_old_extr.reinit(face);
1068 *   phi_old_extr.gather_evaluate(src[1], EvaluationFlags::values);
1069 *   phi_old_press.reinit(face);
1070 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
1071 *   phi.reinit(face);
1072 *  
1073 *   const auto boundary_id = data.get_boundary_id(face); /*--- Get the id in order to impose the proper boundary condition ---*/
1074 *   const auto coef_jump = (boundary_id == 1) ? 0.0 : C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1075 *   const double aux_coeff = (boundary_id == 1) ? 0.0 : 1.0;
1076 *  
1077 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1078 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1079 *   const auto& n_plus = phi.get_normal_vector(q);
1080 *  
1081 *   const auto& grad_u_old = phi_old.get_gradient(q);
1082 *   const auto& tensor_product_u_n = outer_product(phi_old.get_value(q), phi_old_extr.get_value(q));
1083 *   const auto& p_old = phi_old_press.get_value(q);
1084 *   const auto& point_vectorized = phi.quadrature_point(q);
1085 *   auto u_int_m = Tensor<1, dim, VectorizedArray<Number>>();
1086 *   if(boundary_id == 0) {
1087 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1088 *   Point<dim> point; /*--- The point returned by the 'quadrature_point' function is not an instance of Point
1089 *   and so it is not ready to be directly used. We need to pay attention to the
1090 *   vectorization ---*/
1091 *   for(unsigned int d = 0; d < dim; ++d)
1092 *   point[d] = point_vectorized[d][v];
1093 *   for(unsigned int d = 0; d < dim; ++d)
1094 *   u_int_m[d][v] = vel_boundary_inflow.value(point, d);
1095 *   }
1096 *   }
1097 *   const auto tensor_product_u_int_m = outer_product(u_int_m, phi_old_extr.get_value(q));
1098 *   const auto lambda = (boundary_id == 1) ? 0.0 : std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
1099 *  
1100 *   phi.submit_value((a21/Re*grad_u_old - a21*tensor_product_u_n)*n_plus - p_old*n_plus +
1101 *   a22/Re*2.0*coef_jump*u_int_m -
1102 *   aux_coeff*a22*tensor_product_u_int_m*n_plus + a22*lambda*u_int_m, q);
1103 *   phi.submit_normal_derivative(-aux_coeff*theta_v*a22/Re*u_int_m, q); /*--- This is equivalent to multiply to the gradient
1104 *   with outer product and use 'submit_gradient' ---*/
1105 *   }
1106 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1107 *   }
1108 *   }
1109 *   else {
1110 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1112 *   phi_old(data, true, 0),
1113 *   phi_int(data, true, 0),
1114 *   phi_int_extr(data, true, 0);
1116 *  
1117 *   /*--- We loop over the faces in the range ---*/
1118 *   for(unsigned int face = face_range.first; face < face_range.second; ++ face) {
1119 *   phi_old.reinit(face);
1120 *   phi_old.gather_evaluate(src[0], EvaluationFlags::values | EvaluationFlags::gradients);
1121 *   phi_int.reinit(face);
1122 *   phi_int.gather_evaluate(src[1], EvaluationFlags::values | EvaluationFlags::gradients);
1123 *   phi_old_press.reinit(face);
1124 *   phi_old_press.gather_evaluate(src[2], EvaluationFlags::values);
1125 *   phi_int_extr.reinit(face);
1126 *   phi_int_extr.gather_evaluate(src[3], EvaluationFlags::values);
1127 *   phi.reinit(face);
1128 *  
1129 *   const auto boundary_id = data.get_boundary_id(face);
1130 *   const auto coef_jump = (boundary_id == 1) ? 0.0 : C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1131 *   const double aux_coeff = (boundary_id == 1) ? 0.0 : 1.0;
1132 *  
1133 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1134 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1135 *   const auto& n_plus = phi.get_normal_vector(q);
1136 *  
1137 *   const auto& grad_u_old = phi_old.get_gradient(q);
1138 *   const auto& grad_u_int = phi_int.get_gradient(q);
1139 *   const auto& tensor_product_u_n = outer_product(phi_old.get_value(q), phi_old.get_value(q));
1140 *   const auto& tensor_product_u_n_gamma = outer_product(phi_int.get_value(q), phi_int.get_value(q));
1141 *   const auto& p_old = phi_old_press.get_value(q);
1142 *   const auto& point_vectorized = phi.quadrature_point(q);
1143 *   auto u_m = Tensor<1, dim, VectorizedArray<Number>>();
1144 *   if(boundary_id == 0) {
1145 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1146 *   Point<dim> point;
1147 *   for(unsigned int d = 0; d < dim; ++d)
1148 *   point[d] = point_vectorized[d][v];
1149 *   for(unsigned int d = 0; d < dim; ++d)
1150 *   u_m[d][v] = vel_boundary_inflow.value(point, d);
1151 *   }
1152 *   }
1153 *   const auto tensor_product_u_m = outer_product(u_m, phi_int_extr.get_value(q));
1154 *   const auto lambda = (boundary_id == 1) ? 0.0 : std::abs(scalar_product(phi_int_extr.get_value(q), n_plus));
1155 *  
1156 *   phi.submit_value((a31/Re*grad_u_old + a32/Re*grad_u_int -
1157 *   a31*tensor_product_u_n - a32*tensor_product_u_n_gamma)*n_plus - p_old*n_plus +
1158 *   a33/Re*2.0*coef_jump*u_m -
1159 *   aux_coeff*a33*tensor_product_u_m*n_plus + a33*lambda*u_m, q);
1160 *   phi.submit_normal_derivative(-aux_coeff*theta_v*a33/Re*u_m, q);
1161 *   }
1162 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1163 *   }
1164 *   }
1165 *   }
1166 *  
1167 *  
1168 * @endcode
1169 *
1170 * Put together all the previous steps for velocity. This is done automatically by the loop function of 'MatrixFree' class
1171 *
1172
1173 *
1174 *
1175 * @code
1176 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1177 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1178 *   vmult_rhs_velocity(Vec& dst, const std::vector<Vec>& src) const {
1179 *   for(auto& vec : src)
1180 *   vec.update_ghost_values();
1181 *  
1182 *   this->data->loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_velocity,
1183 *   &NavierStokesProjectionOperator::assemble_rhs_face_term_velocity,
1184 *   &NavierStokesProjectionOperator::assemble_rhs_boundary_term_velocity,
1185 *   this, dst, src, true,
1188 *   }
1189 *  
1190 *  
1191 * @endcode
1192 *
1193 * Now we focus on computing the rhs for the projection step for the pressure with the same ratio.
1194 * The following function assembles rhs cell term for the pressure
1195 *
1196
1197 *
1198 *
1199 * @code
1200 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1201 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1202 *   assemble_rhs_cell_term_pressure(const MatrixFree<dim, Number>& data,
1203 *   Vec& dst,
1204 *   const std::vector<Vec>& src,
1205 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1206 *   /*--- We first start by declaring the suitable instances to read already available quantities.
1207 *   The third parameter specifies that we want to use the second quadrature formula stored. ---*/
1209 *   phi_old(data, 1, 1);
1211 *  
1212 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0e6*gamma*dt*gamma*dt : 1.0e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
1213 *  
1214 *   const double coeff_2 = (TR_BDF2_stage == 1) ? gamma*dt : (1.0 - gamma)*dt;
1215 *  
1216 *   /*--- We loop over cells in the range ---*/
1217 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1218 *   phi_proj.reinit(cell);
1219 *   phi_proj.gather_evaluate(src[0], EvaluationFlags::values);
1220 *   phi_old.reinit(cell);
1221 *   phi_old.gather_evaluate(src[1], EvaluationFlags::values);
1222 *   phi.reinit(cell);
1223 *  
1224 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1225 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1226 *   const auto& u_star_star = phi_proj.get_value(q);
1227 *   const auto& p_old = phi_old.get_value(q);
1228 *  
1229 *   phi.submit_value(1.0/coeff*p_old, q);
1230 *   phi.submit_gradient(1.0/coeff_2*u_star_star, q);
1231 *   }
1232 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1233 *   }
1234 *   }
1235 *  
1236 *  
1237 * @endcode
1238 *
1239 * The following function assembles rhs face term for the pressure
1240 *
1241
1242 *
1243 *
1244 * @code
1245 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1246 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1247 *   assemble_rhs_face_term_pressure(const MatrixFree<dim, Number>& data,
1248 *   Vec& dst,
1249 *   const std::vector<Vec>& src,
1250 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1251 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1253 *   phi_m(data, false, 1, 1);
1255 *   phi_proj_m(data, false, 0, 1);
1256 *  
1257 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0/(gamma*dt) : 1.0/((1.0 - gamma)*dt);
1258 *  
1259 *   /*--- We loop over faces in the range ---*/
1260 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1261 *   phi_proj_p.reinit(face);
1262 *   phi_proj_p.gather_evaluate(src[0], EvaluationFlags::values);
1263 *   phi_proj_m.reinit(face);
1264 *   phi_proj_m.gather_evaluate(src[0], EvaluationFlags::values);
1265 *   phi_p.reinit(face);
1266 *   phi_m.reinit(face);
1267 *  
1268 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1269 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1270 *   const auto& n_plus = phi_p.get_normal_vector(q);
1271 *   const auto& avg_u_star_star = 0.5*(phi_proj_p.get_value(q) + phi_proj_m.get_value(q));
1272 *  
1273 *   phi_p.submit_value(-coeff*scalar_product(avg_u_star_star, n_plus), q);
1274 *   phi_m.submit_value(coeff*scalar_product(avg_u_star_star, n_plus), q);
1275 *   }
1276 *   phi_p.integrate_scatter(EvaluationFlags::values, dst);
1277 *   phi_m.integrate_scatter(EvaluationFlags::values, dst);
1278 *   }
1279 *   }
1280 *  
1281 *  
1282 * @endcode
1283 *
1284 * The following function assembles rhs boundary term for the pressure
1285 *
1286
1287 *
1288 *
1289 * @code
1290 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1291 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1292 *   assemble_rhs_boundary_term_pressure(const MatrixFree<dim, Number>& data,
1293 *   Vec& dst,
1294 *   const std::vector<Vec>& src,
1295 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1296 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1299 *  
1300 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0/(gamma*dt) : 1.0/((1.0 - gamma)*dt);
1301 *  
1302 *   /*--- We loop over faces in the range ---*/
1303 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1304 *   phi_proj.reinit(face);
1305 *   phi_proj.gather_evaluate(src[0], EvaluationFlags::values);
1306 *   phi.reinit(face);
1307 *  
1308 *   /*--- Now we loop over all the quadrature points to compute the integrals ---*/
1309 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1310 *   const auto& n_plus = phi.get_normal_vector(q);
1311 *  
1312 *   phi.submit_value(-coeff*scalar_product(phi_proj.get_value(q), n_plus), q);
1313 *   }
1314 *   phi.integrate_scatter(EvaluationFlags::values, dst);
1315 *   }
1316 *   }
1317 *  
1318 *  
1319 * @endcode
1320 *
1321 * Put together all the previous steps for pressure
1322 *
1323
1324 *
1325 *
1326 * @code
1327 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1328 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1329 *   vmult_rhs_pressure(Vec& dst, const std::vector<Vec>& src) const {
1330 *   for(auto& vec : src)
1331 *   vec.update_ghost_values();
1332 *  
1333 *   this->data->loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_pressure,
1334 *   &NavierStokesProjectionOperator::assemble_rhs_face_term_pressure,
1335 *   &NavierStokesProjectionOperator::assemble_rhs_boundary_term_pressure,
1336 *   this, dst, src, true,
1339 *   }
1340 *  
1341 *  
1342 * @endcode
1343 *
1344 * Now we need to build the 'matrices', i.e. the bilinear forms. We start by
1345 * assembling the cell term for the velocity
1346 *
1347
1348 *
1349 *
1350 * @code
1351 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1352 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1353 *   assemble_cell_term_velocity(const MatrixFree<dim, Number>& data,
1354 *   Vec& dst,
1355 *   const Vec& src,
1356 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1357 *   if(TR_BDF2_stage == 1) {
1358 *   /*--- We first start by declaring the suitable instances to read already available quantities. Moreover 'phi' in
1359 *   this case serves for a bilinear form and so it will not used only to submit but also to read the src ---*/
1361 *   phi_old_extr(data, 0);
1362 *  
1363 *   /*--- We loop over all cells in the range ---*/
1364 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1365 *   phi.reinit(cell);
1366 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1367 *   phi_old_extr.reinit(cell);
1368 *   phi_old_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1369 *  
1370 *   /*--- Now we loop over all quadrature points ---*/
1371 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1372 *   const auto& u_int = phi.get_value(q);
1373 *   const auto& grad_u_int = phi.get_gradient(q);
1374 *   const auto& u_n_gamma_ov_2 = phi_old_extr.get_value(q);
1375 *   const auto& tensor_product_u_int = outer_product(u_int, u_n_gamma_ov_2);
1376 *  
1377 *   phi.submit_value(1.0/(gamma*dt)*u_int, q);
1378 *   phi.submit_gradient(-a22*tensor_product_u_int + a22/Re*grad_u_int, q);
1379 *   }
1380 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1381 *   }
1382 *   }
1383 *   else {
1384 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1386 *   phi_int_extr(data, 0);
1387 *  
1388 *   /*--- We loop over all cells in the range ---*/
1389 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1390 *   phi.reinit(cell);
1391 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1392 *   phi_int_extr.reinit(cell);
1393 *   phi_int_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1394 *  
1395 *   /*--- Now we loop over all quadrature points ---*/
1396 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1397 *   const auto& u_curr = phi.get_value(q);
1398 *   const auto& grad_u_curr = phi.get_gradient(q);
1399 *   const auto& u_n1_int = phi_int_extr.get_value(q);
1400 *   const auto& tensor_product_u_curr = outer_product(u_curr, u_n1_int);
1401 *  
1402 *   phi.submit_value(1.0/((1.0 - gamma)*dt)*u_curr, q);
1403 *   phi.submit_gradient(-a33*tensor_product_u_curr + a33/Re*grad_u_curr, q);
1404 *   }
1405 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1406 *   }
1407 *   }
1408 *   }
1409 *  
1410 *  
1411 * @endcode
1412 *
1413 * The following function assembles face term for the velocity
1414 *
1415
1416 *
1417 *
1418 * @code
1419 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1420 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1421 *   assemble_face_term_velocity(const MatrixFree<dim, Number>& data,
1422 *   Vec& dst,
1423 *   const Vec& src,
1424 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1425 *   if(TR_BDF2_stage == 1) {
1426 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1428 *   phi_m(data, false, 0),
1429 *   phi_old_extr_p(data, true, 0),
1430 *   phi_old_extr_m(data, false, 0);
1431 *  
1432 *   /*--- We loop over all faces in the range ---*/
1433 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1434 *   phi_p.reinit(face);
1435 *   phi_p.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1436 *   phi_m.reinit(face);
1437 *   phi_m.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1438 *   phi_old_extr_p.reinit(face);
1439 *   phi_old_extr_p.gather_evaluate(u_extr, EvaluationFlags::values);
1440 *   phi_old_extr_m.reinit(face);
1441 *   phi_old_extr_m.gather_evaluate(u_extr, EvaluationFlags::values);
1442 *  
1443 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
1444 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
1445 *  
1446 *   /*--- Now we loop over all quadrature points ---*/
1447 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1448 *   const auto& n_plus = phi_p.get_normal_vector(q);
1449 *  
1450 *   const auto& avg_grad_u_int = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
1451 *   const auto& jump_u_int = phi_p.get_value(q) - phi_m.get_value(q);
1452 *   const auto& avg_tensor_product_u_int = 0.5*(outer_product(phi_p.get_value(q), phi_old_extr_p.get_value(q)) +
1453 *   outer_product(phi_m.get_value(q), phi_old_extr_m.get_value(q)));
1454 *   const auto lambda = std::max(std::abs(scalar_product(phi_old_extr_p.get_value(q), n_plus)),
1455 *   std::abs(scalar_product(phi_old_extr_m.get_value(q), n_plus)));
1456 *  
1457 *   phi_p.submit_value(a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) +
1458 *   a22*avg_tensor_product_u_int*n_plus + 0.5*a22*lambda*jump_u_int, q);
1459 *   phi_m.submit_value(-a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) -
1460 *   a22*avg_tensor_product_u_int*n_plus - 0.5*a22*lambda*jump_u_int, q);
1461 *   phi_p.submit_normal_derivative(-theta_v*a22/Re*0.5*jump_u_int, q);
1462 *   phi_m.submit_normal_derivative(-theta_v*a22/Re*0.5*jump_u_int, q);
1463 *   }
1464 *   phi_p.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1465 *   phi_m.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1466 *   }
1467 *   }
1468 *   else {
1469 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1471 *   phi_m(data, false, 0),
1472 *   phi_extr_p(data, true, 0),
1473 *   phi_extr_m(data, false, 0);
1474 *  
1475 *   /*--- We loop over all faces in the range ---*/
1476 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1477 *   phi_p.reinit(face);
1478 *   phi_p.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1479 *   phi_m.reinit(face);
1480 *   phi_m.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1481 *   phi_extr_p.reinit(face);
1482 *   phi_extr_p.gather_evaluate(u_extr, EvaluationFlags::values);
1483 *   phi_extr_m.reinit(face);
1484 *   phi_extr_m.gather_evaluate(u_extr, EvaluationFlags::values);
1485 *  
1486 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
1487 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
1488 *  
1489 *   /*--- Now we loop over all quadrature points ---*/
1490 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1491 *   const auto& n_plus = phi_p.get_normal_vector(q);
1492 *  
1493 *   const auto& avg_grad_u = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
1494 *   const auto& jump_u = phi_p.get_value(q) - phi_m.get_value(q);
1495 *   const auto& avg_tensor_product_u = 0.5*(outer_product(phi_p.get_value(q), phi_extr_p.get_value(q)) +
1496 *   outer_product(phi_m.get_value(q), phi_extr_m.get_value(q)));
1497 *   const auto lambda = std::max(std::abs(scalar_product(phi_extr_p.get_value(q), n_plus)),
1498 *   std::abs(scalar_product(phi_extr_m.get_value(q), n_plus)));
1499 *  
1500 *   phi_p.submit_value(a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) +
1501 *   a33*avg_tensor_product_u*n_plus + 0.5*a33*lambda*jump_u, q);
1502 *   phi_m.submit_value(-a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) -
1503 *   a33*avg_tensor_product_u*n_plus - 0.5*a33*lambda*jump_u, q);
1504 *   phi_p.submit_normal_derivative(-theta_v*a33/Re*0.5*jump_u, q);
1505 *   phi_m.submit_normal_derivative(-theta_v*a33/Re*0.5*jump_u, q);
1506 *   }
1507 *   phi_p.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1508 *   phi_m.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1509 *   }
1510 *   }
1511 *   }
1512 *  
1513 *  
1514 * @endcode
1515 *
1516 * The following function assembles boundary term for the velocity
1517 *
1518
1519 *
1520 *
1521 * @code
1522 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1523 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1524 *   assemble_boundary_term_velocity(const MatrixFree<dim, Number>& data,
1525 *   Vec& dst,
1526 *   const Vec& src,
1527 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1528 *   if(TR_BDF2_stage == 1) {
1529 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1531 *   phi_old_extr(data, true, 0);
1532 *  
1533 *   /*--- We loop over all faces in the range ---*/
1534 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1535 *   phi.reinit(face);
1536 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1537 *   phi_old_extr.reinit(face);
1538 *   phi_old_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1539 *  
1540 *   const auto boundary_id = data.get_boundary_id(face);
1541 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1542 *  
1543 *   /*--- The application of the mirror principle is not so trivial because we have a Dirichlet condition
1544 *   on a single component for the outflow; so we distinguish the two cases ---*/
1545 *   if(boundary_id != 1) {
1546 *   const double coef_trasp = 0.0;
1547 *  
1548 *   /*--- Now we loop over all quadrature points ---*/
1549 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1550 *   const auto& n_plus = phi.get_normal_vector(q);
1551 *   const auto& grad_u_int = phi.get_gradient(q);
1552 *   const auto& u_int = phi.get_value(q);
1553 *   const auto& tensor_product_u_int = outer_product(phi.get_value(q), phi_old_extr.get_value(q));
1554 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
1555 *  
1556 *   phi.submit_value(a22/Re*(-grad_u_int*n_plus + 2.0*coef_jump*u_int) +
1557 *   a22*coef_trasp*tensor_product_u_int*n_plus + a22*lambda*u_int, q);
1558 *   phi.submit_normal_derivative(-theta_v*a22/Re*u_int, q);
1559 *   }
1560 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1561 *   }
1562 *   else {
1563 *   /*--- Now we loop over all quadrature points ---*/
1564 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1565 *   const auto& n_plus = phi.get_normal_vector(q);
1566 *   const auto& grad_u_int = phi.get_gradient(q);
1567 *   const auto& u_int = phi.get_value(q);
1568 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
1569 *  
1570 *   const auto& point_vectorized = phi.quadrature_point(q);
1571 *   auto u_int_m = u_int;
1572 *   auto grad_u_int_m = grad_u_int;
1573 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1574 *   Point<dim> point;
1575 *   for(unsigned int d = 0; d < dim; ++d)
1576 *   point[d] = point_vectorized[d][v];
1577 *  
1578 *   u_int_m[1][v] = -u_int_m[1][v];
1579 *  
1580 *   grad_u_int_m[0][0][v] = -grad_u_int_m[0][0][v];
1581 *   grad_u_int_m[0][1][v] = -grad_u_int_m[0][1][v];
1582 *   }
1583 *  
1584 *   phi.submit_value(a22/Re*(-(0.5*(grad_u_int + grad_u_int_m))*n_plus + coef_jump*(u_int - u_int_m)) +
1585 *   a22*outer_product(0.5*(u_int + u_int_m), phi_old_extr.get_value(q))*n_plus +
1586 *   a22*0.5*lambda*(u_int - u_int_m), q);
1587 *   phi.submit_normal_derivative(-theta_v*a22/Re*(u_int - u_int_m), q);
1588 *   }
1589 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1590 *   }
1591 *   }
1592 *   }
1593 *   else {
1594 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1596 *   phi_extr(data, true, 0);
1597 *  
1598 *   /*--- We loop over all faces in the range ---*/
1599 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1600 *   phi.reinit(face);
1601 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1602 *   phi_extr.reinit(face);
1603 *   phi_extr.gather_evaluate(u_extr, EvaluationFlags::values);
1604 *  
1605 *   const auto boundary_id = data.get_boundary_id(face);
1606 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
1607 *  
1608 *   if(boundary_id != 1) {
1609 *   const double coef_trasp = 0.0;
1610 *  
1611 *   /*--- Now we loop over all quadrature points ---*/
1612 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1613 *   const auto& n_plus = phi.get_normal_vector(q);
1614 *   const auto& grad_u = phi.get_gradient(q);
1615 *   const auto& u = phi.get_value(q);
1616 *   const auto& tensor_product_u = outer_product(phi.get_value(q), phi_extr.get_value(q));
1617 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
1618 *  
1619 *   phi.submit_value(a33/Re*(-grad_u*n_plus + 2.0*coef_jump*u) +
1620 *   a33*coef_trasp*tensor_product_u*n_plus + a33*lambda*u, q);
1621 *   phi.submit_normal_derivative(-theta_v*a33/Re*u, q);
1622 *   }
1623 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1624 *   }
1625 *   else {
1626 *   /*--- Now we loop over all quadrature points ---*/
1627 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1628 *   const auto& n_plus = phi.get_normal_vector(q);
1629 *   const auto& grad_u = phi.get_gradient(q);
1630 *   const auto& u = phi.get_value(q);
1631 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
1632 *  
1633 *   const auto& point_vectorized = phi.quadrature_point(q);
1634 *   auto u_m = u;
1635 *   auto grad_u_m = grad_u;
1636 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
1637 *   Point<dim> point;
1638 *   for(unsigned int d = 0; d < dim; ++d)
1639 *   point[d] = point_vectorized[d][v];
1640 *  
1641 *   u_m[1][v] = -u_m[1][v];
1642 *  
1643 *   grad_u_m[0][0][v] = -grad_u_m[0][0][v];
1644 *   grad_u_m[0][1][v] = -grad_u_m[0][1][v];
1645 *   }
1646 *  
1647 *   phi.submit_value(a33/Re*(-(0.5*(grad_u + grad_u_m))*n_plus + coef_jump*(u - u_m)) +
1648 *   a33*outer_product(0.5*(u + u_m), phi_extr.get_value(q))*n_plus + a33*0.5*lambda*(u - u_m), q);
1649 *   phi.submit_normal_derivative(-theta_v*a33/Re*(u - u_m), q);
1650 *   }
1651 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1652 *   }
1653 *   }
1654 *   }
1655 *   }
1656 *  
1657 *  
1658 * @endcode
1659 *
1660 * Next, we focus on 'matrices' to compute the pressure. We first assemble cell term for the pressure
1661 *
1662
1663 *
1664 *
1665 * @code
1666 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1667 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1668 *   assemble_cell_term_pressure(const MatrixFree<dim, Number>& data,
1669 *   Vec& dst,
1670 *   const Vec& src,
1671 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1672 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1674 *  
1675 *   const double coeff = (TR_BDF2_stage == 1) ? 1.0e6*gamma*dt*gamma*dt : 1.0e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
1676 *  
1677 *   /*--- Loop over all cells in the range ---*/
1678 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1679 *   phi.reinit(cell);
1680 *   phi.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1681 *  
1682 *   /*--- Now we loop over all quadrature points ---*/
1683 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1684 *   phi.submit_gradient(phi.get_gradient(q), q);
1685 *   phi.submit_value(1.0/coeff*phi.get_value(q), q);
1686 *   }
1687 *  
1688 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1689 *   }
1690 *   }
1691 *  
1692 *  
1693 * @endcode
1694 *
1695 * The following function assembles face term for the pressure
1696 *
1697
1698 *
1699 *
1700 * @code
1701 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1702 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1703 *   assemble_face_term_pressure(const MatrixFree<dim, Number>& data,
1704 *   Vec& dst,
1705 *   const Vec& src,
1706 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1707 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1709 *   phi_m(data, false, 1, 1);
1710 *  
1711 *   /*--- Loop over all faces in the range ---*/
1712 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1713 *   phi_p.reinit(face);
1714 *   phi_p.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1715 *   phi_m.reinit(face);
1716 *   phi_m.gather_evaluate(src, EvaluationFlags::values | EvaluationFlags::gradients);
1717 *  
1718 *   const auto coef_jump = C_p*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
1719 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
1720 *  
1721 *   /*--- Loop over quadrature points ---*/
1722 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
1723 *   const auto& n_plus = phi_p.get_normal_vector(q);
1724 *  
1725 *   const auto& avg_grad_pres = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
1726 *   const auto& jump_pres = phi_p.get_value(q) - phi_m.get_value(q);
1727 *  
1728 *   phi_p.submit_value(-scalar_product(avg_grad_pres, n_plus) + coef_jump*jump_pres, q);
1729 *   phi_m.submit_value(scalar_product(avg_grad_pres, n_plus) - coef_jump*jump_pres, q);
1730 *   phi_p.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
1731 *   phi_m.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
1732 *   }
1733 *   phi_p.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1734 *   phi_m.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1735 *   }
1736 *   }
1737 *  
1738 *  
1739 * @endcode
1740 *
1741 * The following function assembles boundary term for the pressure
1742 *
1743
1744 *
1745 *
1746 * @code
1747 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1748 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1749 *   assemble_boundary_term_pressure(const MatrixFree<dim, Number>& data,
1750 *   Vec& dst,
1751 *   const Vec& src,
1752 *   const std::pair<unsigned int, unsigned int>& face_range) const {
1754 *  
1755 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
1756 *   const auto boundary_id = data.get_boundary_id(face);
1757 *  
1758 *   if(boundary_id == 1) {
1759 *   phi.reinit(face);
1760 *   phi.gather_evaluate(src, true, true);
1761 *  
1762 *   const auto coef_jump = C_p*std::abs((phi.get_normal_vector(0)*phi.inverse_jacobian(0))[dim - 1]);
1763 *  
1764 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1765 *   const auto& n_plus = phi.get_normal_vector(q);
1766 *  
1767 *   const auto& grad_pres = phi.get_gradient(q);
1768 *   const auto& pres = phi.get_value(q);
1769 *  
1770 *   phi.submit_value(-scalar_product(grad_pres, n_plus) + coef_jump*pres , q);
1771 *   phi.submit_normal_derivative(-theta_p*pres, q);
1772 *   }
1773 *   phi.integrate_scatter(EvaluationFlags::values | EvaluationFlags::gradients, dst);
1774 *   }
1775 *   }
1776 *   }
1777 *  
1778 *  
1779 * @endcode
1780 *
1781 * Before coding the 'apply_add' function, which is the one that will perform the loop, we focus on
1782 * the linear system that arises to project the gradient of the pressure into the velocity space.
1783 * The following function assembles rhs cell term for the projection of gradient of pressure. Since no
1784 * integration by parts is performed, only a cell term contribution is present.
1785 *
1786
1787 *
1788 *
1789 * @code
1790 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1791 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1792 *   assemble_rhs_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
1793 *   Vec& dst,
1794 *   const Vec& src,
1795 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1796 *   /*--- We first start by declaring the suitable instances to read already available quantities. ---*/
1799 *  
1800 *   /*--- Loop over all cells in the range ---*/
1801 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1802 *   phi_pres.reinit(cell);
1803 *   phi_pres.gather_evaluate(src, EvaluationFlags::gradients);
1804 *   phi.reinit(cell);
1805 *  
1806 *   /*--- Loop over quadrature points ---*/
1807 *   for(unsigned int q = 0; q < phi.n_q_points; ++q)
1808 *   phi.submit_value(phi_pres.get_gradient(q), q);
1809 *  
1810 *   phi.integrate_scatter(EvaluationFlags::values, dst);
1811 *   }
1812 *   }
1813 *  
1814 *  
1815 * @endcode
1816 *
1817 * Put together all the previous steps for porjection of pressure gradient. Here we loop only over cells
1818 *
1819
1820 *
1821 *
1822 * @code
1823 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1824 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1825 *   vmult_grad_p_projection(Vec& dst, const Vec& src) const {
1826 *   this->data->cell_loop(&NavierStokesProjectionOperator::assemble_rhs_cell_term_projection_grad_p,
1827 *   this, dst, src, true);
1828 *   }
1829 *  
1830 *  
1831 * @endcode
1832 *
1833 * Assemble now cell term for the projection of gradient of pressure. This is nothing but a mass matrix
1834 *
1835
1836 *
1837 *
1838 * @code
1839 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1840 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1841 *   assemble_cell_term_projection_grad_p(const MatrixFree<dim, Number>& data,
1842 *   Vec& dst,
1843 *   const Vec& src,
1844 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1846 *  
1847 *   /*--- Loop over all cells in the range ---*/
1848 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1849 *   phi.reinit(cell);
1850 *   phi.gather_evaluate(src, EvaluationFlags::values);
1851 *  
1852 *   /*--- Loop over quadrature points ---*/
1853 *   for(unsigned int q = 0; q < phi.n_q_points; ++q)
1854 *   phi.submit_value(phi.get_value(q), q);
1855 *  
1856 *   phi.integrate_scatter(EvaluationFlags::values, dst);
1857 *   }
1858 *   }
1859 *  
1860 *  
1861 * @endcode
1862 *
1863 * Put together all previous steps. This is the overriden function that effectively performs the
1864 * matrix-vector multiplication.
1865 *
1866
1867 *
1868 *
1869 * @code
1870 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1871 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1872 *   apply_add(Vec& dst, const Vec& src) const {
1873 *   if(NS_stage == 1) {
1874 *   this->data->loop(&NavierStokesProjectionOperator::assemble_cell_term_velocity,
1875 *   &NavierStokesProjectionOperator::assemble_face_term_velocity,
1876 *   &NavierStokesProjectionOperator::assemble_boundary_term_velocity,
1877 *   this, dst, src, false,
1880 *   }
1881 *   else if(NS_stage == 2) {
1882 *   this->data->loop(&NavierStokesProjectionOperator::assemble_cell_term_pressure,
1883 *   &NavierStokesProjectionOperator::assemble_face_term_pressure,
1884 *   &NavierStokesProjectionOperator::assemble_boundary_term_pressure,
1885 *   this, dst, src, false,
1888 *   }
1889 *   else if(NS_stage == 3) {
1890 *   this->data->cell_loop(&NavierStokesProjectionOperator::assemble_cell_term_projection_grad_p,
1891 *   this, dst, src, false); /*--- Since we have only a cell term contribution, we use cell_loop ---*/
1892 *   }
1893 *   else
1894 *   Assert(false, ExcNotImplemented());
1895 *   }
1896 *  
1897 *  
1898 * @endcode
1899 *
1900 * Finally, we focus on computing the diagonal for preconditioners and we start by assembling
1901 * the diagonal cell term for the velocity. Since we do not have access to the entries of the matrix,
1902 * in order to compute the element i, we test the matrix against a vector which is equal to 1 in position i and 0 elsewhere.
1903 * This is why 'src' will result as unused.
1904 *
1905
1906 *
1907 *
1908 * @code
1909 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
1910 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
1911 *   assemble_diagonal_cell_term_velocity(const MatrixFree<dim, Number>& data,
1912 *   Vec& dst,
1913 *   const unsigned int& ,
1914 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
1915 *   if(TR_BDF2_stage == 1) {
1917 *   phi_old_extr(data, 0);
1918 *  
1920 *   /*--- Build a vector of ones to be tested (here we will see the velocity as a whole vector, since
1921 *   dof_handler_velocity is vectorial and so the dof values are vectors). ---*/
1923 *   for(unsigned int d = 0; d < dim; ++d)
1924 *   tmp[d] = make_vectorized_array<Number>(1.0);
1925 *  
1926 *   /*--- Loop over cells in the range ---*/
1927 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1928 *   phi_old_extr.reinit(cell);
1929 *   phi_old_extr.gather_evaluate(u_extr, true, false);
1930 *   phi.reinit(cell);
1931 *  
1932 *   /*--- Loop over dofs ---*/
1933 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
1934 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
1935 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j); /*--- Set all dofs to zero ---*/
1936 *   phi.submit_dof_value(tmp, i); /*--- Set dof i equal to one ---*/
1937 *   phi.evaluate(true, true);
1938 *  
1939 *   /*--- Loop over quadrature points ---*/
1940 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1941 *   const auto& u_int = phi.get_value(q);
1942 *   const auto& grad_u_int = phi.get_gradient(q);
1943 *   const auto& u_n_gamma_ov_2 = phi_old_extr.get_value(q);
1944 *   const auto& tensor_product_u_int = outer_product(u_int, u_n_gamma_ov_2);
1945 *  
1946 *   phi.submit_value(1.0/(gamma*dt)*u_int, q);
1947 *   phi.submit_gradient(-a22*tensor_product_u_int + a22/Re*grad_u_int, q);
1948 *   }
1949 *   phi.integrate(true, true);
1950 *   diagonal[i] = phi.get_dof_value(i);
1951 *   }
1952 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
1953 *   phi.submit_dof_value(diagonal[i], i);
1954 *   phi.distribute_local_to_global(dst);
1955 *   }
1956 *   }
1957 *   else {
1959 *   phi_int_extr(data, 0);
1960 *  
1963 *   for(unsigned int d = 0; d < dim; ++d)
1964 *   tmp[d] = make_vectorized_array<Number>(1.0);
1965 *  
1966 *   /*--- Loop over cells in the range ---*/
1967 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
1968 *   phi_int_extr.reinit(cell);
1969 *   phi_int_extr.gather_evaluate(u_extr, true, false);
1970 *   phi.reinit(cell);
1971 *  
1972 *   /*--- Loop over dofs ---*/
1973 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
1974 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
1975 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
1976 *   phi.submit_dof_value(tmp, i);
1977 *   phi.evaluate(true, true);
1978 *  
1979 *   /*--- Loop over quadrature points ---*/
1980 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
1981 *   const auto& u_curr = phi.get_value(q);
1982 *   const auto& grad_u_curr = phi.get_gradient(q);
1983 *   const auto& u_n1_int = phi_int_extr.get_value(q);
1984 *   const auto& tensor_product_u_curr = outer_product(u_curr, u_n1_int);
1985 *  
1986 *   phi.submit_value(1.0/((1.0 - gamma)*dt)*u_curr, q);
1987 *   phi.submit_gradient(-a33*tensor_product_u_curr + a33/Re*grad_u_curr, q);
1988 *   }
1989 *   phi.integrate(true, true);
1990 *   diagonal[i] = phi.get_dof_value(i);
1991 *   }
1992 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
1993 *   phi.submit_dof_value(diagonal[i], i);
1994 *   phi.distribute_local_to_global(dst);
1995 *   }
1996 *   }
1997 *   }
1998 *  
1999 *  
2000 * @endcode
2001 *
2002 * The following function assembles diagonal face term for the velocity
2003 *
2004
2005 *
2006 *
2007 * @code
2008 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
2009 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2010 *   assemble_diagonal_face_term_velocity(const MatrixFree<dim, Number>& data,
2011 *   Vec& dst,
2012 *   const unsigned int& ,
2013 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2014 *   if(TR_BDF2_stage == 1) {
2016 *   phi_m(data, false, 0),
2017 *   phi_old_extr_p(data, true, 0),
2018 *   phi_old_extr_m(data, false, 0);
2019 *  
2020 *   AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component); /*--- We just assert for safety that dimension match,
2021 *   in the sense that we have selected the proper
2022 *   space ---*/
2023 *   AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal_p(phi_p.dofs_per_component),
2024 *   diagonal_m(phi_m.dofs_per_component);
2025 *  
2027 *   for(unsigned int d = 0; d < dim; ++d)
2028 *   tmp[d] = make_vectorized_array<Number>(1.0); /*--- We build the usal vector of ones that we will use as dof value ---*/
2029 *  
2030 *   /*--- Now we loop over faces ---*/
2031 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2032 *   phi_old_extr_p.reinit(face);
2033 *   phi_old_extr_p.gather_evaluate(u_extr, true, false);
2034 *   phi_old_extr_m.reinit(face);
2035 *   phi_old_extr_m.gather_evaluate(u_extr, true, false);
2036 *   phi_p.reinit(face);
2037 *   phi_m.reinit(face);
2038 *  
2039 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
2040 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
2041 *  
2042 *   /*--- Loop over dofs. We will set all equal to zero apart from the current one ---*/
2043 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2044 *   for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
2045 *   phi_p.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2046 *   phi_m.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2047 *   }
2048 *   phi_p.submit_dof_value(tmp, i);
2049 *   phi_p.evaluate(true, true);
2050 *   phi_m.submit_dof_value(tmp, i);
2051 *   phi_m.evaluate(true, true);
2052 *  
2053 *   /*--- Loop over quadrature points to compute the integral ---*/
2054 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
2055 *   const auto& n_plus = phi_p.get_normal_vector(q);
2056 *   const auto& avg_grad_u_int = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
2057 *   const auto& jump_u_int = phi_p.get_value(q) - phi_m.get_value(q);
2058 *   const auto& avg_tensor_product_u_int = 0.5*(outer_product(phi_p.get_value(q), phi_old_extr_p.get_value(q)) +
2059 *   outer_product(phi_m.get_value(q), phi_old_extr_m.get_value(q)));
2060 *   const auto lambda = std::max(std::abs(scalar_product(phi_old_extr_p.get_value(q), n_plus)),
2061 *   std::abs(scalar_product(phi_old_extr_m.get_value(q), n_plus)));
2062 *  
2063 *   phi_p.submit_value(a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) +
2064 *   a22*avg_tensor_product_u_int*n_plus + 0.5*a22*lambda*jump_u_int , q);
2065 *   phi_m.submit_value(-a22/Re*(-avg_grad_u_int*n_plus + coef_jump*jump_u_int) -
2066 *   a22*avg_tensor_product_u_int*n_plus - 0.5*a22*lambda*jump_u_int, q);
2067 *   phi_p.submit_normal_derivative(-theta_v*0.5*a22/Re*jump_u_int, q);
2068 *   phi_m.submit_normal_derivative(-theta_v*0.5*a22/Re*jump_u_int, q);
2069 *   }
2070 *   phi_p.integrate(true, true);
2071 *   diagonal_p[i] = phi_p.get_dof_value(i);
2072 *   phi_m.integrate(true, true);
2073 *   diagonal_m[i] = phi_m.get_dof_value(i);
2074 *   }
2075 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2076 *   phi_p.submit_dof_value(diagonal_p[i], i);
2077 *   phi_m.submit_dof_value(diagonal_m[i], i);
2078 *   }
2079 *   phi_p.distribute_local_to_global(dst);
2080 *   phi_m.distribute_local_to_global(dst);
2081 *   }
2082 *   }
2083 *   else {
2085 *   phi_m(data, false, 0),
2086 *   phi_extr_p(data, true, 0),
2087 *   phi_extr_m(data, false, 0);
2088 *  
2089 *   AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component);
2090 *   AlignedVector<Tensor<1, dim, VectorizedArray<Number>>> diagonal_p(phi_p.dofs_per_component),
2091 *   diagonal_m(phi_m.dofs_per_component);
2093 *   for(unsigned int d = 0; d < dim; ++d)
2094 *   tmp[d] = make_vectorized_array<Number>(1.0);
2095 *  
2096 *   /*--- Now we loop over faces ---*/
2097 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2098 *   phi_extr_p.reinit(face);
2099 *   phi_extr_p.gather_evaluate(u_extr, true, false);
2100 *   phi_extr_m.reinit(face);
2101 *   phi_extr_m.gather_evaluate(u_extr, true, false);
2102 *   phi_p.reinit(face);
2103 *   phi_m.reinit(face);
2104 *  
2105 *   const auto coef_jump = C_u*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
2106 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
2107 *  
2108 *   /*--- Loop over dofs. We will set all equal to zero apart from the current one ---*/
2109 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2110 *   for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
2111 *   phi_p.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2112 *   phi_m.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2113 *   }
2114 *   phi_p.submit_dof_value(tmp, i);
2115 *   phi_p.evaluate(true, true);
2116 *   phi_m.submit_dof_value(tmp, i);
2117 *   phi_m.evaluate(true, true);
2118 *  
2119 *   /*--- Loop over quadrature points to compute the integral ---*/
2120 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
2121 *   const auto& n_plus = phi_p.get_normal_vector(q);
2122 *   const auto& avg_grad_u = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
2123 *   const auto& jump_u = phi_p.get_value(q) - phi_m.get_value(q);
2124 *   const auto& avg_tensor_product_u = 0.5*(outer_product(phi_p.get_value(q), phi_extr_p.get_value(q)) +
2125 *   outer_product(phi_m.get_value(q), phi_extr_m.get_value(q)));
2126 *   const auto lambda = std::max(std::abs(scalar_product(phi_extr_p.get_value(q), n_plus)),
2127 *   std::abs(scalar_product(phi_extr_m.get_value(q), n_plus)));
2128 *  
2129 *   phi_p.submit_value(a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) +
2130 *   a33*avg_tensor_product_u*n_plus + 0.5*a33*lambda*jump_u, q);
2131 *   phi_m.submit_value(-a33/Re*(-avg_grad_u*n_plus + coef_jump*jump_u) -
2132 *   a33*avg_tensor_product_u*n_plus - 0.5*a33*lambda*jump_u, q);
2133 *   phi_p.submit_normal_derivative(-theta_v*0.5*a33/Re*jump_u, q);
2134 *   phi_m.submit_normal_derivative(-theta_v*0.5*a33/Re*jump_u, q);
2135 *   }
2136 *   phi_p.integrate(true, true);
2137 *   diagonal_p[i] = phi_p.get_dof_value(i);
2138 *   phi_m.integrate(true, true);
2139 *   diagonal_m[i] = phi_m.get_dof_value(i);
2140 *   }
2141 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2142 *   phi_p.submit_dof_value(diagonal_p[i], i);
2143 *   phi_m.submit_dof_value(diagonal_m[i], i);
2144 *   }
2145 *   phi_p.distribute_local_to_global(dst);
2146 *   phi_m.distribute_local_to_global(dst);
2147 *   }
2148 *   }
2149 *   }
2150 *  
2151 *  
2152 * @endcode
2153 *
2154 * The following function assembles boundary term for the velocity
2155 *
2156
2157 *
2158 *
2159 * @code
2160 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
2161 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2162 *   assemble_diagonal_boundary_term_velocity(const MatrixFree<dim, Number>& data,
2163 *   Vec& dst,
2164 *   const unsigned int& ,
2165 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2166 *   if(TR_BDF2_stage == 1) {
2168 *   phi_old_extr(data, true, 0);
2169 *  
2172 *   for(unsigned int d = 0; d < dim; ++d)
2173 *   tmp[d] = make_vectorized_array<Number>(1.0);
2174 *  
2175 *   /*--- Loop over all faces in the range ---*/
2176 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2177 *   phi_old_extr.reinit(face);
2178 *   phi_old_extr.gather_evaluate(u_extr, true, false);
2179 *   phi.reinit(face);
2180 *  
2181 *   const auto boundary_id = data.get_boundary_id(face);
2182 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
2183 *  
2184 *   if(boundary_id != 1) {
2185 *   const double coef_trasp = 0.0;
2186 *  
2187 *   /*--- Loop over all dofs ---*/
2188 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2189 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2190 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2191 *   phi.submit_dof_value(tmp, i);
2192 *   phi.evaluate(true, true);
2193 *  
2194 *   /*--- Loop over quadrature points to compute the integral ---*/
2195 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2196 *   const auto& n_plus = phi.get_normal_vector(q);
2197 *   const auto& grad_u_int = phi.get_gradient(q);
2198 *   const auto& u_int = phi.get_value(q);
2199 *   const auto& tensor_product_u_int = outer_product(phi.get_value(q), phi_old_extr.get_value(q));
2200 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
2201 *  
2202 *   phi.submit_value(a22/Re*(-grad_u_int*n_plus + 2.0*coef_jump*u_int) +
2203 *   a22*coef_trasp*tensor_product_u_int*n_plus + a22*lambda*u_int, q);
2204 *   phi.submit_normal_derivative(-theta_v*a22/Re*u_int, q);
2205 *   }
2206 *   phi.integrate(true, true);
2207 *   diagonal[i] = phi.get_dof_value(i);
2208 *   }
2209 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2210 *   phi.submit_dof_value(diagonal[i], i);
2211 *   phi.distribute_local_to_global(dst);
2212 *   }
2213 *   else {
2214 *   /*--- Loop over all dofs ---*/
2215 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2216 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2217 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2218 *   phi.submit_dof_value(tmp, i);
2219 *   phi.evaluate(true, true);
2220 *  
2221 *   /*--- Loop over quadrature points to compute the integral ---*/
2222 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2223 *   const auto& n_plus = phi.get_normal_vector(q);
2224 *   const auto& grad_u_int = phi.get_gradient(q);
2225 *   const auto& u_int = phi.get_value(q);
2226 *   const auto& lambda = std::abs(scalar_product(phi_old_extr.get_value(q), n_plus));
2227 *  
2228 *   const auto& point_vectorized = phi.quadrature_point(q);
2229 *   auto u_int_m = u_int;
2230 *   auto grad_u_int_m = grad_u_int;
2231 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
2232 *   Point<dim> point;
2233 *   for(unsigned int d = 0; d < dim; ++d)
2234 *   point[d] = point_vectorized[d][v];
2235 *  
2236 *   u_int_m[1][v] = -u_int_m[1][v];
2237 *  
2238 *   grad_u_int_m[0][0][v] = -grad_u_int_m[0][0][v];
2239 *   grad_u_int_m[0][1][v] = -grad_u_int_m[0][1][v];
2240 *   }
2241 *  
2242 *   phi.submit_value(a22/Re*(-(0.5*(grad_u_int + grad_u_int_m))*n_plus + coef_jump*(u_int - u_int_m)) +
2243 *   a22*outer_product(0.5*(u_int + u_int_m), phi_old_extr.get_value(q))*n_plus +
2244 *   a22*0.5*lambda*(u_int - u_int_m), q);
2245 *   phi.submit_normal_derivative(-theta_v*a22/Re*(u_int - u_int_m), q);
2246 *   }
2247 *   phi.integrate(true, true);
2248 *   diagonal[i] = phi.get_dof_value(i);
2249 *   }
2250 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2251 *   phi.submit_dof_value(diagonal[i], i);
2252 *   phi.distribute_local_to_global(dst);
2253 *   }
2254 *   }
2255 *   }
2256 *   else {
2258 *   phi_extr(data, true, 0);
2259 *  
2262 *   for(unsigned int d = 0; d < dim; ++d)
2263 *   tmp[d] = make_vectorized_array<Number>(1.0);
2264 *  
2265 *   /*--- Loop over all faces in the range ---*/
2266 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2267 *   phi_extr.reinit(face);
2268 *   phi_extr.gather_evaluate(u_extr, true, false);
2269 *   phi.reinit(face);
2270 *  
2271 *   const auto boundary_id = data.get_boundary_id(face);
2272 *   const auto coef_jump = C_u*std::abs((phi.get_normal_vector(0) * phi.inverse_jacobian(0))[dim - 1]);
2273 *  
2274 *   if(boundary_id != 1) {
2275 *   const double coef_trasp = 0.0;
2276 *  
2277 *   /*--- Loop over all dofs ---*/
2278 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2279 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2280 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2281 *   phi.submit_dof_value(tmp, i);
2282 *   phi.evaluate(true, true);
2283 *  
2284 *   /*--- Loop over quadrature points to compute the integral ---*/
2285 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2286 *   const auto& n_plus = phi.get_normal_vector(q);
2287 *   const auto& grad_u = phi.get_gradient(q);
2288 *   const auto& u = phi.get_value(q);
2289 *   const auto& tensor_product_u = outer_product(phi.get_value(q), phi_extr.get_value(q));
2290 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
2291 *  
2292 *   phi.submit_value(a33/Re*(-grad_u*n_plus + 2.0*coef_jump*u) +
2293 *   a33*coef_trasp*tensor_product_u*n_plus + a33*lambda*u, q);
2294 *   phi.submit_normal_derivative(-theta_v*a33/Re*u, q);
2295 *   }
2296 *   phi.integrate(true, true);
2297 *   diagonal[i] = phi.get_dof_value(i);
2298 *   }
2299 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2300 *   phi.submit_dof_value(diagonal[i], i);
2301 *   phi.distribute_local_to_global(dst);
2302 *   }
2303 *   else {
2304 *   /*--- Loop over all dofs ---*/
2305 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2306 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2307 *   phi.submit_dof_value(Tensor<1, dim, VectorizedArray<Number>>(), j);
2308 *   phi.submit_dof_value(tmp, i);
2309 *   phi.evaluate(true, true);
2310 *  
2311 *   /*--- Loop over quadrature points to compute the integral ---*/
2312 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2313 *   const auto& n_plus = phi.get_normal_vector(q);
2314 *   const auto& grad_u = phi.get_gradient(q);
2315 *   const auto& u = phi.get_value(q);
2316 *   const auto& lambda = std::abs(scalar_product(phi_extr.get_value(q), n_plus));
2317 *  
2318 *   const auto& point_vectorized = phi.quadrature_point(q);
2319 *   auto u_m = u;
2320 *   auto grad_u_m = grad_u;
2321 *   for(unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v) {
2322 *   Point<dim> point;
2323 *   for(unsigned int d = 0; d < dim; ++d)
2324 *   point[d] = point_vectorized[d][v];
2325 *  
2326 *   u_m[1][v] = -u_m[1][v];
2327 *  
2328 *   grad_u_m[0][0][v] = -grad_u_m[0][0][v];
2329 *   grad_u_m[0][1][v] = -grad_u_m[0][1][v];
2330 *   }
2331 *  
2332 *   phi.submit_value(a33/Re*(-(0.5*(grad_u + grad_u_m))*n_plus + coef_jump*(u - u_m)) +
2333 *   a33*outer_product(0.5*(u + u_m), phi_extr.get_value(q))*n_plus +
2334 *   a33*0.5*lambda*(u - u_m), q);
2335 *   phi.submit_normal_derivative(-theta_v*a33/Re*(u - u_m), q);
2336 *   }
2337 *   phi.integrate(true, true);
2338 *   diagonal[i] = phi.get_dof_value(i);
2339 *   }
2340 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2341 *   phi.submit_dof_value(diagonal[i], i);
2342 *   phi.distribute_local_to_global(dst);
2343 *   }
2344 *   }
2345 *   }
2346 *   }
2347 *  
2348 *  
2349 * @endcode
2350 *
2351 * Now we consider the pressure related bilinear forms. We first assemble diagonal cell term for the pressure
2352 *
2353
2354 *
2355 *
2356 * @code
2357 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
2358 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2359 *   assemble_diagonal_cell_term_pressure(const MatrixFree<dim, Number>& data,
2360 *   Vec& dst,
2361 *   const unsigned int& ,
2362 *   const std::pair<unsigned int, unsigned int>& cell_range) const {
2364 *  
2365 *   AlignedVector<VectorizedArray<Number>> diagonal(phi.dofs_per_component); /*--- Here we are using dofs_per_component but
2366 *   it coincides with dofs_per_cell since it is
2367 *   scalar finite element space ---*/
2368 *  
2369 *   const double coeff = (TR_BDF2_stage == 1) ? 1e6*gamma*dt*gamma*dt : 1e6*(1.0 - gamma)*dt*(1.0 - gamma)*dt;
2370 *  
2371 *   /*--- Loop over all cells in the range ---*/
2372 *   for(unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) {
2373 *   phi.reinit(cell);
2374 *  
2375 *   /*--- Loop over all dofs ---*/
2376 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2377 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2378 *   phi.submit_dof_value(VectorizedArray<Number>(), j); /*--- We set all dofs to zero ---*/
2379 *   phi.submit_dof_value(make_vectorized_array<Number>(1.0), i); /*--- Now we set the current one to 1; since it is scalar,
2380 *   we can directly use 'make_vectorized_array' without
2381 *   relying on 'Tensor' ---*/
2383 *  
2384 *   /*--- Loop over quadrature points ---*/
2385 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2386 *   phi.submit_value(1.0/coeff*phi.get_value(q), q);
2387 *   phi.submit_gradient(phi.get_gradient(q), q);
2388 *   }
2390 *   diagonal[i] = phi.get_dof_value(i);
2391 *   }
2392 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2393 *   phi.submit_dof_value(diagonal[i], i);
2394 *  
2395 *   phi.distribute_local_to_global(dst);
2396 *   }
2397 *   }
2398 *  
2399 *  
2400 * @endcode
2401 *
2402 * The following function assembles diagonal face term for the pressure
2403 *
2404
2405 *
2406 *
2407 * @code
2408 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
2409 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2410 *   assemble_diagonal_face_term_pressure(const MatrixFree<dim, Number>& data,
2411 *   Vec& dst,
2412 *   const unsigned int& ,
2413 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2415 *   phi_m(data, false, 1, 1);
2416 *  
2417 *   AssertDimension(phi_p.dofs_per_component, phi_m.dofs_per_component);
2418 *   AlignedVector<VectorizedArray<Number>> diagonal_p(phi_p.dofs_per_component),
2419 *   diagonal_m(phi_m.dofs_per_component); /*--- Again, we just assert for safety that dimension
2420 *   match, in the sense that we have selected
2421 *   the proper space ---*/
2422 *  
2423 *   /*--- Loop over all faces ---*/
2424 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2425 *   phi_p.reinit(face);
2426 *   phi_m.reinit(face);
2427 *  
2428 *   const auto coef_jump = C_p*0.5*(std::abs((phi_p.get_normal_vector(0)*phi_p.inverse_jacobian(0))[dim - 1]) +
2429 *   std::abs((phi_m.get_normal_vector(0)*phi_m.inverse_jacobian(0))[dim - 1]));
2430 *  
2431 *   /*--- Loop over all dofs ---*/
2432 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2433 *   for(unsigned int j = 0; j < phi_p.dofs_per_component; ++j) {
2434 *   phi_p.submit_dof_value(VectorizedArray<Number>(), j);
2435 *   phi_m.submit_dof_value(VectorizedArray<Number>(), j);
2436 *   }
2437 *   phi_p.submit_dof_value(make_vectorized_array<Number>(1.0), i);
2438 *   phi_m.submit_dof_value(make_vectorized_array<Number>(1.0), i);
2441 *  
2442 *   /*--- Loop over all quadrature points to compute the integral ---*/
2443 *   for(unsigned int q = 0; q < phi_p.n_q_points; ++q) {
2444 *   const auto& n_plus = phi_p.get_normal_vector(q);
2445 *  
2446 *   const auto& avg_grad_pres = 0.5*(phi_p.get_gradient(q) + phi_m.get_gradient(q));
2447 *   const auto& jump_pres = phi_p.get_value(q) - phi_m.get_value(q);
2448 *  
2449 *   phi_p.submit_value(-scalar_product(avg_grad_pres, n_plus) + coef_jump*jump_pres, q);
2450 *   phi_m.submit_value(scalar_product(avg_grad_pres, n_plus) - coef_jump*jump_pres, q);
2451 *   phi_p.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
2452 *   phi_m.submit_gradient(-theta_p*0.5*jump_pres*n_plus, q);
2453 *   }
2455 *   diagonal_p[i] = phi_p.get_dof_value(i);
2457 *   diagonal_m[i] = phi_m.get_dof_value(i);
2458 *   }
2459 *   for(unsigned int i = 0; i < phi_p.dofs_per_component; ++i) {
2460 *   phi_p.submit_dof_value(diagonal_p[i], i);
2461 *   phi_m.submit_dof_value(diagonal_m[i], i);
2462 *   }
2463 *   phi_p.distribute_local_to_global(dst);
2464 *   phi_m.distribute_local_to_global(dst);
2465 *   }
2466 *   }
2467 *  
2468 *  
2469 * @endcode
2470 *
2471 * Eventually, we assemble diagonal boundary term for the pressure
2472 *
2473
2474 *
2475 *
2476 * @code
2477 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
2478 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2479 *   assemble_diagonal_boundary_term_pressure(const MatrixFree<dim, Number>& data,
2480 *   Vec& dst,
2481 *   const unsigned int& ,
2482 *   const std::pair<unsigned int, unsigned int>& face_range) const {
2484 *  
2485 *   AlignedVector<VectorizedArray<Number>> diagonal(phi.dofs_per_component);
2486 *  
2487 *   for(unsigned int face = face_range.first; face < face_range.second; ++face) {
2488 *   const auto boundary_id = data.get_boundary_id(face);
2489 *  
2490 *   if(boundary_id == 1) {
2491 *   phi.reinit(face);
2492 *  
2493 *   const auto coef_jump = C_p*std::abs((phi.get_normal_vector(0)*phi.inverse_jacobian(0))[dim - 1]);
2494 *  
2495 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i) {
2496 *   for(unsigned int j = 0; j < phi.dofs_per_component; ++j)
2497 *   phi.submit_dof_value(VectorizedArray<Number>(), j);
2498 *   phi.submit_dof_value(make_vectorized_array<Number>(1.0), i);
2500 *  
2501 *   for(unsigned int q = 0; q < phi.n_q_points; ++q) {
2502 *   const auto& n_plus = phi.get_normal_vector(q);
2503 *  
2504 *   const auto& grad_pres = phi.get_gradient(q);
2505 *   const auto& pres = phi.get_value(q);
2506 *  
2507 *   phi.submit_value(-scalar_product(grad_pres, n_plus) + 2.0*coef_jump*pres , q);
2508 *   phi.submit_normal_derivative(-theta_p*pres, q);
2509 *   }
2511 *   diagonal[i] = phi.get_dof_value(i);
2512 *   }
2513 *   for(unsigned int i = 0; i < phi.dofs_per_component; ++i)
2514 *   phi.submit_dof_value(diagonal[i], i);
2515 *   phi.distribute_local_to_global(dst);
2516 *   }
2517 *   }
2518 *   }
2519 *  
2520 *  
2521 * @endcode
2522 *
2523 * Put together all previous steps. We create a dummy auxliary vector that serves for the src input argument in
2524 * the previous functions that as we have seen before is unused. Then everything is done by the 'loop' function
2525 * and it is saved in the field 'inverse_diagonal_entries' already present in the base class. Anyway since there is
2526 * only one field, we need to resize properly depending on whether we are considering the velocity or the pressure.
2527 *
2528
2529 *
2530 *
2531 * @code
2532 *   template<int dim, int fe_degree_p, int fe_degree_v, int n_q_points_1d_p, int n_q_points_1d_v, typename Vec>
2533 *   void NavierStokesProjectionOperator<dim, fe_degree_p, fe_degree_v, n_q_points_1d_p, n_q_points_1d_v, Vec>::
2534 *   compute_diagonal() {
2535 *   Assert(NS_stage == 1 || NS_stage == 2, ExcInternalError());
2536 *  
2537 *   this->inverse_diagonal_entries.reset(new DiagonalMatrix<Vec>());
2538 *   auto& inverse_diagonal = this->inverse_diagonal_entries->get_vector();
2539 *  
2540 *   if(NS_stage == 1) {
2541 *   ::MatrixFreeTools::compute_diagonal<dim, Number, VectorizedArray<Number>>
2542 *   (*(this->data),
2543 *   inverse_diagonal,
2544 *   [&](const auto& data, auto& dst, const auto& src, const auto& cell_range) {
2545 *   (this->assemble_diagonal_cell_term_velocity)(data, dst, src, cell_range);
2546 *   },
2547 *   [&](const auto& data, auto& dst, const auto& src, const auto& face_range) {
2548 *   (this->assemble_diagonal_face_term_velocity)(data, dst, src, face_range);
2549 *   },
2550 *   [&](const auto& data, auto& dst, const auto& src, const auto& boundary_range) {
2551 *   (this->assemble_diagonal_boundary_term_velocity)(data, dst, src, boundary_range);
2552 *   },
2553 *   0);
2554 *   }
2555 *   else if(NS_stage == 2) {
2556 *   ::MatrixFreeTools::compute_diagonal<dim, Number, VectorizedArray<Number>>
2557 *   (*(this->data),
2558 *   inverse_diagonal,
2559 *   [&](const auto& data, auto& dst, const auto& src, const auto& cell_range) {
2560 *   (this->assemble_diagonal_cell_term_pressure)(data, dst, src, cell_range);
2561 *   },
2562 *   [&](const auto& data, auto& dst, const auto& src, const auto& face_range) {
2563 *   (this->assemble_diagonal_face_term_pressure)(data, dst, src, face_range);
2564 *   },
2565 *   [&](const auto& data, auto& dst, const auto& src, const auto& boundary_range) {
2566 *   (this->assemble_diagonal_boundary_term_pressure)(data, dst, src, boundary_range);
2567 *   },
2568 *   1);
2569 *   }
2570 *  
2571 *   for(unsigned int i = 0; i < inverse_diagonal.locally_owned_size(); ++i) {
2572 *   Assert(inverse_diagonal.local_element(i) != 0.0,
2573 *   ExcMessage("No diagonal entry in a definite operator should be zero"));
2574 *   inverse_diagonal.local_element(i) = 1.0/inverse_diagonal.local_element(i);
2575 *   }
2576 *   }
2577 *  
2578 *  
2579 * @endcode
2580 *
2581 *
2582 * <a name=""></a>
2583 * @sect{The <code>NavierStokesProjection</code> class}
2584 *
2585
2586 *
2587 * Now we are ready for the main class of the program. It implements the calls to the various steps
2588 * of the projection method for Navier-Stokes equations.
2589 *
2590
2591 *
2592 *
2593 * @code
2594 *   template<int dim>
2595 *   class NavierStokesProjection {
2596 *   public:
2597 *   NavierStokesProjection(RunTimeParameters::Data_Storage& data);
2598 *  
2599 *   void run(const bool verbose = false, const unsigned int output_interval = 10);
2600 *  
2601 *   protected:
2602 *   const double t_0;
2603 *   const double T;
2604 *   const double gamma; //--- TR-BDF2 parameter
2605 *   unsigned int TR_BDF2_stage; //--- Flag to check at which current stage of TR-BDF2 are
2606 *   const double Re;
2607 *   double dt;
2608 *  
2609 *   EquationData::Velocity<dim> vel_init;
2610 *   EquationData::Pressure<dim> pres_init; /*--- Instance of 'Velocity' and 'Pressure' classes to initialize. ---*/
2611 *  
2613 *  
2614 *   /*--- Finite Element spaces ---*/
2615 *   FESystem<dim> fe_velocity;
2616 *   FESystem<dim> fe_pressure;
2617 *  
2618 *   /*--- Handler for dofs ---*/
2619 *   DoFHandler<dim> dof_handler_velocity;
2620 *   DoFHandler<dim> dof_handler_pressure;
2621 *  
2622 *   /*--- Quadrature formulas for velocity and pressure, respectively ---*/
2623 *   QGauss<dim> quadrature_pressure;
2624 *   QGauss<dim> quadrature_velocity;
2625 *  
2626 *   /*--- Now we define all the vectors for the solution. We start from the pressure
2627 *   with p^n, p^(n+gamma) and a vector for rhs ---*/
2631 *  
2632 *   /*--- Next, we move to the velocity, with u^n, u^(n-1), u^(n+gamma/2),
2633 *   u^(n+gamma) and other two auxiliary vectors as well as the rhs ---*/
2642 *  
2643 *   Vector<double> Linfty_error_per_cell_vel;
2644 *  
2645 *   DeclException2(ExcInvalidTimeStep,
2646 *   double,
2647 *   double,
2648 *   << " The time step " << arg1 << " is out of range."
2649 *   << std::endl
2650 *   << " The permitted range is (0," << arg2 << "]");
2651 *  
2652 *   void create_triangulation(const unsigned int n_refines);
2653 *  
2654 *   void setup_dofs();
2655 *  
2656 *   void initialize();
2657 *  
2658 *   void interpolate_velocity();
2659 *  
2660 *   void diffusion_step();
2661 *  
2662 *   void projection_step();
2663 *  
2664 *   void project_grad(const unsigned int flag);
2665 *  
2666 *   double get_maximal_velocity();
2667 *  
2668 *   double get_maximal_difference_velocity();
2669 *  
2670 *   void output_results(const unsigned int step);
2671 *  
2672 *   void refine_mesh();
2673 *  
2674 *   void interpolate_max_res(const unsigned int level);
2675 *  
2676 *   void save_max_res();
2677 *  
2678 *   private:
2679 *   void compute_lift_and_drag();
2680 *  
2681 *   /*--- Technical member to handle the various steps ---*/
2682 *   std::shared_ptr<MatrixFree<dim, double>> matrix_free_storage;
2683 *  
2684 *   /*--- Now we need an instance of the class implemented before with the weak form ---*/
2685 *   NavierStokesProjectionOperator<dim, EquationData::degree_p, EquationData::degree_p + 1,
2686 *   EquationData::degree_p + 1, EquationData::degree_p + 2,
2687 *   LinearAlgebra::distributed::Vector<double>> navier_stokes_matrix;
2688 *  
2689 *   /*--- This is an instance for geometric multigrid preconditioner ---*/
2690 *   MGLevelObject<NavierStokesProjectionOperator<dim, EquationData::degree_p, EquationData::degree_p + 1,
2691 *   EquationData::degree_p + 1, EquationData::degree_p + 2,
2693 *  
2694 *   /*--- Here we define two 'AffineConstraints' instance, one for each finite element space.
2695 *   This is just a technical issue, due to MatrixFree requirements. In general
2696 *   this class is used to impose boundary conditions (or any kind of constraints), but in this case, since
2697 *   we are using a weak imposition of bcs, everything is already in the weak forms and so these instances
2698 *   will be default constructed ---*/
2699 *   AffineConstraints<double> constraints_velocity,
2700 *   constraints_pressure;
2701 *  
2702 *   /*--- Now a bunch of variables handled by 'ParamHandler' introduced at the beginning of the code ---*/
2703 *   unsigned int max_its;
2704 *   double eps;
2705 *  
2706 *   unsigned int max_loc_refinements;
2707 *   unsigned int min_loc_refinements;
2708 *   unsigned int refinement_iterations;
2709 *  
2710 *   std::string saving_dir;
2711 *  
2712 *   /*--- Finally, some output related streams ---*/
2713 *   ConditionalOStream pcout;
2714 *  
2715 *   std::ofstream time_out;
2716 *   ConditionalOStream ptime_out;
2717 *   TimerOutput time_table;
2718 *  
2719 *   std::ofstream output_n_dofs_velocity;
2720 *   std::ofstream output_n_dofs_pressure;
2721 *  
2722 *   std::ofstream output_lift;
2723 *   std::ofstream output_drag;
2724 *   };
2725 *  
2726 *  
2727 * @endcode
2728 *
2729 * In the constructor, we just read all the data from the
2730 * <code>Data_Storage</code> object that is passed as an argument, verify that
2731 * the data we read are reasonable and, finally, create the triangulation and
2732 * load the initial data.
2733 *
2734
2735 *
2736 *
2737 * @code
2738 *   template<int dim>
2739 *   NavierStokesProjection<dim>::NavierStokesProjection(RunTimeParameters::Data_Storage& data):
2740 *   t_0(data.initial_time),
2741 *   T(data.final_time),
2742 *   gamma(2.0 - std::sqrt(2.0)), //--- Save also in the NavierStokes class the TR-BDF2 parameter value
2743 *   TR_BDF2_stage(1), //--- Initialize the flag for the TR_BDF2 stage
2744 *   Re(data.Reynolds),
2745 *   dt(data.dt),
2746 *   vel_init(data.initial_time),
2747 *   pres_init(data.initial_time),
2748 *   triangulation(MPI_COMM_WORLD, parallel::distributed::Triangulation<dim>::limit_level_difference_at_vertices,
2749 *   parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy),
2750 *   fe_velocity(FE_DGQ<dim>(EquationData::degree_p + 1), dim),
2751 *   fe_pressure(FE_DGQ<dim>(EquationData::degree_p), 1),
2752 *   dof_handler_velocity(triangulation),
2753 *   dof_handler_pressure(triangulation),
2754 *   quadrature_pressure(EquationData::degree_p + 1),
2755 *   quadrature_velocity(EquationData::degree_p + 2),
2756 *   navier_stokes_matrix(data),
2757 *   max_its(data.max_iterations),
2758 *   eps(data.eps),
2759 *   max_loc_refinements(data.max_loc_refinements),
2760 *   min_loc_refinements(data.min_loc_refinements),
2761 *   refinement_iterations(data.refinement_iterations),
2762 *   saving_dir(data.dir),
2763 *   pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
2764 *   time_out("./" + data.dir + "/time_analysis_" +
2765 *   Utilities::int_to_string(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)) + "proc.dat"),
2766 *   ptime_out(time_out, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
2767 *   time_table(ptime_out, TimerOutput::summary, TimerOutput::cpu_and_wall_times),
2768 *   output_n_dofs_velocity("./" + data.dir + "/n_dofs_velocity.dat", std::ofstream::out),
2769 *   output_n_dofs_pressure("./" + data.dir + "/n_dofs_pressure.dat", std::ofstream::out),
2770 *   output_lift("./" + data.dir + "/lift.dat", std::ofstream::out),
2771 *   output_drag("./" + data.dir + "/drag.dat", std::ofstream::out) {
2772 *   if(EquationData::degree_p < 1) {
2773 *   pcout
2774 *   << " WARNING: The chosen pair of finite element spaces is not stable."
2775 *   << std::endl
2776 *   << " The obtained results will be nonsense" << std::endl;
2777 *   }
2778 *  
2779 *   AssertThrow(!((dt <= 0.0) || (dt > 0.5*T)), ExcInvalidTimeStep(dt, 0.5*T));
2780 *  
2781 *   matrix_free_storage = std::make_shared<MatrixFree<dim, double>>();
2782 *  
2783 *   create_triangulation(data.n_refines);
2784 *   setup_dofs();
2785 *   initialize();
2786 *   }
2787 *  
2788 *  
2789 * @endcode
2790 *
2791 * The method that creates the triangulation and refines it the needed number
2792 * of times.
2793 *
2794
2795 *
2796 *
2797 * @code
2798 *   template<int dim>
2799 *   void NavierStokesProjection<dim>::create_triangulation(const unsigned int n_refines) {
2800 *   TimerOutput::Scope t(time_table, "Create triangulation");
2801 *  
2802 *   GridGenerator::plate_with_a_hole(triangulation, 0.5, 1.0, 1.0, 1.1, 1.0, 19.0, Point<2>(2.0, 2.0), 0, 1, 1.0, 2, true);
2803 *   /*--- We strongly advice to check the documentation to verify the meaning of all input parameters. ---*/
2804 *  
2805 *   pcout << "Number of refines = " << n_refines << std::endl;
2806 *   triangulation.refine_global(n_refines);
2807 *   }
2808 *  
2809 *  
2810 * @endcode
2811 *
2812 * After creating the triangulation, it creates the mesh dependent
2813 * data, i.e. it distributes degrees of freedom, and
2814 * initializes the vectors that we will use.
2815 *
2816
2817 *
2818 *
2819 * @code
2820 *   template<int dim>
2821 *   void NavierStokesProjection<dim>::setup_dofs() {
2822 *   pcout << "Number of active cells: " << triangulation.n_global_active_cells() << std::endl;
2823 *   pcout << "Number of levels: " << triangulation.n_global_levels() << std::endl;
2824 *  
2825 *   /*--- Distribute dofs and prepare for multigrid ---*/
2826 *   dof_handler_velocity.distribute_dofs(fe_velocity);
2827 *   dof_handler_pressure.distribute_dofs(fe_pressure);
2828 *  
2829 *   pcout << "dim (X_h) = " << dof_handler_velocity.n_dofs()
2830 *   << std::endl
2831 *   << "dim (M_h) = " << dof_handler_pressure.n_dofs()
2832 *   << std::endl
2833 *   << "Re = " << Re << std::endl
2834 *   << std::endl;
2835 *  
2836 *   if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
2837 *   output_n_dofs_velocity << dof_handler_velocity.n_dofs() << std::endl;
2838 *   output_n_dofs_pressure << dof_handler_pressure.n_dofs() << std::endl;
2839 *   }
2840 *  
2841 *   typename MatrixFree<dim, double>::AdditionalData additional_data;
2842 *   additional_data.mapping_update_flags = (update_gradients | update_JxW_values |
2844 *   additional_data.mapping_update_flags_inner_faces = (update_gradients | update_JxW_values | update_quadrature_points |
2846 *   additional_data.mapping_update_flags_boundary_faces = (update_gradients | update_JxW_values | update_quadrature_points |
2848 *   additional_data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none;
2849 *  
2850 *   std::vector<const DoFHandler<dim>*> dof_handlers; /*--- Vector of dof_handlers to feed the 'MatrixFree'. Here the order
2851 *   counts and enters into the game as parameter of FEEvaluation and
2852 *   FEFaceEvaluation in the previous class ---*/
2853 *   dof_handlers.push_back(&dof_handler_velocity);
2854 *   dof_handlers.push_back(&dof_handler_pressure);
2855 *  
2856 *   constraints_velocity.clear();
2857 *   constraints_velocity.close();
2858 *   constraints_pressure.clear();
2859 *   constraints_pressure.close();
2860 *   std::vector<const AffineConstraints<double>*> constraints;
2861 *   constraints.push_back(&constraints_velocity);
2862 *   constraints.push_back(&constraints_pressure);
2863 *  
2864 *   std::vector<QGauss<1>> quadratures; /*--- We cannot directly use 'quadrature_velocity' and 'quadrature_pressure',
2865 *   because the 'MatrixFree' structure wants a quadrature formula for 1D
2866 *   (this is way the template parameter of the previous class was called 'n_q_points_1d_p'
2867 *   and 'n_q_points_1d_v' and the reason of '1' as QGauss template parameter). ---*/
2868 *   quadratures.push_back(QGauss<1>(EquationData::degree_p + 2));
2869 *   quadratures.push_back(QGauss<1>(EquationData::degree_p + 1));
2870 *  
2871 *   /*--- Initialize the matrix-free structure and size properly the vectors. Here again the
2872 *   second input argument of the 'initialize_dof_vector' method depends on the order of 'dof_handlers' ---*/
2873 *   matrix_free_storage->reinit(MappingQ1<dim>(),dof_handlers, constraints, quadratures, additional_data);
2874 *   matrix_free_storage->initialize_dof_vector(u_star, 0);
2875 *   matrix_free_storage->initialize_dof_vector(rhs_u, 0);
2876 *   matrix_free_storage->initialize_dof_vector(u_n, 0);
2877 *   matrix_free_storage->initialize_dof_vector(u_extr, 0);
2878 *   matrix_free_storage->initialize_dof_vector(u_n_minus_1, 0);
2879 *   matrix_free_storage->initialize_dof_vector(u_n_gamma, 0);
2880 *   matrix_free_storage->initialize_dof_vector(u_tmp, 0);
2881 *   matrix_free_storage->initialize_dof_vector(grad_pres_int, 0);
2882 *  
2883 *   matrix_free_storage->initialize_dof_vector(pres_int, 1);
2884 *   matrix_free_storage->initialize_dof_vector(pres_n, 1);
2885 *   matrix_free_storage->initialize_dof_vector(rhs_p, 1);
2886 *  
2887 *   /*--- Initialize the multigrid structure. We dedicate ad hoc 'dof_handlers_mg' and 'constraints_mg' because
2888 *   we use float as type. Moreover we can initialize already with the index of the finite element of the pressure;
2889 *   anyway we need by requirement to declare also structures for the velocity for coherence (basically because
2890 *   the index of finite element space has to be the same, so the pressure has to be the second).---*/
2891 *   mg_matrices.clear_elements();
2892 *   dof_handler_velocity.distribute_mg_dofs();
2893 *   dof_handler_pressure.distribute_mg_dofs();
2894 *  
2895 *   const unsigned int nlevels = triangulation.n_global_levels();
2896 *   mg_matrices.resize(0, nlevels - 1);
2897 *   for(unsigned int level = 0; level < nlevels; ++level) {
2898 *   typename MatrixFree<dim, float>::AdditionalData additional_data_mg;
2900 *   additional_data_mg.mapping_update_flags = (update_gradients | update_JxW_values);
2901 *   additional_data_mg.mapping_update_flags_inner_faces = (update_gradients | update_JxW_values);
2902 *   additional_data_mg.mapping_update_flags_boundary_faces = (update_gradients | update_JxW_values);
2903 *   additional_data_mg.mg_level = level;
2904 *  
2905 *   std::vector<const DoFHandler<dim>*> dof_handlers_mg;
2906 *   dof_handlers_mg.push_back(&dof_handler_velocity);
2907 *   dof_handlers_mg.push_back(&dof_handler_pressure);
2908 *   std::vector<const AffineConstraints<float>*> constraints_mg;
2909 *   AffineConstraints<float> constraints_velocity_mg;
2910 *   constraints_velocity_mg.clear();
2911 *   constraints_velocity_mg.close();
2912 *   constraints_mg.push_back(&constraints_velocity_mg);
2913 *   AffineConstraints<float> constraints_pressure_mg;
2914 *   constraints_pressure_mg.clear();
2915 *   constraints_pressure_mg.close();
2916 *   constraints_mg.push_back(&constraints_pressure_mg);
2917 *  
2918 *   std::shared_ptr<MatrixFree<dim, float>> mg_mf_storage_level(new MatrixFree<dim, float>());
2919 *   mg_mf_storage_level->reinit(MappingQ1<dim>(),dof_handlers_mg, constraints_mg, quadratures, additional_data_mg);
2920 *   const std::vector<unsigned int> tmp = {1};
2921 *   mg_matrices[level].initialize(mg_mf_storage_level, tmp, tmp);
2922 *   mg_matrices[level].set_dt(dt);
2923 *   mg_matrices[level].set_NS_stage(2);
2924 *   }
2925 *  
2926 *   Linfty_error_per_cell_vel.reinit(triangulation.n_active_cells());
2927 *   }
2928 *  
2929 *  
2930 * @endcode
2931 *
2932 * This method loads the initial data. It simply uses the class <code>Pressure</code> instance for the pressure
2933 * and the class <code>Velocity</code> instance for the velocity.
2934 *
2935
2936 *
2937 *
2938 * @code
2939 *   template<int dim>
2940 *   void NavierStokesProjection<dim>::initialize() {
2941 *   TimerOutput::Scope t(time_table, "Initialize pressure and velocity");
2942 *  
2943 *   VectorTools::interpolate(dof_handler_pressure, pres_init, pres_n);
2944 *  
2945 *   VectorTools::interpolate(dof_handler_velocity, vel_init, u_n_minus_1);
2946 *   VectorTools::interpolate(dof_handler_velocity, vel_init, u_n);
2947 *   }
2948 *  
2949 *  
2950 * @endcode
2951 *
2952 * This function computes the extrapolated velocity to be used in the momentum predictor
2953 *
2954
2955 *
2956 *
2957 * @code
2958 *   template<int dim>
2959 *   void NavierStokesProjection<dim>::interpolate_velocity() {
2960 *   TimerOutput::Scope t(time_table, "Interpolate velocity");
2961 *  
2962 * @endcode
2963 *
2964 * --- TR-BDF2 first step
2965 *
2966 * @code
2967 *   if(TR_BDF2_stage == 1) {
2968 *   u_extr.equ(1.0 + gamma/(2.0*(1.0 - gamma)), u_n);
2969 *   u_tmp.equ(gamma/(2.0*(1.0 - gamma)), u_n_minus_1);
2970 *   u_extr -= u_tmp;
2971 *   }
2972 * @endcode
2973 *
2974 * --- TR-BDF2 second step
2975 *
2976 * @code
2977 *   else {
2978 *   u_extr.equ(1.0 + (1.0 - gamma)/gamma, u_n_gamma);
2979 *   u_tmp.equ((1.0 - gamma)/gamma, u_n);
2980 *   u_extr -= u_tmp;
2981 *   }
2982 *   }
2983 *  
2984 *  
2985 * @endcode
2986 *
2987 * We are finally ready to solve the diffusion step.
2988 *
2989
2990 *
2991 *
2992 * @code
2993 *   template<int dim>
2994 *   void NavierStokesProjection<dim>::diffusion_step() {
2995 *   TimerOutput::Scope t(time_table, "Diffusion step");
2996 *  
2997 *   /*--- We first speicify that we want to deal with velocity dof_handler (index 0, since it is the first one
2998 *   in the 'dof_handlers' vector) ---*/
2999 *   const std::vector<unsigned int> tmp = {0};
3000 *   navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
3001 *  
3002 *   /*--- Next, we specify at we are at stage 1, namely the diffusion step ---*/
3003 *   navier_stokes_matrix.set_NS_stage(1);
3004 *  
3005 *   /*--- Now, we compute the right-hand side and we set the convective velocity. The necessity of 'set_u_extr' is
3006 *   that this quantity is required in the bilinear forms and we can't use a vector of src like on the right-hand side,
3007 *   so it has to be available ---*/
3008 *   if(TR_BDF2_stage == 1) {
3009 *   navier_stokes_matrix.vmult_rhs_velocity(rhs_u, {u_n, u_extr, pres_n});
3010 *   navier_stokes_matrix.set_u_extr(u_extr);
3011 *   u_star = u_extr;
3012 *   }
3013 *   else {
3014 *   navier_stokes_matrix.vmult_rhs_velocity(rhs_u, {u_n, u_n_gamma, pres_int, u_extr});
3015 *   navier_stokes_matrix.set_u_extr(u_extr);
3016 *   u_star = u_extr;
3017 *   }
3018 *  
3019 *   /*--- Build the linear solver; in this case we specifiy the maximum number of iterations and residual ---*/
3020 *   SolverControl solver_control(max_its, eps*rhs_u.l2_norm());
3022 *  
3023 *   /*--- Build a Jacobi preconditioner and solve ---*/
3024 *   PreconditionJacobi<NavierStokesProjectionOperator<dim,
3025 *   EquationData::degree_p,
3026 *   EquationData::degree_p + 1,
3027 *   EquationData::degree_p + 1,
3028 *   EquationData::degree_p + 2,
3029 *   LinearAlgebra::distributed::Vector<double>>> preconditioner;
3030 *   navier_stokes_matrix.compute_diagonal();
3031 *   preconditioner.initialize(navier_stokes_matrix);
3032 *  
3033 *   gmres.solve(navier_stokes_matrix, u_star, rhs_u, preconditioner);
3034 *   }
3035 *  
3036 *  
3037 * @endcode
3038 *
3039 * Next, we solve the projection step.
3040 *
3041
3042 *
3043 *
3044 * @code
3045 *   template<int dim>
3046 *   void NavierStokesProjection<dim>::projection_step() {
3047 *   TimerOutput::Scope t(time_table, "Projection step pressure");
3048 *  
3049 *   /*--- We start in the same way of 'diffusion_step': we first reinitialize with the index of FE space,
3050 *   we specify that this is the second stage and we compute the right-hand side ---*/
3051 *   const std::vector<unsigned int> tmp = {1};
3052 *   navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
3053 *  
3054 *   navier_stokes_matrix.set_NS_stage(2);
3055 *  
3056 *   if(TR_BDF2_stage == 1)
3057 *   navier_stokes_matrix.vmult_rhs_pressure(rhs_p, {u_star, pres_n});
3058 *   else
3059 *   navier_stokes_matrix.vmult_rhs_pressure(rhs_p, {u_star, pres_int});
3060 *  
3061 *   /*--- Build the linear solver (Conjugate Gradient in this case) ---*/
3062 *   SolverControl solver_control(max_its, eps*rhs_p.l2_norm());
3064 *  
3065 *   /*--- Build the preconditioner (as in @ref step_37 "step-37") ---*/
3066 *   MGTransferMatrixFree<dim, float> mg_transfer;
3067 *   mg_transfer.build(dof_handler_pressure);
3068 *  
3069 *   using SmootherType = PreconditionChebyshev<NavierStokesProjectionOperator<dim,
3070 *   EquationData::degree_p,
3071 *   EquationData::degree_p + 1,
3072 *   EquationData::degree_p + 1,
3073 *   EquationData::degree_p + 2,
3078 *   smoother_data.resize(0, triangulation.n_global_levels() - 1);
3079 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level) {
3080 *   if(level > 0) {
3081 *   smoother_data[level].smoothing_range = 15.0;
3082 *   smoother_data[level].degree = 3;
3083 *   smoother_data[level].eig_cg_n_iterations = 10;
3084 *   }
3085 *   else {
3086 *   smoother_data[0].smoothing_range = 2e-2;
3087 *   smoother_data[0].degree = numbers::invalid_unsigned_int;
3088 *   smoother_data[0].eig_cg_n_iterations = mg_matrices[0].m();
3089 *   }
3090 *   mg_matrices[level].compute_diagonal();
3091 *   smoother_data[level].preconditioner = mg_matrices[level].get_matrix_diagonal_inverse();
3092 *   }
3093 *   mg_smoother.initialize(mg_matrices, smoother_data);
3094 *  
3096 *   SolverCG<LinearAlgebra::distributed::Vector<float>> cg_mg(solver_control);
3099 *   NavierStokesProjectionOperator<dim,
3100 *   EquationData::degree_p,
3101 *   EquationData::degree_p + 1,
3102 *   EquationData::degree_p + 1,
3103 *   EquationData::degree_p + 2,
3105 *   PreconditionIdentity> mg_coarse(cg_mg, mg_matrices[0], identity);
3106 *  
3107 *   mg::Matrix<LinearAlgebra::distributed::Vector<float>> mg_matrix(mg_matrices);
3108 *  
3109 *   Multigrid<LinearAlgebra::distributed::Vector<float>> mg(mg_matrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
3110 *  
3111 *   PreconditionMG<dim,
3113 *   MGTransferMatrixFree<dim, float>> preconditioner(dof_handler_pressure, mg, mg_transfer);
3114 *  
3115 *   /*--- Solve the linear system ---*/
3116 *   if(TR_BDF2_stage == 1) {
3117 *   pres_int = pres_n;
3118 *   cg.solve(navier_stokes_matrix, pres_int, rhs_p, preconditioner);
3119 *   }
3120 *   else {
3121 *   pres_n = pres_int;
3122 *   cg.solve(navier_stokes_matrix, pres_n, rhs_p, preconditioner);
3123 *   }
3124 *   }
3125 *  
3126 *  
3127 * @endcode
3128 *
3129 * This implements the projection step for the gradient of pressure
3130 *
3131
3132 *
3133 *
3134 * @code
3135 *   template<int dim>
3136 *   void NavierStokesProjection<dim>::project_grad(const unsigned int flag) {
3137 *   TimerOutput::Scope t(time_table, "Gradient of pressure projection");
3138 *  
3139 *   /*--- The input parameter flag is used just to specify where we want to save the result ---*/
3140 *   AssertIndexRange(flag, 3);
3141 *   Assert(flag > 0, ExcInternalError());
3142 *  
3143 *   /*--- We need to select the dof handler related to the velocity since the result lives there ---*/
3144 *   const std::vector<unsigned int> tmp = {0};
3145 *   navier_stokes_matrix.initialize(matrix_free_storage, tmp, tmp);
3146 *  
3147 *   if(flag == 1)
3148 *   navier_stokes_matrix.vmult_grad_p_projection(rhs_u, pres_n);
3149 *   else if(flag == 2)
3150 *   navier_stokes_matrix.vmult_grad_p_projection(rhs_u, pres_int);
3151 *  
3152 *   /*--- We conventionally decide that the this corresponds to third stage ---*/
3153 *   navier_stokes_matrix.set_NS_stage(3);
3154 *  
3155 *   /*--- Solve the system ---*/
3156 *   SolverControl solver_control(max_its, 1e-12*rhs_u.l2_norm());
3158 *   cg.solve(navier_stokes_matrix, u_tmp, rhs_u, PreconditionIdentity());
3159 *   }
3160 *  
3161 *  
3162 * @endcode
3163 *
3164 * The following function is used in determining the maximal velocity
3165 * in order to compute the Courant number.
3166 *
3167
3168 *
3169 *
3170 * @code
3171 *   template<int dim>
3172 *   double NavierStokesProjection<dim>::get_maximal_velocity() {
3173 *   return u_n.linfty_norm();
3174 *   }
3175 *  
3176 *  
3177 * @endcode
3178 *
3179 * The following function is used in determining the maximal nodal difference
3180 * between old and current velocity value in order to see if we have reched steady-state.
3181 *
3182
3183 *
3184 *
3185 * @code
3186 *   template<int dim>
3187 *   double NavierStokesProjection<dim>::get_maximal_difference_velocity() {
3188 *   u_tmp = u_n;
3189 *   u_tmp -= u_n_minus_1;
3190 *  
3191 *   return u_tmp.linfty_norm();
3192 *   }
3193 *  
3194 *  
3195 * @endcode
3196 *
3197 * This method plots the current solution. The main difficulty is that we want
3198 * to create a single output file that contains the data for all velocity
3199 * components and the pressure. On the other hand, velocities and the pressure
3200 * live on separate DoFHandler objects, so we need to pay attention when we use
3201 * 'add_data_vector' to select the proper space.
3202 *
3203
3204 *
3205 *
3206 * @code
3207 *   template<int dim>
3208 *   void NavierStokesProjection<dim>::output_results(const unsigned int step) {
3209 *   TimerOutput::Scope t(time_table, "Output results");
3210 *  
3211 *   DataOut<dim> data_out;
3212 *  
3213 *   std::vector<std::string> velocity_names(dim, "v");
3214 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
3215 *   component_interpretation_velocity(dim, DataComponentInterpretation::component_is_part_of_vector);
3216 *   u_n.update_ghost_values();
3217 *   data_out.add_data_vector(dof_handler_velocity, u_n, velocity_names, component_interpretation_velocity);
3218 *   pres_n.update_ghost_values();
3219 *   data_out.add_data_vector(dof_handler_pressure, pres_n, "p", {DataComponentInterpretation::component_is_scalar});
3220 *  
3221 *   std::vector<std::string> velocity_names_old(dim, "v_old");
3222 *   u_n_minus_1.update_ghost_values();
3223 *   data_out.add_data_vector(dof_handler_velocity, u_n_minus_1, velocity_names_old, component_interpretation_velocity);
3224 *  
3225 *   /*--- Here we rely on the postprocessor we have built ---*/
3226 *   PostprocessorVorticity<dim> postprocessor;
3227 *   data_out.add_data_vector(dof_handler_velocity, u_n, postprocessor);
3228 *  
3229 *   data_out.build_patches(MappingQ1<dim>(), 1, DataOut<dim>::curved_inner_cells);
3230 *  
3231 *   const std::string output = "./" + saving_dir + "/solution-" + Utilities::int_to_string(step, 5) + ".vtu";
3232 *   data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD);
3233 *   }
3234 *  
3235 *  
3236 * @endcode
3237 *
3238 *
3239 * <a name=""></a>
3240 * @sect{<code>NavierStokesProjection::compute_lift_and_drag</code>}
3241 *
3242
3243 *
3244 * This routine computes the lift and the drag forces in a non-dimensional framework
3245 * (so basically for the classical coefficients, it is necessary to multiply by a factor 2).
3246 *
3247
3248 *
3249 *
3250 * @code
3251 *   template<int dim>
3252 *   void NavierStokesProjection<dim>::compute_lift_and_drag() {
3253 *   QGauss<dim - 1> face_quadrature_formula(EquationData::degree_p + 2);
3254 *   const int n_q_points = face_quadrature_formula.size();
3255 *  
3256 *   std::vector<double> pressure_values(n_q_points);
3257 *   std::vector<std::vector<Tensor<1, dim>>> velocity_gradients(n_q_points, std::vector<Tensor<1, dim>>(dim));
3258 *  
3259 *   Tensor<1, dim> normal_vector;
3260 *   Tensor<2, dim> fluid_stress;
3261 *   Tensor<2, dim> fluid_pressure;
3262 *   Tensor<1, dim> forces;
3263 *  
3264 *   /*--- We need to compute the integral over the cylinder boundary, so we need to use 'FEFaceValues' instances.
3265 *   For the velocity we need the gradients, for the pressure the values. ---*/
3266 *   FEFaceValues<dim> fe_face_values_velocity(fe_velocity, face_quadrature_formula,
3269 *   FEFaceValues<dim> fe_face_values_pressure(fe_pressure, face_quadrature_formula, update_values);
3270 *  
3271 *   double local_drag = 0.0;
3272 *   double local_lift = 0.0;
3273 *  
3274 *   /*--- We need to perform a unique loop because the whole stress tensor takes into account contributions of
3275 *   velocity and pressure obviously. However, the two dof_handlers are different, so we neede to create an ad-hoc
3276 *   iterator for the pressure that we update manually. It is guaranteed that the cells are visited in the same order
3277 *   (see the documentation) ---*/
3278 *   auto tmp_cell = dof_handler_pressure.begin_active();
3279 *   for(const auto& cell : dof_handler_velocity.active_cell_iterators()) {
3280 *   if(cell->is_locally_owned()) {
3281 *   for(unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell; ++face) {
3282 *   if(cell->face(face)->at_boundary() && cell->face(face)->boundary_id() == 4) {
3283 *   fe_face_values_velocity.reinit(cell, face);
3284 *   fe_face_values_pressure.reinit(tmp_cell, face);
3285 *  
3286 *   fe_face_values_velocity.get_function_gradients(u_n, velocity_gradients); /*--- velocity gradients ---*/
3287 *   fe_face_values_pressure.get_function_values(pres_n, pressure_values); /*--- pressure values ---*/
3288 *  
3289 *   for(int q = 0; q < n_q_points; q++) {
3290 *   normal_vector = -fe_face_values_velocity.normal_vector(q);
3291 *  
3292 *   for(unsigned int d = 0; d < dim; ++ d) {
3293 *   fluid_pressure[d][d] = pressure_values[q];
3294 *   for(unsigned int k = 0; k < dim; ++k)
3295 *   fluid_stress[d][k] = 1.0/Re*velocity_gradients[q][d][k];
3296 *   }
3297 *   fluid_stress = fluid_stress - fluid_pressure;
3298 *  
3299 *   forces = fluid_stress*normal_vector*fe_face_values_velocity.JxW(q);
3300 *  
3301 *   local_drag += forces[0];
3302 *   local_lift += forces[1];
3303 *   }
3304 *   }
3305 *   }
3306 *   }
3307 *   ++tmp_cell;
3308 *   }
3309 *  
3310 *   /*--- At the end, each processor has computed the contribution to the boundary cells it owns and, therefore,
3311 *   we need to sum up all the contributions. ---*/
3312 *   const double lift = Utilities::MPI::sum(local_lift, MPI_COMM_WORLD);
3313 *   const double drag = Utilities::MPI::sum(local_drag, MPI_COMM_WORLD);
3314 *   if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) {
3315 *   output_lift << lift << std::endl;
3316 *   output_drag << drag << std::endl;
3317 *   }
3318 *   }
3319 *  
3320 *  
3321 * @endcode
3322 *
3323 *
3324 * <a name=""></a>
3325 * @sect{ <code>NavierStokesProjection::refine_mesh</code>}
3326 *
3327
3328 *
3329 * After finding a good initial guess on the coarse mesh, we hope to
3330 * decrease the error through refining the mesh. We also need to transfer the current solution to the
3331 * next mesh using the SolutionTransfer class.
3332 *
3333
3334 *
3335 *
3336 * @code
3337 *   template <int dim>
3338 *   void NavierStokesProjection<dim>::refine_mesh() {
3339 *   TimerOutput::Scope t(time_table, "Refine mesh");
3340 *  
3341 *   /*--- We first create a proper vector for computing estimator ---*/
3342 *   IndexSet locally_relevant_dofs;
3343 *   DoFTools::extract_locally_relevant_dofs(dof_handler_velocity, locally_relevant_dofs);
3345 *   tmp_velocity.reinit(dof_handler_velocity.locally_owned_dofs(), locally_relevant_dofs, MPI_COMM_WORLD);
3346 *   tmp_velocity = u_n;
3347 *   tmp_velocity.update_ghost_values();
3348 *  
3349 *   using Iterator = typename DoFHandler<dim>::active_cell_iterator;
3350 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
3351 *  
3352 *   /*--- This is basically the indicator per cell computation (see @ref step_50 "step-50"). Since it is not so complciated
3353 *   we implement it through a lambda expression ---*/
3354 *   const auto cell_worker = [&](const Iterator& cell,
3355 *   ScratchData<dim>& scratch_data,
3356 *   CopyData& copy_data) {
3357 *   FEValues<dim>& fe_values = scratch_data.fe_values; /*--- Here we finally use the 'FEValues' inside ScratchData ---*/
3358 *   fe_values.reinit(cell);
3359 *  
3360 *   /*--- Compute the gradients for all quadrature points ---*/
3361 *   std::vector<std::vector<Tensor<1, dim>>> gradients(fe_values.n_quadrature_points, std::vector<Tensor<1, dim>>(dim));
3362 *   fe_values.get_function_gradients(tmp_velocity, gradients);
3363 *   copy_data.cell_index = cell->active_cell_index();
3364 *   double vorticity_norm_square = 0.0;
3365 *   /*--- Loop over quadrature points and evaluate the integral multiplying the vorticty
3366 *   by the weights and the determinant of the Jacobian (which are included in 'JxW') ---*/
3367 *   for(unsigned k = 0; k < fe_values.n_quadrature_points; ++k) {
3368 *   const double vorticity = gradients[k][1][0] - gradients[k][0][1];
3369 *   vorticity_norm_square += vorticity*vorticity*fe_values.JxW(k);
3370 *   }
3371 *   copy_data.value = cell->diameter()*cell->diameter()*vorticity_norm_square;
3372 *   };
3373 *  
3375 *  
3376 *   const auto copier = [&](const CopyData &copy_data) {
3377 *   if(copy_data.cell_index != numbers::invalid_unsigned_int)
3378 *   estimated_error_per_cell[copy_data.cell_index] += copy_data.value;
3379 *   };
3380 *  
3381 *   /*--- Now everything is 'automagically' handled by 'mesh_loop' ---*/
3382 *   ScratchData<dim> scratch_data(fe_velocity, EquationData::degree_p + 2, cell_flags);
3383 *   CopyData copy_data;
3384 *   MeshWorker::mesh_loop(dof_handler_velocity.begin_active(),
3385 *   dof_handler_velocity.end(),
3386 *   cell_worker,
3387 *   copier,
3388 *   scratch_data,
3389 *   copy_data,
3391 *  
3392 *   /*--- Refine grid. In case the refinement level is above a certain value (or the coarsening level is below)
3393 *   we clear the flags. ---*/
3395 *   for(const auto& cell: triangulation.active_cell_iterators()) {
3396 *   if(cell->refine_flag_set() && static_cast<unsigned int>(cell->level()) == max_loc_refinements)
3397 *   cell->clear_refine_flag();
3398 *   if(cell->coarsen_flag_set() && static_cast<unsigned int>(cell->level()) == min_loc_refinements)
3399 *   cell->clear_coarsen_flag();
3400 *   }
3401 *   triangulation.prepare_coarsening_and_refinement();
3402 *  
3403 *   /*--- Now we prepare the object for transfering, basically saving the old quantities using SolutionTransfer.
3404 *   Since the 'prepare_for_coarsening_and_refinement' method can be called only once, but we have two vectors
3405 *   for dof_handler_velocity, we need to put them in an auxiliary vector. ---*/
3406 *   std::vector<const LinearAlgebra::distributed::Vector<double>*> velocities;
3407 *   velocities.push_back(&u_n);
3408 *   velocities.push_back(&u_n_minus_1);
3410 *   solution_transfer_velocity(dof_handler_velocity);
3411 *   solution_transfer_velocity.prepare_for_coarsening_and_refinement(velocities);
3413 *   solution_transfer_pressure(dof_handler_pressure);
3414 *   solution_transfer_pressure.prepare_for_coarsening_and_refinement(pres_n);
3415 *  
3416 *   triangulation.execute_coarsening_and_refinement(); /*--- Effectively perform the remeshing ---*/
3417 *  
3418 *   /*--- First DoFHandler objects are set up within the new grid ----*/
3419 *   setup_dofs();
3420 *  
3421 *   /*--- Interpolate current solutions to new mesh. This is done using auxliary vectors just for safety,
3422 *   but the new u_n or pres_n could be used. Again, the only point is that the function 'interpolate'
3423 *   can be called once and so the vectors related to 'dof_handler_velocity' have to collected in an auxiliary vector. ---*/
3424 *   LinearAlgebra::distributed::Vector<double> transfer_velocity,
3425 *   transfer_velocity_minus_1,
3426 *   transfer_pressure;
3427 *   transfer_velocity.reinit(u_n);
3428 *   transfer_velocity.zero_out_ghost_values();
3429 *   transfer_velocity_minus_1.reinit(u_n_minus_1);
3430 *   transfer_velocity_minus_1.zero_out_ghost_values();
3431 *   transfer_pressure.reinit(pres_n);
3432 *   transfer_pressure.zero_out_ghost_values();
3433 *  
3434 *   std::vector<LinearAlgebra::distributed::Vector<double>*> transfer_velocities;
3435 *   transfer_velocities.push_back(&transfer_velocity);
3436 *   transfer_velocities.push_back(&transfer_velocity_minus_1);
3437 *   solution_transfer_velocity.interpolate(transfer_velocities);
3438 *   transfer_velocity.update_ghost_values();
3439 *   transfer_velocity_minus_1.update_ghost_values();
3440 *   solution_transfer_pressure.interpolate(transfer_pressure);
3441 *   transfer_pressure.update_ghost_values();
3442 *  
3443 *   u_n = transfer_velocity;
3444 *   u_n_minus_1 = transfer_velocity_minus_1;
3445 *   pres_n = transfer_pressure;
3446 *   }
3447 *  
3448 *  
3449 * @endcode
3450 *
3451 * Interpolate the locally refined solution to a mesh with maximal resolution
3452 * and transfer velocity and pressure.
3453 *
3454
3455 *
3456 *
3457 * @code
3458 *   template<int dim>
3459 *   void NavierStokesProjection<dim>::interpolate_max_res(const unsigned int level) {
3461 *   solution_transfer_velocity(dof_handler_velocity);
3462 *   std::vector<const LinearAlgebra::distributed::Vector<double>*> velocities;
3463 *   velocities.push_back(&u_n);
3464 *   velocities.push_back(&u_n_minus_1);
3465 *   solution_transfer_velocity.prepare_for_coarsening_and_refinement(velocities);
3466 *  
3468 *   solution_transfer_pressure(dof_handler_pressure);
3469 *   solution_transfer_pressure.prepare_for_coarsening_and_refinement(pres_n);
3470 *  
3471 *   for(const auto& cell: triangulation.active_cell_iterators_on_level(level)) {
3472 *   if(cell->is_locally_owned())
3473 *   cell->set_refine_flag();
3474 *   }
3475 *   triangulation.execute_coarsening_and_refinement();
3476 *  
3477 *   setup_dofs();
3478 *  
3479 *   LinearAlgebra::distributed::Vector<double> transfer_velocity, transfer_velocity_minus_1,
3480 *   transfer_pressure;
3481 *  
3482 *   transfer_velocity.reinit(u_n);
3483 *   transfer_velocity.zero_out_ghost_values();
3484 *   transfer_velocity_minus_1.reinit(u_n_minus_1);
3485 *   transfer_velocity_minus_1.zero_out_ghost_values();
3486 *  
3487 *   transfer_pressure.reinit(pres_n);
3488 *   transfer_pressure.zero_out_ghost_values();
3489 *  
3490 *   std::vector<LinearAlgebra::distributed::Vector<double>*> transfer_velocities;
3491 *  
3492 *   transfer_velocities.push_back(&transfer_velocity);
3493 *   transfer_velocities.push_back(&transfer_velocity_minus_1);
3494 *   solution_transfer_velocity.interpolate(transfer_velocities);
3495 *   transfer_velocity.update_ghost_values();
3496 *   transfer_velocity_minus_1.update_ghost_values();
3497 *  
3498 *   solution_transfer_pressure.interpolate(transfer_pressure);
3499 *   transfer_pressure.update_ghost_values();
3500 *  
3501 *   u_n = transfer_velocity;
3502 *   u_n_minus_1 = transfer_velocity_minus_1;
3503 *   pres_n = transfer_pressure;
3504 *   }
3505 *  
3506 *  
3507 * @endcode
3508 *
3509 * Save maximum resolution to a mesh adapted.
3510 *
3511
3512 *
3513 *
3514 * @code
3515 *   template<int dim>
3516 *   void NavierStokesProjection<dim>::save_max_res() {
3517 *   parallel::distributed::Triangulation<dim> triangulation_tmp(MPI_COMM_WORLD);
3518 *   GridGenerator::plate_with_a_hole(triangulation_tmp, 0.5, 1.0, 1.0, 1.1, 1.0, 19.0, Point<2>(2.0, 2.0), 0, 1, 1.0, 2, true);
3519 *   triangulation_tmp.refine_global(triangulation.n_global_levels() - 1);
3520 *  
3521 *   DoFHandler<dim> dof_handler_velocity_tmp(triangulation_tmp);
3522 *   DoFHandler<dim> dof_handler_pressure_tmp(triangulation_tmp);
3523 *   dof_handler_velocity_tmp.distribute_dofs(fe_velocity);
3524 *   dof_handler_pressure_tmp.distribute_dofs(fe_pressure);
3525 *  
3527 *   pres_n_tmp;
3528 *   u_n_tmp.reinit(dof_handler_velocity_tmp.n_dofs());
3529 *   pres_n_tmp.reinit(dof_handler_pressure_tmp.n_dofs());
3530 *  
3531 *   DataOut<dim> data_out;
3532 *   std::vector<std::string> velocity_names(dim, "v");
3533 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
3534 *   component_interpretation_velocity(dim, DataComponentInterpretation::component_is_part_of_vector);
3535 *   VectorTools::interpolate_to_different_mesh(dof_handler_velocity, u_n, dof_handler_velocity_tmp, u_n_tmp);
3536 *   u_n_tmp.update_ghost_values();
3537 *   data_out.add_data_vector(dof_handler_velocity_tmp, u_n_tmp, velocity_names, component_interpretation_velocity);
3538 *   VectorTools::interpolate_to_different_mesh(dof_handler_pressure, pres_n, dof_handler_pressure_tmp, pres_n_tmp);
3539 *   pres_n_tmp.update_ghost_values();
3540 *   data_out.add_data_vector(dof_handler_pressure_tmp, pres_n_tmp, "p", {DataComponentInterpretation::component_is_scalar});
3541 *   PostprocessorVorticity<dim> postprocessor;
3542 *   data_out.add_data_vector(dof_handler_velocity_tmp, u_n_tmp, postprocessor);
3543 *  
3544 *   data_out.build_patches(MappingQ1<dim>(), 1, DataOut<dim>::curved_inner_cells);
3545 *   const std::string output = "./" + saving_dir + "/solution_max_res_end.vtu";
3546 *   data_out.write_vtu_in_parallel(output, MPI_COMM_WORLD);
3547 *   }
3548 *  
3549 *  
3550 * @endcode
3551 *
3552 *
3553 * <a name=""></a>
3554 * @sect{ <code>NavierStokesProjection::run</code> }
3555 *
3556
3557 *
3558 * This is the time marching function, which starting at <code>t_0</code>
3559 * advances in time using the projection method with time step <code>dt</code>
3560 * until <code>T</code>.
3561 *
3562
3563 *
3564 * Its second parameter, <code>verbose</code> indicates whether the function
3565 * should output information what it is doing at any given moment:
3566 * we use the ConditionalOStream class to do that for us.
3567 *
3568
3569 *
3570 *
3571 * @code
3572 *   template<int dim>
3573 *   void NavierStokesProjection<dim>::run(const bool verbose, const unsigned int output_interval) {
3574 *   ConditionalOStream verbose_cout(std::cout, verbose && Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0);
3575 *  
3576 *   output_results(1);
3577 *   double time = t_0 + dt;
3578 *   unsigned int n = 1;
3579 *   while(std::abs(T - time) > 1e-10) {
3580 *   time += dt;
3581 *   n++;
3582 *   pcout << "Step = " << n << " Time = " << time << std::endl;
3583 *  
3584 *   /*--- First stage of TR-BDF2 and we start by setting the proper flag ---*/
3585 *   TR_BDF2_stage = 1;
3586 *   navier_stokes_matrix.set_TR_BDF2_stage(TR_BDF2_stage);
3587 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
3588 *   mg_matrices[level].set_TR_BDF2_stage(TR_BDF2_stage);
3589 *  
3590 *   verbose_cout << " Interpolating the velocity stage 1" << std::endl;
3591 *   interpolate_velocity();
3592 *  
3593 *   verbose_cout << " Diffusion Step stage 1 " << std::endl;
3594 *   diffusion_step();
3595 *  
3596 *   verbose_cout << " Projection Step stage 1" << std::endl;
3597 *   project_grad(1);
3598 *   u_tmp.equ(gamma*dt, u_tmp);
3599 *   u_star += u_tmp; /*--- In the rhs of the projection step we need u_star + gamma*dt*grad(pres_n) and we save it into u_star ---*/
3600 *   projection_step();
3601 *  
3602 *   verbose_cout << " Updating the Velocity stage 1" << std::endl;
3603 *   u_n_gamma.equ(1.0, u_star);
3604 *   project_grad(2);
3605 *   grad_pres_int.equ(1.0, u_tmp); /*--- We save grad(pres_int), because we will need it soon ---*/
3606 *   u_tmp.equ(-gamma*dt, u_tmp);
3607 *   u_n_gamma += u_tmp; /*--- u_n_gamma = u_star - gamma*dt*grad(pres_int) ---*/
3608 *   u_n_minus_1 = u_n;
3609 *  
3610 *   /*--- Second stage of TR-BDF2 ---*/
3611 *   TR_BDF2_stage = 2;
3612 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
3613 *   mg_matrices[level].set_TR_BDF2_stage(TR_BDF2_stage);
3614 *   navier_stokes_matrix.set_TR_BDF2_stage(TR_BDF2_stage);
3615 *  
3616 *   verbose_cout << " Interpolating the velocity stage 2" << std::endl;
3617 *   interpolate_velocity();
3618 *  
3619 *   verbose_cout << " Diffusion Step stage 2 " << std::endl;
3620 *   diffusion_step();
3621 *  
3622 *   verbose_cout << " Projection Step stage 2" << std::endl;
3623 *   u_tmp.equ((1.0 - gamma)*dt, grad_pres_int);
3624 *   u_star += u_tmp; /*--- In the rhs of the projection step we need u_star + (1 - gamma)*dt*grad(pres_int) ---*/
3625 *   projection_step();
3626 *  
3627 *   verbose_cout << " Updating the Velocity stage 2" << std::endl;
3628 *   u_n.equ(1.0, u_star);
3629 *   project_grad(1);
3630 *   u_tmp.equ((gamma - 1.0)*dt, u_tmp);
3631 *   u_n += u_tmp; /*--- u_n = u_star - (1 - gamma)*dt*grad(pres_n) ---*/
3632 *  
3633 *   const double max_vel = get_maximal_velocity();
3634 *   pcout<< "Maximal velocity = " << max_vel << std::endl;
3635 *   /*--- The Courant number is computed taking into account the polynomial degree for the velocity ---*/
3636 *   pcout << "CFL = " << dt*max_vel*(EquationData::degree_p + 1)*
3638 *   compute_lift_and_drag();
3639 *   if(n % output_interval == 0) {
3640 *   verbose_cout << "Plotting Solution final" << std::endl;
3641 *   output_results(n);
3642 *   }
3643 *   /*--- In case dt is not a multiple of T, we reduce dt in order to end up at T ---*/
3644 *   if(T - time < dt && T - time > 1e-10) {
3645 *   dt = T - time;
3646 *   navier_stokes_matrix.set_dt(dt);
3647 *   for(unsigned int level = 0; level < triangulation.n_global_levels(); ++level)
3648 *   mg_matrices[level].set_dt(dt);
3649 *   }
3650 *   /*--- Perform the refinement if desired ---*/
3651 *   if(refinement_iterations > 0 && n % refinement_iterations == 0) {
3652 *   verbose_cout << "Refining mesh" << std::endl;
3653 *   refine_mesh();
3654 *   }
3655 *   }
3656 *   if(n % output_interval != 0) {
3657 *   verbose_cout << "Plotting Solution final" << std::endl;
3658 *   output_results(n);
3659 *   }
3660 *   if(refinement_iterations > 0) {
3661 *   for(unsigned int lev = 0; lev < triangulation.n_global_levels() - 1; ++ lev)
3662 *   interpolate_max_res(lev);
3663 *   save_max_res();
3664 *   }
3665 *   }
3666 *  
3667 *   } // namespace NS_TRBDF2
3668 *  
3669 *  
3670 * @endcode
3671 *
3672 *
3673 * <a name=""></a>
3674 * @sect{ The main function }
3675 *
3676
3677 *
3678 * The main function looks very much like in all the other tutorial programs. We first initialize MPI,
3679 * we initialize the class 'NavierStokesProjection' with the dimension as template parameter and then
3680 * let the method 'run' do the job.
3681 *
3682
3683 *
3684 *
3685 * @code
3686 *   int main(int argc, char *argv[]) {
3687 *   try {
3688 *   using namespace NS_TRBDF2;
3689 *  
3690 *   RunTimeParameters::Data_Storage data;
3691 *   data.read_data("parameter-file.prm");
3692 *  
3693 *   Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, -1);
3694 *  
3695 *   const auto& curr_rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
3696 *   deallog.depth_console(data.verbose && curr_rank == 0 ? 2 : 0);
3697 *  
3698 *   NavierStokesProjection<2> test(data);
3699 *   test.run(data.verbose, data.output_interval);
3700 *  
3701 *   if(curr_rank == 0)
3702 *   std::cout << "----------------------------------------------------"
3703 *   << std::endl
3704 *   << "Apparently everything went fine!" << std::endl
3705 *   << "Don't forget to brush your teeth :-)" << std::endl
3706 *   << std::endl;
3707 *  
3708 *   return 0;
3709 *   }
3710 *   catch(std::exception &exc) {
3711 *   std::cerr << std::endl
3712 *   << std::endl
3713 *   << "----------------------------------------------------"
3714 *   << std::endl;
3715 *   std::cerr << "Exception on processing: " << std::endl
3716 *   << exc.what() << std::endl
3717 *   << "Aborting!" << std::endl
3718 *   << "----------------------------------------------------"
3719 *   << std::endl;
3720 *   return 1;
3721 *   }
3722 *   catch(...) {
3723 *   std::cerr << std::endl
3724 *   << std::endl
3725 *   << "----------------------------------------------------"
3726 *   << std::endl;
3727 *   std::cerr << "Unknown exception!" << std::endl
3728 *   << "Aborting!" << std::endl
3729 *   << "----------------------------------------------------"
3730 *   << std::endl;
3731 *   return 1;
3732 *   }
3733 *  
3734 *   }
3735 * @endcode
3736
3737
3738<a name="ann-runtime_parameters.h"></a>
3739<h1>Annotated version of runtime_parameters.h</h1>
3740 *
3741 *
3742 *
3743 * We start by including all the necessary deal.II header files
3744 *
3745
3746 *
3747 *
3748 * @code
3749 *   #include <deal.II/base/parameter_handler.h>
3750 *  
3751 * @endcode
3752 *
3753 *
3754 * <a name=""></a>
3755 * @sect{Run time parameters}
3756 *
3757
3758 *
3759 * Since our method has several parameters that can be fine-tuned we put them
3760 * into an external file, so that they can be determined at run-time.
3761 *
3762
3763 *
3764 *
3765 * @code
3766 *   namespace RunTimeParameters {
3767 *   using namespace dealii;
3768 *  
3769 *   class Data_Storage {
3770 *   public:
3771 *   Data_Storage();
3772 *  
3773 *   void read_data(const std::string& filename);
3774 *  
3775 *   double initial_time;
3776 *   double final_time;
3777 *  
3778 *   double Reynolds;
3779 *   double dt;
3780 *  
3781 *   unsigned int n_refines; /*--- Number of refinements ---*/
3782 *   unsigned int max_loc_refinements; /*--- Number of maximum local refinements allowed ---*/
3783 *   unsigned int min_loc_refinements; /*--- Number of minimum local refinements allowed
3784 *   once reached that level ---*/
3785 *  
3786 *   /*--- Parameters related to the linear solver ---*/
3787 *   unsigned int max_iterations;
3788 *   double eps;
3789 *  
3790 *   bool verbose;
3791 *   unsigned int output_interval;
3792 *  
3793 *   std::string dir; /*--- Auxiliary string variable for output storage ---*/
3794 *  
3795 *   unsigned int refinement_iterations; /*--- Auxiliary variable about how many steps perform remeshing ---*/
3796 *  
3797 *   protected:
3798 *   ParameterHandler prm;
3799 *   };
3800 *  
3801 * @endcode
3802 *
3803 * In the constructor of this class we declare all the parameters in suitable (but arbitrary) subsections.
3804 *
3805
3806 *
3807 *
3808 * @code
3809 *   Data_Storage::Data_Storage(): initial_time(0.0),
3810 *   final_time(1.0),
3811 *   Reynolds(1.0),
3812 *   dt(5e-4),
3813 *   n_refines(0),
3814 *   max_loc_refinements(0),
3815 *   min_loc_refinements(0),
3816 *   max_iterations(1000),
3817 *   eps(1e-12),
3818 *   verbose(true),
3819 *   output_interval(15),
3820 *   refinement_iterations(0) {
3821 *   prm.enter_subsection("Physical data");
3822 *   {
3823 *   prm.declare_entry("initial_time",
3824 *   "0.0",
3825 *   Patterns::Double(0.0),
3826 *   " The initial time of the simulation. ");
3827 *   prm.declare_entry("final_time",
3828 *   "1.0",
3829 *   Patterns::Double(0.0),
3830 *   " The final time of the simulation. ");
3831 *   prm.declare_entry("Reynolds",
3832 *   "1.0",
3833 *   Patterns::Double(0.0),
3834 *   " The Reynolds number. ");
3835 *   }
3836 *   prm.leave_subsection();
3837 *  
3838 *   prm.enter_subsection("Time step data");
3839 *   {
3840 *   prm.declare_entry("dt",
3841 *   "5e-4",
3842 *   Patterns::Double(0.0),
3843 *   " The time step size. ");
3844 *   }
3845 *   prm.leave_subsection();
3846 *  
3847 *   prm.enter_subsection("Space discretization");
3848 *   {
3849 *   prm.declare_entry("n_of_refines",
3850 *   "100",
3851 *   Patterns::Integer(0, 1500),
3852 *   " The number of cells we want on each direction of the mesh. ");
3853 *   prm.declare_entry("max_loc_refinements",
3854 *   "4",
3855 *   Patterns::Integer(0, 10),
3856 *   " The number of maximum local refinements. ");
3857 *   prm.declare_entry("min_loc_refinements",
3858 *   "2",
3859 *   Patterns::Integer(0, 10),
3860 *   " The number of minimum local refinements. ");
3861 *   }
3862 *   prm.leave_subsection();
3863 *  
3864 *   prm.enter_subsection("Data solve");
3865 *   {
3866 *   prm.declare_entry("max_iterations",
3867 *   "1000",
3868 *   Patterns::Integer(1, 30000),
3869 *   " The maximal number of iterations linear solvers must make. ");
3870 *   prm.declare_entry("eps",
3871 *   "1e-12",
3872 *   Patterns::Double(0.0),
3873 *   " The stopping criterion. ");
3874 *   }
3875 *   prm.leave_subsection();
3876 *  
3877 *   prm.declare_entry("refinement_iterations",
3878 *   "0",
3879 *   Patterns::Integer(0),
3880 *   " This number indicates how often we need to "
3881 *   "refine the mesh");
3882 *  
3883 *   prm.declare_entry("saving directory", "SimTest");
3884 *  
3885 *   prm.declare_entry("verbose",
3886 *   "true",
3887 *   Patterns::Bool(),
3888 *   " This indicates whether the output of the solution "
3889 *   "process should be verbose. ");
3890 *  
3891 *   prm.declare_entry("output_interval",
3892 *   "1",
3893 *   Patterns::Integer(1),
3894 *   " This indicates between how many time steps we print "
3895 *   "the solution. ");
3896 *   }
3897 *  
3898 * @endcode
3899 *
3900 * We need now a routine to read all declared parameters in the constructor
3901 *
3902
3903 *
3904 *
3905 * @code
3906 *   void Data_Storage::read_data(const std::string& filename) {
3907 *   std::ifstream file(filename);
3908 *   AssertThrow(file, ExcFileNotOpen(filename));
3909 *  
3910 *   prm.parse_input(file);
3911 *  
3912 *   prm.enter_subsection("Physical data");
3913 *   {
3914 *   initial_time = prm.get_double("initial_time");
3915 *   final_time = prm.get_double("final_time");
3916 *   Reynolds = prm.get_double("Reynolds");
3917 *   }
3918 *   prm.leave_subsection();
3919 *  
3920 *   prm.enter_subsection("Time step data");
3921 *   {
3922 *   dt = prm.get_double("dt");
3923 *   }
3924 *   prm.leave_subsection();
3925 *  
3926 *   prm.enter_subsection("Space discretization");
3927 *   {
3928 *   n_refines = prm.get_integer("n_of_refines");
3929 *   max_loc_refinements = prm.get_integer("max_loc_refinements");
3930 *   min_loc_refinements = prm.get_integer("min_loc_refinements");
3931 *   }
3932 *   prm.leave_subsection();
3933 *  
3934 *   prm.enter_subsection("Data solve");
3935 *   {
3936 *   max_iterations = prm.get_integer("max_iterations");
3937 *   eps = prm.get_double("eps");
3938 *   }
3939 *   prm.leave_subsection();
3940 *  
3941 *   dir = prm.get("saving directory");
3942 *  
3943 *   refinement_iterations = prm.get_integer("refinement_iterations");
3944 *  
3945 *   verbose = prm.get_bool("verbose");
3946 *  
3947 *   output_interval = prm.get_integer("output_interval");
3948 *   }
3949 *  
3950 *   } // namespace RunTimeParameters
3951 * @endcode
3952
3953
3954*/
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
virtual void vector_value(const Point< dim > &p, Vector< RangeNumberType > &values) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:350
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&...args)
void build(const DoFHandler< dim, dim > &dof_handler, const std::vector< std::shared_ptr< const Utilities::MPI::Partitioner > > &external_partitioners=std::vector< std::shared_ptr< const Utilities::MPI::Partitioner > >())
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
void clear()
unsigned int mg_level
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
Definition point.h:112
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
unsigned int cell_index
__global__ void set(Number *val, const Number s, const size_type N)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
typename ActiveSelector::active_cell_iterator active_cell_iterator
void loop(ITERATOR begin, std_cxx20::type_identity_t< ITERATOR > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:439
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:282
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.
LogStream deallog
Definition logstream.cc:37
const Event initial
Definition event.cc:65
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > &center=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
double minimal_cell_diameter(const Triangulation< dim, spacedim > &triangulation, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
void compute_diagonal(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, VectorType &diagonal_global, const std::function< void(FEEvaluation< dim, fe_degree, n_q_points_1d, n_components, Number, VectorizedArrayType > &)> &local_vmult, const unsigned int dof_no=0, const unsigned int quad_no=0, const unsigned int first_selected_component=0)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
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)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:150
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:161
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:471
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask=ComponentMask())
void interpolate_to_different_mesh(const DoFHandler< dim, spacedim > &dof1, const VectorType &u1, const DoFHandler< dim, spacedim > &dof2, VectorType &u2)
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 >(0)), 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)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition loop.h:71
Definition mg.h:82
static const unsigned int invalid_unsigned_int
Definition types.h:213
void refine_and_coarsen_fixed_number(::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
STL namespace.
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int boundary_id
Definition types.h:141
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector<::Vector< double > > solution_values
std::vector< std::vector< Tensor< 1, spacedim > > > solution_gradients
TasksParallelScheme tasks_parallel_scheme
DEAL_II_HOST constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
std_cxx20::type_identity< T > identity
const ::Triangulation< dim, spacedim > & tria