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