Loading [MathJax]/extensions/TeX/AMSsymbols.js
 Reference documentation for deal.II version 9.4.1
\(\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\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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>
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);
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);
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);
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,
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;
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);
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,
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") ---*/
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 *
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 *
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 *
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. ---*/
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 *
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 add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1064
const unsigned int n_quadrature_points
Definition: fe_values.h:2432
void get_function_gradients(const InputVector &fe_function, std::vector< Tensor< 1, spacedim, typename InputVector::value_type > > &gradients) const
Definition: fe_values.cc:3495
double JxW(const unsigned int quadrature_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
Definition: fe_dgq.h:111
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
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()
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:111
Definition: tensor.h:503
Definition: vector.h:109
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
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.
Point< 2 > second
Definition: grid_out.cc:4604
Point< 2 > first
Definition: grid_out.cc:4603
unsigned int level
Definition: grid_out.cc:4606
unsigned int cell_index
Definition: grid_tools.cc:1129
__global__ void set(Number *val, const Number s, const size_type N)
void write_vtu_in_parallel(const std::string &filename, const MPI_Comm &comm) const
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition: dof_handler.h:438
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
void loop(ITERATOR begin, typename identity< ITERATOR >::type 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 reinit(const size_type size, const bool omit_zeroing_entries=false)
LogStream deallog
Definition: logstream.cc:37
const Event initial
Definition: event.cc:65
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
Definition: dof_tools.cc:1144
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 >()))
Definition: grid_tools.cc:4390
@ 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:190
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)
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:151
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:140
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:473
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 typename identity< Iterator >::type &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)
Definition: work_stream.h:474
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:201
void refine_and_coarsen_fixed_number(parallel::distributed::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:129
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
Definition: matrix_free.h:343
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:407
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:387
UpdateFlags mapping_update_flags
Definition: matrix_free.h:367
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
const ::Triangulation< dim, spacedim > & tria