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