deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-55.h
Go to the documentation of this file.
1
234 *  
235 *   namespace LA
236 *   {
239 *   using namespace dealii::LinearAlgebraPETSc;
242 *   using namespace dealii::LinearAlgebraTrilinos;
243 *   #else
245 *   #endif
246 *   } // namespace LA
247 *  
248 *   #include <deal.II/lac/vector.h>
249 *   #include <deal.II/lac/full_matrix.h>
250 *   #include <deal.II/lac/solver_cg.h>
251 *   #include <deal.II/lac/solver_gmres.h>
253 *   #include <deal.II/lac/affine_constraints.h>
255 *  
257 *   #include <deal.II/lac/petsc_vector.h>
258 *   #include <deal.II/lac/petsc_solver.h>
260 *  
261 *   #include <deal.II/grid/grid_generator.h>
262 *   #include <deal.II/grid/manifold_lib.h>
263 *   #include <deal.II/grid/grid_tools.h>
264 *   #include <deal.II/dofs/dof_handler.h>
265 *   #include <deal.II/dofs/dof_renumbering.h>
266 *   #include <deal.II/dofs/dof_tools.h>
267 *   #include <deal.II/fe/fe_values.h>
268 *   #include <deal.II/fe/fe_q.h>
269 *   #include <deal.II/fe/fe_system.h>
270 *   #include <deal.II/numerics/vector_tools.h>
271 *   #include <deal.II/numerics/data_out.h>
272 *   #include <deal.II/numerics/error_estimator.h>
273 *  
274 *   #include <deal.II/base/utilities.h>
275 *   #include <deal.II/base/conditional_ostream.h>
276 *   #include <deal.II/base/index_set.h>
278 *   #include <deal.II/distributed/tria.h>
279 *   #include <deal.II/distributed/grid_refinement.h>
280 *  
284 *  
285 *   namespace Step55
286 *   {
287 *   using namespace dealii;
288 *  
289 * @endcode
290 *
291 *
292 * <a name="step_55-Linearsolversandpreconditioners"></a>
294 *
295
296 *
297 * We need a few helper classes to represent our solver strategy
299 *
300
301 *
302 *
303 * @code
304 *   namespace LinearSolvers
305 *   {
306 * @endcode
307 *
308 * This class exposes the action of applying the inverse of a
309 * giving matrix via the function
310 * InverseMatrix::vmult(). Internally, the inverse is not formed
311 * explicitly. Instead, a linear solver with CG is performed. This
312 * class extends the InverseMatrix class in @ref step_22 "step-22" with an option
313 * to specify a preconditioner, and to allow for different vector
314 * types in the vmult function. We use the same mechanism as in
315 * @ref step_31 "step-31" to convert a run-time exception into a failed assertion
316 * should the inner solver not converge.
317 *
318 * @code
319 *   template <class Matrix, class Preconditioner>
320 *   class InverseMatrix : public EnableObserverPointer
321 *   {
322 *   public:
323 *   InverseMatrix(const Matrix &m, const Preconditioner &preconditioner);
324 *  
325 *   template <typename VectorType>
326 *   void vmult(VectorType &dst, const VectorType &src) const;
327 *  
328 *   private:
330 *   const Preconditioner &preconditioner;
331 *   };
332 *  
333 *  
334 *   template <class Matrix, class Preconditioner>
336 *   const Matrix &m,
337 *   const Preconditioner &preconditioner)
338 *   : matrix(&m)
339 *   , preconditioner(preconditioner)
340 *   {}
341 *  
342 *  
343 *  
344 *   template <class Matrix, class Preconditioner>
345 *   template <typename VectorType>
346 *   void
348 *   const VectorType &src) const
349 *   {
350 *   SolverControl solver_control(src.size(), 1e-8 * src.l2_norm());
351 *   SolverCG<VectorType> cg(solver_control);
352 *   dst = 0;
353 *  
354 *   try
355 *   {
356 *   cg.solve(*matrix, dst, src, preconditioner);
357 *   }
358 *   catch (std::exception &e)
359 *   {
360 *   Assert(false, ExcMessage(e.what()));
361 *   }
362 *   }
363 *  
364 *  
365 * @endcode
366 *
367 * The class A template class for a simple block diagonal preconditioner
368 * for 2x2 matrices.
369 *
370 * @code
371 *   template <class PreconditionerA, class PreconditionerS>
373 *   {
374 *   public:
377 *  
378 *   void vmult(LA::MPI::BlockVector &dst,
379 *   const LA::MPI::BlockVector &src) const;
380 *  
381 *   private:
384 *   };
385 *  
386 *   template <class PreconditionerA, class PreconditionerS>
392 *   {}
393 *  
394 *  
395 *   template <class PreconditionerA, class PreconditionerS>
397 *   LA::MPI::BlockVector &dst,
398 *   const LA::MPI::BlockVector &src) const
399 *   {
400 *   preconditioner_A.vmult(dst.block(0), src.block(0));
401 *   preconditioner_S.vmult(dst.block(1), src.block(1));
402 *   }
403 *  
404 *   } // namespace LinearSolvers
405 *  
406 * @endcode
407 *
408 *
409 * <a name="step_55-Problemsetup"></a>
410 * <h3>Problem setup</h3>
411 *
412
413 *
415 * solution for the test problem.
416 *
417
418 *
419 *
420 * @code
421 *   template <int dim>
422 *   class RightHandSide : public Function<dim>
423 *   {
424 *   public:
425 *   RightHandSide()
426 *   : Function<dim>(dim + 1)
427 *   {}
428 *  
429 *   virtual void vector_value(const Point<dim> &p,
430 *   Vector<double> &value) const override;
431 *   };
432 *  
433 *  
434 *   template <int dim>
436 *   Vector<double> &values) const
437 *   {
438 *   const double R_x = p[0];
439 *   const double R_y = p[1];
440 *  
441 *   constexpr double pi = numbers::PI;
442 *   constexpr double pi2 = numbers::PI * numbers::PI;
443 *  
444 * @endcode
445 *
446 * velocity
447 *
448 * @code
449 *   values[0] = -1.0L / 2.0L * (-2 * std::sqrt(25.0 + 4 * pi2) + 10.0) *
450 *   std::exp(R_x * (-2 * std::sqrt(25.0 + 4 * pi2) + 10.0)) -
451 *   0.4 * pi2 * std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
452 *   std::cos(2 * R_y * pi) +
453 *   0.1 *
454 *   Utilities::fixed_power<2>(-std::sqrt(25.0 + 4 * pi2) + 5.0) *
455 *   std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
456 *   std::cos(2 * R_y * pi);
457 *   values[1] = 0.2 * pi * (-std::sqrt(25.0 + 4 * pi2) + 5.0) *
458 *   std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
459 *   std::sin(2 * R_y * pi) -
460 *   0.05 *
461 *   Utilities::fixed_power<3>(-std::sqrt(25.0 + 4 * pi2) + 5.0) *
462 *   std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
463 *   std::sin(2 * R_y * pi) / pi;
464 *  
465 * @endcode
466 *
467 * pressure
468 *
469 * @code
470 *   values[dim] = 0;
471 *   }
472 *  
473 *  
474 *   template <int dim>
475 *   class ExactSolution : public Function<dim>
476 *   {
477 *   public:
478 *   ExactSolution()
479 *   : Function<dim>(dim + 1)
480 *   {}
481 *  
482 *   virtual void vector_value(const Point<dim> &p,
483 *   Vector<double> &values) const override;
484 *   };
485 *  
486 *   template <int dim>
488 *   Vector<double> &values) const
489 *   {
490 *   const double R_x = p[0];
491 *   const double R_y = p[1];
492 *  
493 *   constexpr double pi = numbers::PI;
494 *   constexpr double pi2 = numbers::PI * numbers::PI;
495 *  
496 * @endcode
497 *
498 * velocity
499 *
500 * @code
501 *   values[0] = -std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
502 *   std::cos(2 * R_y * pi) +
503 *   1;
504 *   values[1] = (1.0L / 2.0L) * (-std::sqrt(25.0 + 4 * pi2) + 5.0) *
505 *   std::exp(R_x * (-std::sqrt(25.0 + 4 * pi2) + 5.0)) *
506 *   std::sin(2 * R_y * pi) / pi;
507 *  
508 * @endcode
509 *
510 * pressure
511 *
512 * @code
513 *   values[dim] =
514 *   -1.0L / 2.0L * std::exp(R_x * (-2 * std::sqrt(25.0 + 4 * pi2) + 10.0)) -
515 *   2.0 *
516 *   (-6538034.74494422 +
517 *   0.0134758939981709 * std::exp(4 * std::sqrt(25.0 + 4 * pi2))) /
518 *   (-80.0 * std::exp(3 * std::sqrt(25.0 + 4 * pi2)) +
519 *   16.0 * std::sqrt(25.0 + 4 * pi2) *
520 *   std::exp(3 * std::sqrt(25.0 + 4 * pi2))) -
521 *   1634508.68623606 * std::exp(-3.0 * std::sqrt(25.0 + 4 * pi2)) /
522 *   (-10.0 + 2.0 * std::sqrt(25.0 + 4 * pi2)) +
523 *   (-0.00673794699908547 * std::exp(std::sqrt(25.0 + 4 * pi2)) +
524 *   3269017.37247211 * std::exp(-3 * std::sqrt(25.0 + 4 * pi2))) /
525 *   (-8 * std::sqrt(25.0 + 4 * pi2) + 40.0) +
526 *   0.00336897349954273 * std::exp(1.0 * std::sqrt(25.0 + 4 * pi2)) /
527 *   (-10.0 + 2.0 * std::sqrt(25.0 + 4 * pi2));
528 *   }
529 *  
530 *  
531 *  
532 * @endcode
533 *
534 *
535 * <a name="step_55-Themainprogram"></a>
536 * <h3>The main program</h3>
537 *
538
539 *
540 * The main class is very similar to @ref step_40 "step-40", except that matrices and
541 * vectors are now block versions, and we store a std::vector<IndexSet>
543 * exactly two IndexSets, one for all velocity unknowns and one for all
545 *
546 * @code
547 *   template <int dim>
548 *   class StokesProblem
549 *   {
550 *   public:
551 *   StokesProblem(unsigned int velocity_degree);
552 *  
553 *   void run();
554 *  
555 *   private:
556 *   void make_grid();
557 *   void setup_system();
558 *   void assemble_system();
559 *   void solve();
560 *   void refine_grid();
561 *   void output_results(const unsigned int cycle);
562 *  
563 *   const unsigned int velocity_degree;
564 *   const double viscosity;
565 *   MPI_Comm mpi_communicator;
566 *  
567 *   const FESystem<dim> fe;
569 *   DoFHandler<dim> dof_handler;
570 *  
571 *   std::vector<IndexSet> owned_partitioning;
572 *   std::vector<IndexSet> relevant_partitioning;
573 *  
574 *   AffineConstraints<double> constraints;
575 *  
576 *   LA::MPI::BlockSparseMatrix system_matrix;
577 *   LA::MPI::BlockSparseMatrix preconditioner_matrix;
578 *   LA::MPI::BlockVector locally_relevant_solution;
579 *   LA::MPI::BlockVector system_rhs;
580 *  
583 *   };
584 *  
585 *  
586 *  
587 *   template <int dim>
590 *   , viscosity(0.1)
591 *   , mpi_communicator(MPI_COMM_WORLD)
593 *   , triangulation(mpi_communicator,
597 *   , dof_handler(triangulation)
598 *   , pcout(std::cout,
599 *   (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
600 *   , computing_timer(mpi_communicator,
601 *   pcout,
604 *   {}
605 *  
606 *  
607 * @endcode
608 *
609 * The Kovasznay flow is defined on the domain [-0.5, 1.5]^2, which we
611 *
612 * @code
613 *   template <int dim>
615 *   {
617 *   triangulation.refine_global(3);
618 *   }
619 *  
620 * @endcode
621 *
622 *
623 * <a name="step_55-SystemSetup"></a>
624 * <h3>System Setup</h3>
625 *
626
627 *
628 * The construction of the block matrices and vectors is new compared to
629 * @ref step_40 "step-40" and is different compared to serial codes like @ref step_22 "step-22", because
631 *
632 * @code
633 *   template <int dim>
635 *   {
637 *  
638 *   dof_handler.distribute_dofs(fe);
639 *  
640 * @endcode
641 *
642 * Put all dim velocities into block 0 and the pressure into block 1,
644 * we have per block.
645 *
646 * @code
647 *   std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
648 *   stokes_sub_blocks[dim] = 1;
650 *  
651 *   const std::vector<types::global_dof_index> dofs_per_block =
653 *  
654 *   const unsigned int n_u = dofs_per_block[0];
655 *   const unsigned int n_p = dofs_per_block[1];
656 *  
657 *   pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() << " ("
658 *   << n_u << '+' << n_p << ')' << std::endl;
659 *  
660 * @endcode
661 *
663 * into two IndexSets based on how we want to create the block matrices
664 * and vectors.
665 *
666 * @code
667 *   const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs();
668 *   owned_partitioning = {locally_owned_dofs.get_view(0, n_u),
669 *   locally_owned_dofs.get_view(n_u, n_u + n_p)};
670 *  
674 *   locally_relevant_dofs.get_view(n_u, n_u + n_p)};
675 *  
676 * @endcode
677 *
678 * Setting up the constraints for boundary conditions and hanging nodes
679 * is identical to @ref step_40 "step-40". Even though we don't have any hanging nodes
680 * because we only perform global refinement, it is still a good idea
681 * to put this function call in, in case adaptive refinement gets
682 * introduced later.
683 *
684 * @code
685 *   {
686 *   constraints.reinit(locally_owned_dofs, locally_relevant_dofs);
687 *  
688 *   const FEValuesExtractors::Vector velocities(0);
689 *   DoFTools::make_hanging_node_constraints(dof_handler, constraints);
690 *   VectorTools::interpolate_boundary_values(dof_handler,
691 *   0,
692 *   ExactSolution<dim>(),
693 *   constraints,
694 *   fe.component_mask(velocities));
695 *   constraints.close();
696 *   }
697 *  
698 * @endcode
699 *
700 * Now we create the system matrix based on a BlockDynamicSparsityPattern.
701 * We know that we won't have coupling between different velocity
702 * components (because we use the laplace and not the deformation tensor)
706 *
707 * @code
708 *   {
709 *   system_matrix.clear();
710 *  
711 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
712 *   for (unsigned int c = 0; c < dim + 1; ++c)
713 *   for (unsigned int d = 0; d < dim + 1; ++d)
714 *   if (c == dim && d == dim)
715 *   coupling[c][d] = DoFTools::none;
716 *   else if (c == dim || d == dim || c == d)
718 *   else
719 *   coupling[c][d] = DoFTools::none;
720 *  
722 *  
724 *   dof_handler, coupling, dsp, constraints, false);
725 *  
727 *   dsp,
728 *   dof_handler.locally_owned_dofs(),
729 *   mpi_communicator,
731 *  
732 *   system_matrix.reinit(owned_partitioning, dsp, mpi_communicator);
733 *   }
734 *  
735 * @endcode
736 *
737 * The preconditioner matrix has a different coupling (we only fill in
738 * the 1,1 block with the @ref GlossMassMatrix "mass matrix"), otherwise this code is identical
739 * to the construction of the system_matrix above.
740 *
741 * @code
742 *   {
744 *  
745 *   Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
746 *   for (unsigned int c = 0; c < dim + 1; ++c)
747 *   for (unsigned int d = 0; d < dim + 1; ++d)
748 *   if (c == dim && d == dim)
750 *   else
751 *   coupling[c][d] = DoFTools::none;
752 *  
754 *  
756 *   dof_handler, coupling, dsp, constraints, false);
758 *   dsp,
759 *   Utilities::MPI::all_gather(mpi_communicator,
760 *   dof_handler.locally_owned_dofs()),
761 *   mpi_communicator,
763 *   preconditioner_matrix.reinit(owned_partitioning, dsp, mpi_communicator);
764 *   }
765 *  
766 * @endcode
767 *
768 * Finally, we construct the block vectors with the right sizes. The
769 * function call with two std::vector<IndexSet> will create a ghosted
770 * vector.
771 *
772 * @code
775 *   mpi_communicator);
776 *   system_rhs.reinit(owned_partitioning, mpi_communicator);
777 *   }
778 *  
779 *  
780 *  
781 * @endcode
782 *
783 *
784 * <a name="step_55-Assembly"></a>
785 * <h3>Assembly</h3>
786 *
787
788 *
789 * This function assembles the system matrix, the preconditioner matrix,
790 * and the right hand side. The code is pretty standard.
791 *
792 * @code
793 *   template <int dim>
795 *   {
796 *   TimerOutput::Scope t(computing_timer, "assembly");
797 *  
798 *   system_matrix = 0;
800 *   system_rhs = 0;
801 *  
803 *  
804 *   FEValues<dim> fe_values(fe,
808 *  
809 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
810 *   const unsigned int n_q_points = quadrature_formula.size();
811 *  
812 *   FullMatrix<double> system_cell_matrix(dofs_per_cell, dofs_per_cell);
813 *   FullMatrix<double> preconditioner_cell_matrix(dofs_per_cell, dofs_per_cell);
814 *   Vector<double> cell_rhs(dofs_per_cell);
815 *  
817 *   std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
818 *  
819 *   std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
820 *   std::vector<double> div_phi_u(dofs_per_cell);
821 *   std::vector<double> phi_p(dofs_per_cell);
822 *  
823 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
826 *  
827 *   for (const auto &cell : dof_handler.active_cell_iterators())
828 *   if (cell->is_locally_owned())
829 *   {
832 *   cell_rhs = 0;
833 *  
834 *   fe_values.reinit(cell);
835 *   right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
836 *   rhs_values);
837 *   for (unsigned int q = 0; q < n_q_points; ++q)
838 *   {
839 *   for (unsigned int k = 0; k < dofs_per_cell; ++k)
840 *   {
841 *   grad_phi_u[k] = fe_values[velocities].gradient(k, q);
842 *   div_phi_u[k] = fe_values[velocities].divergence(k, q);
843 *   phi_p[k] = fe_values[pressure].value(k, q);
844 *   }
845 *  
846 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
847 *   {
848 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
849 *   {
850 *   system_cell_matrix(i, j) +=
851 *   (viscosity *
852 *   scalar_product(grad_phi_u[i], grad_phi_u[j]) -
853 *   div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
854 *   fe_values.JxW(q);
855 *  
856 *   preconditioner_cell_matrix(i, j) += 1.0 / viscosity *
857 *   phi_p[i] * phi_p[j] *
858 *   fe_values.JxW(q);
859 *   }
860 *  
861 *   const unsigned int component_i =
862 *   fe.system_to_component_index(i).first;
863 *   cell_rhs(i) += fe_values.shape_value(i, q) *
864 *   rhs_values[q](component_i) * fe_values.JxW(q);
865 *   }
866 *   }
867 *  
868 *  
869 *   cell->get_dof_indices(local_dof_indices);
870 *   constraints.distribute_local_to_global(system_cell_matrix,
871 *   cell_rhs,
872 *   local_dof_indices,
873 *   system_matrix,
874 *   system_rhs);
875 *  
876 *   constraints.distribute_local_to_global(preconditioner_cell_matrix,
877 *   local_dof_indices,
879 *   }
880 *  
881 *   system_matrix.compress(VectorOperation::add);
884 *   }
885 *  
886 *  
887 *  
888 * @endcode
889 *
890 *
891 * <a name="step_55-Solving"></a>
892 * <h3>Solving</h3>
893 *
894
895 *
896 * This function solves the linear system with MINRES with a block diagonal
897 * preconditioner and AMG for the two diagonal blocks as described in the
898 * introduction. The preconditioner applies a v cycle to the 0,0 block
899 * and a CG with the mass matrix for the 1,1 block (the Schur complement).
900 *
901 * @code
902 *   template <int dim>
904 *   {
906 *  
907 *   LA::MPI::PreconditionAMG prec_A;
908 *   {
909 *   LA::MPI::PreconditionAMG::AdditionalData data;
910 *  
912 *   data.symmetric_operator = true;
913 *   #endif
914 *   prec_A.initialize(system_matrix.block(0, 0), data);
915 *   }
916 *  
917 *   LA::MPI::PreconditionAMG prec_S;
918 *   {
919 *   LA::MPI::PreconditionAMG::AdditionalData data;
920 *  
922 *   data.symmetric_operator = true;
923 *   #endif
924 *   prec_S.initialize(preconditioner_matrix.block(1, 1), data);
925 *   }
926 *  
927 * @endcode
928 *
929 * The InverseMatrix is used to solve for the mass matrix:
930 *
931 * @code
932 *   using mp_inverse_t = LinearSolvers::InverseMatrix<LA::MPI::SparseMatrix,
933 *   LA::MPI::PreconditionAMG>;
935 *  
936 * @endcode
937 *
938 * This constructs the block preconditioner based on the preconditioners
939 * for the individual blocks defined above.
940 *
941 * @code
942 *   const LinearSolvers::BlockDiagonalPreconditioner<LA::MPI::PreconditionAMG,
944 *   preconditioner(prec_A, mp_inverse);
945 *  
946 * @endcode
947 *
948 * With that, we can finally set up the linear solver and solve the system:
949 *
950 * @code
951 *   SolverControl solver_control(system_matrix.m(),
952 *   1e-10 * system_rhs.l2_norm());
953 *  
954 *   SolverMinRes<LA::MPI::BlockVector> solver(solver_control);
955 *  
956 *   LA::MPI::BlockVector distributed_solution(owned_partitioning,
957 *   mpi_communicator);
958 *  
959 *   constraints.set_zero(distributed_solution);
960 *  
961 *   solver.solve(system_matrix,
963 *   system_rhs,
964 *   preconditioner);
965 *  
966 *   pcout << " Solved in " << solver_control.last_step() << " iterations."
967 *   << std::endl;
968 *  
969 *   constraints.distribute(distributed_solution);
970 *  
971 * @endcode
972 *
973 * Like in @ref step_56 "step-56", we subtract the mean pressure to allow error
974 * computations against our reference solution, which has a mean value
975 * of zero.
976 *
977 * @code
979 *   const double mean_pressure =
983 *   dim);
984 *   distributed_solution.block(1).add(-mean_pressure);
986 *   }
987 *  
988 *  
989 *  
990 * @endcode
991 *
992 *
993 * <a name="step_55-Therest"></a>
994 * <h3>The rest</h3>
995 *
996
997 *
1000 *
1001 * @code
1002 *   template <int dim>
1004 *   {
1005 *   TimerOutput::Scope t(computing_timer, "refine");
1006 *  
1007 *   triangulation.refine_global();
1008 *   }
1009 *  
1010 *  
1011 *  
1012 *   template <int dim>
1013 *   void StokesProblem<dim>::output_results(const unsigned int cycle)
1014 *   {
1015 *   TimerOutput::Scope t(computing_timer, "output");
1016 *  
1017 *   {
1018 *   const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
1019 *   const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
1020 *   dim + 1);
1021 *  
1022 *   Vector<double> cellwise_errors(triangulation.n_active_cells());
1023 *   const QGauss<dim> quadrature(velocity_degree + 2);
1024 *  
1029 *   quadrature,
1031 *   &velocity_mask);
1032 *  
1033 *   const double error_u_l2 =
1037 *  
1042 *   quadrature,
1044 *   &pressure_mask);
1045 *  
1046 *   const double error_p_l2 =
1050 *  
1051 *   pcout << "error: u_0: " << error_u_l2 << " p_0: " << error_p_l2
1052 *   << std::endl;
1053 *   }
1054 *  
1055 *  
1056 *   std::vector<std::string> solution_names(dim, "velocity");
1057 *   solution_names.emplace_back("pressure");
1058 *   std::vector<DataComponentInterpretation::DataComponentInterpretation>
1059 *   data_component_interpretation(
1061 *   data_component_interpretation.push_back(
1063 *  
1064 *   DataOut<dim> data_out;
1065 *   data_out.attach_dof_handler(dof_handler);
1066 *   data_out.add_data_vector(locally_relevant_solution,
1069 *   data_component_interpretation);
1070 *  
1071 *   LA::MPI::BlockVector interpolated;
1074 *  
1075 *   LA::MPI::BlockVector interpolated_relevant(owned_partitioning,
1077 *   MPI_COMM_WORLD);
1079 *   {
1080 *   std::vector<std::string> solution_names(dim, "ref_u");
1081 *   solution_names.emplace_back("ref_p");
1082 *   data_out.add_data_vector(interpolated_relevant,
1085 *   data_component_interpretation);
1086 *   }
1087 *  
1088 *  
1089 *   Vector<float> subdomain(triangulation.n_active_cells());
1090 *   for (unsigned int i = 0; i < subdomain.size(); ++i)
1091 *   subdomain(i) = triangulation.locally_owned_subdomain();
1092 *   data_out.add_data_vector(subdomain, "subdomain");
1093 *  
1094 *   data_out.build_patches();
1095 *  
1096 *   data_out.write_vtu_with_pvtu_record(
1097 *   "./", "solution", cycle, mpi_communicator, 2);
1098 *   }
1099 *  
1100 *  
1101 *  
1102 *   template <int dim>
1104 *   {
1106 *   pcout << "Running using PETSc." << std::endl;
1107 *   #else
1108 *   pcout << "Running using Trilinos." << std::endl;
1109 *   #endif
1110 *   const unsigned int n_cycles = 5;
1111 *   for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1112 *   {
1113 *   pcout << "Cycle " << cycle << ':' << std::endl;
1114 *  
1115 *   if (cycle == 0)
1116 *   make_grid();
1117 *   else
1118 *   refine_grid();
1119 *  
1120 *   setup_system();
1121 *  
1122 *   assemble_system();
1123 *   solve();
1124 *  
1125 *   if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32)
1126 *   output_results(cycle);
1127 *  
1128 *   computing_timer.print_summary();
1129 *   computing_timer.reset();
1130 *  
1131 *   pcout << std::endl;
1132 *   }
1133 *   }
1134 *   } // namespace Step55
1135 *  
1136 *  
1137 *  
1138 *   int main(int argc, char *argv[])
1139 *   {
1140 *   try
1141 *   {
1142 *   using namespace dealii;
1143 *   using namespace Step55;
1144 *  
1146 *  
1148 *   problem.run();
1149 *   }
1150 *   catch (std::exception &exc)
1151 *   {
1152 *   std::cerr << std::endl
1153 *   << std::endl
1154 *   << "----------------------------------------------------"
1155 *   << std::endl;
1156 *   std::cerr << "Exception on processing: " << std::endl
1157 *   << exc.what() << std::endl
1158 *   << "Aborting!" << std::endl
1159 *   << "----------------------------------------------------"
1160 *   << std::endl;
1161 *  
1162 *   return 1;
1163 *   }
1164 *   catch (...)
1165 *   {
1166 *   std::cerr << std::endl
1167 *   << std::endl
1168 *   << "----------------------------------------------------"
1169 *   << std::endl;
1170 *   std::cerr << "Unknown exception!" << std::endl
1171 *   << "Aborting!" << std::endl
1172 *   << "----------------------------------------------------"
1173 *   << std::endl;
1174 *   return 1;
1175 *   }
1176 *  
1177 *   return 0;
1178 *   }
1179 * @endcode
1180<a name="step_55-Results"></a><h1>Results</h1>
1181
1182
1185
1186<table>
1187<tr>
1188 <th colspan="2">PETSc</th>
1189 <th colspan="8">number of processors</th>
1190</tr>
1191<tr>
1192 <th>cycle</th>
1193 <th>dofs</th>
1194 <th>1</th>
1195 <th>2</th>
1196 <th>4</th>
1197 <th>8</th>
1198 <th>16</th>
1199 <th>32</th>
1200 <th>64</th>
1201 <th>128</th>
1202</tr>
1203<tr>
1204 <td>0</td>
1205 <td>659</td>
1206 <td>49</td>
1207 <td>49</td>
1208 <td>49</td>
1209 <td>51</td>
1210 <td>51</td>
1211 <td>51</td>
1212 <td>49</td>
1213 <td>49</td>
1214</tr>
1215<tr>
1216 <td>1</td>
1217 <td>2467</td>
1218 <td>52</td>
1219 <td>52</td>
1220 <td>52</td>
1221 <td>52</td>
1222 <td>52</td>
1223 <td>54</td>
1224 <td>54</td>
1225 <td>53</td>
1226</tr>
1227<tr>
1228 <td>2</td>
1229 <td>9539</td>
1230 <td>56</td>
1231 <td>56</td>
1232 <td>56</td>
1233 <td>54</td>
1234 <td>56</td>
1235 <td>56</td>
1236 <td>54</td>
1237 <td>56</td>
1238</tr>
1239<tr>
1240 <td>3</td>
1241 <td>37507</td>
1242 <td>57</td>
1243 <td>57</td>
1244 <td>57</td>
1245 <td>57</td>
1246 <td>57</td>
1247 <td>56</td>
1248 <td>57</td>
1249 <td>56</td>
1250</tr>
1251<tr>
1252 <td>4</td>
1253 <td>148739</td>
1254 <td>58</td>
1255 <td>59</td>
1256 <td>57</td>
1257 <td>59</td>
1258 <td>57</td>
1259 <td>57</td>
1260 <td>57</td>
1261 <td>57</td>
1262</tr>
1263<tr>
1264 <td>5</td>
1265 <td>592387</td>
1266 <td>60</td>
1267 <td>60</td>
1268 <td>59</td>
1269 <td>59</td>
1270 <td>59</td>
1271 <td>59</td>
1272 <td>59</td>
1273 <td>59</td>
1274</tr>
1275<tr>
1276 <td>6</td>
1277 <td>2364419</td>
1278 <td>62</td>
1279 <td>62</td>
1280 <td>61</td>
1281 <td>61</td>
1282 <td>61</td>
1283 <td>61</td>
1284 <td>61</td>
1285 <td>61</td>
1286</tr>
1287</table>
1288
1289<table>
1290<tr>
1291 <th colspan="2">Trilinos</th>
1292 <th colspan="8">number of processors</th>
1293</tr>
1294<tr>
1295 <th>cycle</th>
1296 <th>dofs</th>
1297 <th>1</th>
1298 <th>2</th>
1299 <th>4</th>
1300 <th>8</th>
1301 <th>16</th>
1302 <th>32</th>
1303 <th>64</th>
1304 <th>128</th>
1305</tr>
1306<tr>
1307 <td>0</td>
1308 <td>659</td>
1309 <td>37</td>
1310 <td>37</td>
1311 <td>37</td>
1312 <td>37</td>
1313 <td>37</td>
1314 <td>37</td>
1315 <td>37</td>
1316 <td>37</td>
1317</tr>
1318<tr>
1319 <td>1</td>
1320 <td>2467</td>
1321 <td>92</td>
1322 <td>89</td>
1323 <td>89</td>
1324 <td>82</td>
1325 <td>86</td>
1326 <td>81</td>
1327 <td>78</td>
1328 <td>78</td>
1329</tr>
1330<tr>
1331 <td>2</td>
1332 <td>9539</td>
1333 <td>102</td>
1334 <td>99</td>
1335 <td>96</td>
1336 <td>95</td>
1337 <td>95</td>
1338 <td>88</td>
1339 <td>83</td>
1340 <td>95</td>
1341</tr>
1342<tr>
1343 <td>3</td>
1344 <td>37507</td>
1345 <td>107</td>
1346 <td>105</td>
1347 <td>104</td>
1348 <td>99</td>
1349 <td>100</td>
1350 <td>96</td>
1351 <td>96</td>
1352 <td>90</td>
1353</tr>
1354<tr>
1355 <td>4</td>
1356 <td>148739</td>
1357 <td>112</td>
1358 <td>112</td>
1359 <td>111</td>
1360 <td>111</td>
1361 <td>127</td>
1362 <td>126</td>
1363 <td>115</td>
1364 <td>117</td>
1365</tr>
1366<tr>
1367 <td>5</td>
1368 <td>592387</td>
1369 <td>116</td>
1370 <td>115</td>
1371 <td>114</td>
1372 <td>112</td>
1373 <td>118</td>
1374 <td>120</td>
1375 <td>131</td>
1376 <td>130</td>
1377</tr>
1378<tr>
1379 <td>6</td>
1380 <td>2364419</td>
1381 <td>130</td>
1382 <td>126</td>
1383 <td>120</td>
1384 <td>120</td>
1385 <td>121</td>
1386 <td>122</td>
1387 <td>121</td>
1388 <td>123</td>
1389</tr>
1390</table>
1391
1392While the PETSc results show a constant number of iterations, the iterations
1393increase when using Trilinos. This is likely because of the different settings
1394used for the AMG preconditioner. For performance reasons we do not allow
1396solve (we are using LU by default), a change in number of levels will
1397influence the quality of a V-cycle. Therefore, a V-cycle is closer to an exact
1398solver for smaller problem sizes.
1399
1400<a name="step-55-extensions"></a>
1401<a name="step_55-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1402
1403
1404<a name="step_55-InvestigateTrilinositerations"></a><h4>Investigate Trilinos iterations</h4>
1405
1406
1407Play with the smoothers, smoothing steps, and other properties for the
1408Trilinos AMG to achieve an optimal preconditioner.
1409
1410<a name="step_55-SolvetheOseenprobleminsteadoftheStokessystem"></a><h4>Solve the Oseen problem instead of the Stokes system</h4>
1411
1412
1415
1416You can prescribe the exact flow solution as @f$b@f$ in the convective term @f$b
1418if you set the right hand side to zero.
1419
1420<a name="step_55-Adaptiverefinement"></a><h4>Adaptive refinement</h4>
1421
1422
1424Replacing the code in StokesProblem::refine_grid() by something like
1425@code
1426Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1427
1429KellyErrorEstimator<dim>::estimate(
1430 dof_handler,
1431 QGauss<dim - 1>(fe.degree + 1),
1432 std::map<types::boundary_id, const Function<dim> *>(),
1435 fe.component_mask(velocities));
1436parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number(
1438triangulation.execute_coarsening_and_refinement();
1439@endcode
1441 *
1442 *
1443<a name="step_55-PlainProg"></a>
1444<h1> The plain program</h1>
1445@include "step-55.cc"
1446*/
constexpr void clear()
friend class Tensor
Definition tensor.h:865
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:851
@ wall_times
Definition timer.h:651
#define Assert(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Definition loop.h:564
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Definition mpi.cc:740
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
IndexSet extract_locally_relevant_dofs(const DoFHandler< dim, spacedim > &dof_handler)
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
constexpr char V
constexpr types::blas_int zero
constexpr char A
constexpr types::blas_int one
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)
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm mpi_comm, const IndexSet &locally_relevant_rows)
std::vector< unsigned int > serial(const std::vector< unsigned int > &targets, const std::function< RequestType(const unsigned int)> &create_request, const std::function< AnswerType(const unsigned int, const RequestType &)> &answer_request, const std::function< void(const unsigned int, const AnswerType &)> &process_answer, const MPI_Comm comm)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
Definition mpi.cc:97
std::vector< T > all_gather(const MPI_Comm comm, const T &object_to_send)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:112
double compute_global_error(const Triangulation< dim, spacedim > &tria, const InVector &cellwise_error, const NormType &norm, const double exponent=2.)
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={}, const unsigned int level=numbers::invalid_unsigned_int)
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const ReadVector< Number > &fe_function, const Function< spacedim, Number > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
Number compute_mean_value(const hp::MappingCollection< dim, spacedim > &mapping_collection, const DoFHandler< dim, spacedim > &dof, const hp::QCollection< dim > &q_collection, const ReadVector< Number > &v, const unsigned int component)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
constexpr double PI
Definition numbers.h:241
STL namespace.
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition types.h:32
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation