Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-72.h
Go to the documentation of this file.
1) const
491 *   {
492 *   return std::sin(2 * numbers::PI * (p[0] + p[1]));
493 *   }
494 *  
495 *  
496 * @endcode
497 *
498 *
499 * <a name="ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
500 * <h3>The <code>MinimalSurfaceProblem</code> class implementation</h3>
501 *
502
503 *
504 *
505 * <a name="MinimalSurfaceProblemMinimalSurfaceProblem"></a>
506 * <h4>MinimalSurfaceProblem::MinimalSurfaceProblem</h4>
507 *
508
509 *
510 * There have been no changes made to the class constructor.
511 *
512 * @code
513 *   template <int dim>
514 *   MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
515 *   : dof_handler(triangulation)
516 *   , fe(2)
517 *   , quadrature_formula(fe.degree + 1)
518 *   {}
519 *  
520 *  
521 * @endcode
522 *
523 *
524 * <a name="MinimalSurfaceProblemsetup_system"></a>
525 * <h4>MinimalSurfaceProblem::setup_system</h4>
526 *
527
528 *
529 * There have been no changes made to the function that sets up the class
530 * data structures, namely the DoFHandler, the hanging node constraints
531 * applied to the problem, and the linear system.
532 *
533 * @code
534 *   template <int dim>
535 *   void MinimalSurfaceProblem<dim>::setup_system(const bool initial_step)
536 *   {
537 *   if (initial_step)
538 *   {
539 *   dof_handler.distribute_dofs(fe);
540 *   current_solution.reinit(dof_handler.n_dofs());
541 *  
542 *   hanging_node_constraints.clear();
544 *   hanging_node_constraints);
545 *   hanging_node_constraints.close();
546 *   }
547 *  
548 *   newton_update.reinit(dof_handler.n_dofs());
549 *   system_rhs.reinit(dof_handler.n_dofs());
550 *  
551 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
552 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
553 *  
554 *   hanging_node_constraints.condense(dsp);
555 *  
556 *   sparsity_pattern.copy_from(dsp);
557 *   system_matrix.reinit(sparsity_pattern);
558 *   }
559 *  
560 * @endcode
561 *
562 *
563 * <a name="Assemblingthelinearsystem"></a>
564 * <h4>Assembling the linear system</h4>
565 *
566
567 *
568 *
569 * <a name="Manualassembly"></a>
570 * <h5>Manual assembly</h5>
571 *
572
573 *
574 * The assembly functions are the interesting contributions to this tutorial.
575 * The assemble_system_unassisted() method implements exactly the same
576 * assembly function as is detailed in @ref step_15 "step-15", but in this instance we
577 * use the MeshWorker::mesh_loop() function to multithread the assembly
578 * process. The reason for doing this is quite simple: When using
579 * automatic differentiation, we know that there is to be some additional
580 * computational overhead incurred. In order to mitigate this performance
581 * loss, we'd like to take advantage of as many (easily available)
582 * computational resources as possible. The MeshWorker::mesh_loop() concept
583 * makes this a relatively straightforward task. At the same time, for the
584 * purposes of fair comparison, we need to do the same to the implementation
585 * that uses no assistance when computing the residual or its linearization.
586 * (The MeshWorker::mesh_loop() function is first discussed in @ref step_12 "step-12" and
587 * @ref step_16 "step-16", if you'd like to read up on it.)
588 *
589
590 *
591 * The steps required to implement the multithreading are the same between the
592 * three functions, so we'll use the assemble_system_unassisted() function
593 * as an opportunity to focus on the multithreading itself.
594 *
595 * @code
596 *   template <int dim>
597 *   void MinimalSurfaceProblem<dim>::assemble_system_unassisted()
598 *   {
599 *   system_matrix = 0;
600 *   system_rhs = 0;
601 *  
602 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
603 *  
604 * @endcode
605 *
606 * The MeshWorker::mesh_loop() expects that we provide two exemplar data
607 * structures. The first, `ScratchData`, is to store all large data that
608 * is to be reused between threads. The `CopyData` will hold the
609 * contributions to the linear system that come from each cell. These
610 * independent matrix-vector pairs must be accumulated into the
611 * global linear system sequentially. Since we don't need anything
612 * on top of what the MeshWorker::ScratchData and MeshWorker::CopyData
613 * classes already provide, we use these exact class definitions for
614 * our problem. Note that we only require a single instance of a local
615 * matrix, local right-hand side vector, and cell degree of freedom index
616 * vector -- the MeshWorker::CopyData therefore has `1` for all three
617 * of its template arguments.
618 *
619 * @code
620 *   using ScratchData = MeshWorker::ScratchData<dim>;
621 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
622 *  
623 * @endcode
624 *
625 * We also need to know what type of iterator we'll be working with
626 * during assembly. For simplicity, we just ask the compiler to work
627 * this out for us using the decltype() specifier, knowing that we'll
628 * be iterating over active cells owned by the @p dof_handler.
629 *
630 * @code
631 *   using CellIteratorType = decltype(dof_handler.begin_active());
632 *  
633 * @endcode
634 *
635 * Here we initialize the exemplar data structures. Since we know that
636 * we need to compute the shape function gradients, weighted Jacobian,
637 * and the position of the quadrate points in real space, we pass these
638 * flags into the class constructor.
639 *
640 * @code
641 *   const ScratchData sample_scratch_data(fe,
642 *   quadrature_formula,
643 *   update_gradients |
645 *   update_JxW_values);
646 *   const CopyData sample_copy_data(dofs_per_cell);
647 *  
648 * @endcode
649 *
650 * Now we define a lambda function that will perform the assembly on
651 * a single cell. The three arguments are those that will be expected by
652 * MeshWorker::mesh_loop(), due to the arguments that we'll pass to that
653 * final call. We also capture the @p this pointer, which means that we'll
654 * have access to "this" (i.e., the current `MinimalSurfaceProblem<dim>`)
655 * class instance, and its private member data (since the lambda function is
656 * defined within a MinimalSurfaceProblem<dim> method).
657 *
658
659 *
660 * At the top of the function, we initialize the data structures
661 * that are dependent on the cell for which the work is being
662 * performed. Observe that the reinitialization call actually
663 * returns an instance to an FEValues object that is initialized
664 * and stored within (and, therefore, reused by) the
665 * `scratch_data` object.
666 *
667
668 *
669 * Similarly, we get aliases to the local matrix, local RHS
670 * vector, and local cell DoF indices from the `copy_data`
671 * instance that MeshWorker::mesh_loop() provides. We then
672 * initialize the cell DoF indices, knowing that the local matrix
673 * and vector are already correctly sized.
674 *
675 * @code
676 *   const auto cell_worker = [this](const CellIteratorType &cell,
677 *   ScratchData & scratch_data,
678 *   CopyData & copy_data) {
679 *   const auto &fe_values = scratch_data.reinit(cell);
680 *  
681 *   FullMatrix<double> & cell_matrix = copy_data.matrices[0];
682 *   Vector<double> & cell_rhs = copy_data.vectors[0];
683 *   std::vector<types::global_dof_index> &local_dof_indices =
684 *   copy_data.local_dof_indices[0];
685 *   cell->get_dof_indices(local_dof_indices);
686 *  
687 * @endcode
688 *
689 * For Newton's method, we require the gradient of the solution at the
690 * point about which the problem is being linearized.
691 *
692
693 *
694 * Once we have that, we can perform assembly for this cell in
695 * the usual way. One minor difference to @ref step_15 "step-15" is that we've
696 * used the (rather convenient) range-based loops to iterate
697 * over all quadrature points and degrees-of-freedom.
698 *
699 * @code
700 *   std::vector<Tensor<1, dim>> old_solution_gradients(
701 *   fe_values.n_quadrature_points);
702 *   fe_values.get_function_gradients(current_solution,
703 *   old_solution_gradients);
704 *  
705 *   for (const unsigned int q : fe_values.quadrature_point_indices())
706 *   {
707 *   const double coeff =
708 *   1.0 / std::sqrt(1.0 + old_solution_gradients[q] *
709 *   old_solution_gradients[q]);
710 *  
711 *   for (const unsigned int i : fe_values.dof_indices())
712 *   {
713 *   for (const unsigned int j : fe_values.dof_indices())
714 *   cell_matrix(i, j) +=
715 *   (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
716 *   * coeff // * a_n
717 *   * fe_values.shape_grad(j, q)) // * \nabla \phi_j)
718 *   - // -
719 *   (fe_values.shape_grad(i, q) // (\nabla \phi_i
720 *   * coeff * coeff * coeff // * a_n^3
721 *   * (fe_values.shape_grad(j, q) // * (\nabla \phi_j
722 *   * old_solution_gradients[q]) // * \nabla u_n)
723 *   * old_solution_gradients[q])) // * \nabla u_n)))
724 *   * fe_values.JxW(q)); // * dx
725 *  
726 *   cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
727 *   * coeff // * a_n
728 *   * old_solution_gradients[q] // * \nabla u_n
729 *   * fe_values.JxW(q)); // * dx
730 *   }
731 *   }
732 *   };
733 *  
734 * @endcode
735 *
736 * The second lambda function that MeshWorker::mesh_loop() requires is
737 * one that performs the task of accumulating the local contributions
738 * in the global linear system. That is precisely what this one does,
739 * and the details of the implementation have been seen before. The
740 * primary point to recognize is that the local contributions are stored
741 * in the `copy_data` instance that is passed into this function. This
742 * `copy_data` has been filled with data during @a some call to the
743 * `cell_worker`.
744 *
745 * @code
746 *   const auto copier = [dofs_per_cell, this](const CopyData &copy_data) {
747 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
748 *   const Vector<double> & cell_rhs = copy_data.vectors[0];
749 *   const std::vector<types::global_dof_index> &local_dof_indices =
750 *   copy_data.local_dof_indices[0];
751 *  
752 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
753 *   {
754 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
755 *   system_matrix.add(local_dof_indices[i],
756 *   local_dof_indices[j],
757 *   cell_matrix(i, j));
758 *  
759 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
760 *   }
761 *   };
762 *  
763 * @endcode
764 *
765 * We have all of the required functions definitions in place, so
766 * now we call the MeshWorker::mesh_loop() to perform the actual
767 * assembly. We pass a flag as the last parameter which states
768 * that we only want to perform the assembly on the
769 * cells. Internally, MeshWorker::mesh_loop() then distributes the
770 * available work to different threads, making efficient use of
771 * the multiple cores almost all of today's processors have to
772 * offer.
773 *
774 * @code
775 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
776 *   cell_worker,
777 *   copier,
778 *   sample_scratch_data,
779 *   sample_copy_data,
780 *   MeshWorker::assemble_own_cells);
781 *  
782 * @endcode
783 *
784 * And finally, as is done in @ref step_15 "step-15", we remove hanging nodes from the
785 * system and apply zero boundary values to the linear system that defines
786 * the Newton updates @f$\delta u^n@f$.
787 *
788 * @code
789 *   hanging_node_constraints.condense(system_matrix);
790 *   hanging_node_constraints.condense(system_rhs);
791 *  
792 *   std::map<types::global_dof_index, double> boundary_values;
793 *   VectorTools::interpolate_boundary_values(dof_handler,
794 *   0,
795 *   Functions::ZeroFunction<dim>(),
796 *   boundary_values);
797 *   MatrixTools::apply_boundary_values(boundary_values,
798 *   system_matrix,
799 *   newton_update,
800 *   system_rhs);
801 *   }
802 *  
803 * @endcode
804 *
805 *
806 * <a name="Assemblyviadifferentiationoftheresidualvector"></a>
807 * <h5>Assembly via differentiation of the residual vector</h5>
808 *
809
810 *
811 * As outlined in the introduction, what we need to do for this
812 * second approach is implement the local contributions @f$F(U)^K@f$
813 * from cell @f$K@f$ to the residual vector, and then let the
814 * AD machinery deal with how to compute the
815 * derivatives @f$J(U)_{ij}^K=\frac{\partial F(U)^K_i}{\partial U_j}@f$
816 * from it.
817 *
818
819 *
820 * For the following, recall that
821 * @f[
822 * F(U)_i^K \dealcoloneq
823 * \int\limits_K\nabla \varphi_i \cdot \left[ \frac{1}{\sqrt{1+|\nabla
824 * u|^{2}}} \nabla u \right] \, dV ,
825 * @f]
826 * where @f$u(\mathbf x)=\sum_j U_j \varphi_j(\mathbf x)@f$.
827 *
828
829 *
830 * Let us see how this is implemented in practice:
831 *
832 * @code
833 *   template <int dim>
834 *   void MinimalSurfaceProblem<dim>::assemble_system_with_residual_linearization()
835 *   {
836 *   system_matrix = 0;
837 *   system_rhs = 0;
838 *  
839 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
840 *  
841 *   using ScratchData = MeshWorker::ScratchData<dim>;
842 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
843 *   using CellIteratorType = decltype(dof_handler.begin_active());
844 *  
845 *   const ScratchData sample_scratch_data(fe,
846 *   quadrature_formula,
847 *   update_gradients |
849 *   update_JxW_values);
850 *   const CopyData sample_copy_data(dofs_per_cell);
851 *  
852 * @endcode
853 *
854 * We'll define up front the AD data structures that we'll be using,
855 * utilizing the techniques shown in @ref step_71 "step-71".
856 * In this case, we choose the helper class that will automatically compute
857 * the linearization of the finite element residual using Sacado forward
858 * automatic differentiation types. These number types can be used to
859 * compute first derivatives only. This is exactly what we want, because we
860 * know that we'll only be linearizing the residual, which means that we
861 * only need to compute first-order derivatives. The return values from the
862 * calculations are to be of type `double`.
863 *
864
865 *
866 * We also need an extractor to retrieve some data related to the field
867 * solution to the problem.
868 *
869 * @code
870 *   using ADHelper = Differentiation::AD::ResidualLinearization<
871 *   Differentiation::AD::NumberTypes::sacado_dfad,
872 *   double>;
873 *   using ADNumberType = typename ADHelper::ad_type;
874 *  
875 *   const FEValuesExtractors::Scalar u_fe(0);
876 *  
877 * @endcode
878 *
879 * With this, let us define the lambda function that will be used
880 * to compute the cell contributions to the Jacobian matrix and
881 * the right hand side:
882 *
883 * @code
884 *   const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
885 *   ScratchData & scratch_data,
886 *   CopyData & copy_data) {
887 *   const auto & fe_values = scratch_data.reinit(cell);
888 *   const unsigned int dofs_per_cell = fe_values.get_fe().n_dofs_per_cell();
889 *  
890 *   FullMatrix<double> & cell_matrix = copy_data.matrices[0];
891 *   Vector<double> & cell_rhs = copy_data.vectors[0];
892 *   std::vector<types::global_dof_index> &local_dof_indices =
893 *   copy_data.local_dof_indices[0];
894 *   cell->get_dof_indices(local_dof_indices);
895 *  
896 * @endcode
897 *
898 * We'll now create and initialize an instance of the AD helper class.
899 * To do this, we need to specify how many independent variables and
900 * dependent variables there are. The independent variables will be the
901 * number of local degrees of freedom that our solution vector has,
902 * i.e., the number @f$j@f$ in the per-element representation of the
903 * discretized solution vector
904 * @f$u (\mathbf{x})|_K = \sum\limits_{j} U^K_i \varphi_j(\mathbf{x})@f$
905 * that indicates how many solution coefficients are associated with
906 * each finite element. In deal.II, this equals
907 * FiniteElement::dofs_per_cell. The number of dependent variables will be
908 * the number of entries in the local residual vector that we will be
909 * forming. In this particular problem (like many others that employ the
910 * [standard Galerkin
911 * method](https://en.wikipedia.org/wiki/Galerkin_method)) the number of
912 * local solution coefficients matches the number of local residual
913 * equations.
914 *
915 * @code
916 *   const unsigned int n_independent_variables = local_dof_indices.size();
917 *   const unsigned int n_dependent_variables = dofs_per_cell;
918 *   ADHelper ad_helper(n_independent_variables, n_dependent_variables);
919 *  
920 * @endcode
921 *
922 * Next we inform the helper of the values of the solution, i.e., the
923 * actual values for @f$U_j@f$ about which we
924 * wish to linearize. As this is done on each element individually, we
925 * have to extract the solution coefficients from the global solution
926 * vector. In other words, we define all of those coefficients @f$U_j@f$
927 * where @f$j@f$ is a local degree of freedom as the independent variables
928 * that enter the computation of the vector @f$F(U)^{K}@f$ (the dependent
929 * function).
930 *
931
932 *
933 * Then we get the complete set of degree of freedom values as
934 * represented by auto-differentiable numbers. The operations
935 * performed with these variables are tracked by the AD library
936 * from this point until the object goes out of scope. So it is
937 * <em>precisely these variables</em> with respect to which we will
938 * compute derivatives of the residual entries.
939 *
940 * @code
941 *   ad_helper.register_dof_values(current_solution, local_dof_indices);
942 *  
943 *   const std::vector<ADNumberType> &dof_values_ad =
944 *   ad_helper.get_sensitive_dof_values();
945 *  
946 * @endcode
947 *
948 * Then we do some problem specific tasks, the first being to
949 * compute all values, (spatial) gradients, and the like based on
950 * "sensitive" AD degree of freedom values. In this instance we are
951 * retrieving the solution gradients at each quadrature point. Observe
952 * that the solution gradients are now sensitive
953 * to the values of the degrees of freedom as they use the @p ADNumberType
954 * as the scalar type and the @p dof_values_ad vector provides the local
955 * DoF values.
956 *
957 * @code
958 *   std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(
959 *   fe_values.n_quadrature_points);
960 *   fe_values[u_fe].get_function_gradients_from_local_dof_values(
961 *   dof_values_ad, old_solution_gradients);
962 *  
963 * @endcode
964 *
965 * The next variable that we declare will store the cell residual vector
966 * contributions. This is rather self-explanatory, save for one
967 * <b>very important</b> detail:
968 * Note that each entry in the vector is hand-initialized with a value
969 * of zero. This is a <em>highly recommended</em> practice, as some AD
970 * libraries appear not to safely initialize the internal data
971 * structures of these number types. Not doing so could lead to some
972 * very hard to understand or detect bugs (appreciate that the author
973 * of this program mentions this out of, generally bad, experience). So
974 * out of an abundance of caution it's worthwhile zeroing the initial
975 * value explicitly. After that, apart from a sign change the residual
976 * assembly looks much the same as we saw for the cell RHS vector before:
977 * We loop over all quadrature points, ensure that the coefficient now
978 * encodes its dependence on the (sensitive) finite element DoF values by
979 * using the correct `ADNumberType`, and finally we assemble the
980 * components of the residual vector. For complete clarity, the finite
981 * element shape functions (and their gradients, etc.) as well as the
982 * "JxW" values remain scalar
983 * valued, but the @p coeff and the @p old_solution_gradients at each
984 * quadrature point are computed in terms of the independent
985 * variables.
986 *
987 * @code
988 *   std::vector<ADNumberType> residual_ad(n_dependent_variables,
989 *   ADNumberType(0.0));
990 *   for (const unsigned int q : fe_values.quadrature_point_indices())
991 *   {
992 *   const ADNumberType coeff =
993 *   1.0 / std::sqrt(1.0 + old_solution_gradients[q] *
994 *   old_solution_gradients[q]);
995 *  
996 *   for (const unsigned int i : fe_values.dof_indices())
997 *   {
998 *   residual_ad[i] += (fe_values.shape_grad(i, q) // \nabla \phi_i
999 *   * coeff // * a_n
1000 *   * old_solution_gradients[q]) // * \nabla u_n
1001 *   * fe_values.JxW(q); // * dx
1002 *   }
1003 *   }
1004 *  
1005 * @endcode
1006 *
1007 * Once we have the full cell residual vector computed, we can register
1008 * it with the helper class.
1009 *
1010
1011 *
1012 * Thereafter, we compute the residual values (basically,
1013 * extracting the real values from what we already computed) and
1014 * their Jacobian (the linearization of each residual component
1015 * with respect to all cell DoFs) at the evaluation point. For
1016 * the purposes of assembly into the global linear system, we
1017 * have to respect the sign difference between the residual and
1018 * the RHS contribution: For Newton's method, the right hand
1019 * side vector needs to be equal to the *negative* residual
1020 * vector.
1021 *
1022 * @code
1023 *   ad_helper.register_residual_vector(residual_ad);
1024 *  
1025 *   ad_helper.compute_residual(cell_rhs);
1026 *   cell_rhs *= -1.0;
1027 *  
1028 *   ad_helper.compute_linearization(cell_matrix);
1029 *   };
1030 *  
1031 * @endcode
1032 *
1033 * The remainder of the function equals what we had previously:
1034 *
1035 * @code
1036 *   const auto copier = [dofs_per_cell, this](const CopyData &copy_data) {
1037 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1038 *   const Vector<double> & cell_rhs = copy_data.vectors[0];
1039 *   const std::vector<types::global_dof_index> &local_dof_indices =
1040 *   copy_data.local_dof_indices[0];
1041 *  
1042 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1043 *   {
1044 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1045 *   system_matrix.add(local_dof_indices[i],
1046 *   local_dof_indices[j],
1047 *   cell_matrix(i, j));
1048 *  
1049 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
1050 *   }
1051 *   };
1052 *  
1053 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
1054 *   cell_worker,
1055 *   copier,
1056 *   sample_scratch_data,
1057 *   sample_copy_data,
1059 *  
1060 *   hanging_node_constraints.condense(system_matrix);
1061 *   hanging_node_constraints.condense(system_rhs);
1062 *  
1063 *   std::map<types::global_dof_index, double> boundary_values;
1065 *   0,
1067 *   boundary_values);
1068 *   MatrixTools::apply_boundary_values(boundary_values,
1069 *   system_matrix,
1070 *   newton_update,
1071 *   system_rhs);
1072 *   }
1073 *  
1074 * @endcode
1075 *
1076 *
1077 * <a name="Assemblyviadifferentiationoftheenergyfunctional"></a>
1078 * <h5>Assembly via differentiation of the energy functional</h5>
1079 *
1080
1081 *
1082 * In this third approach, we compute residual and Jacobian as first
1083 * and second derivatives of the local energy functional
1084 * @f[
1085 * E\left( U \right)^K
1086 * \dealcoloneq \int\limits_{K} \Psi \left( u \right) \, dV
1087 * \approx \sum\limits_{q}^{n_{\textrm{q-points}}} \Psi \left( u \left(
1088 * \mathbf{X}_{q} \right) \right) \underbrace{\vert J_{q} \vert \times
1089 * W_{q}}_{\text{JxW(q)}}
1090 * @f]
1091 * with the energy density given by
1092 * @f[
1093 * \Psi \left( u \right) = \sqrt{1+|\nabla u|^{2}} .
1094 * @f]
1095 *
1096
1097 *
1098 * Let us again see how this is done:
1099 *
1100 * @code
1101 *   template <int dim>
1102 *   void MinimalSurfaceProblem<dim>::assemble_system_using_energy_functional()
1103 *   {
1104 *   system_matrix = 0;
1105 *   system_rhs = 0;
1106 *  
1107 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1108 *  
1109 *   using ScratchData = MeshWorker::ScratchData<dim>;
1110 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
1111 *   using CellIteratorType = decltype(dof_handler.begin_active());
1112 *  
1113 *   const ScratchData sample_scratch_data(fe,
1114 *   quadrature_formula,
1115 *   update_gradients |
1117 *   update_JxW_values);
1118 *   const CopyData sample_copy_data(dofs_per_cell);
1119 *  
1120 * @endcode
1121 *
1122 * In this implementation of the assembly process, we choose the helper
1123 * class that will automatically compute both the residual and its
1124 * linearization from the cell contribution to an energy functional using
1125 * nested Sacado forward automatic differentiation types.
1126 * The selected number types can be used to compute both first and
1127 * second derivatives. We require this, as the residual defined as the
1128 * sensitivity of the potential energy with respect to the DoF values (i.e.
1129 * its gradient). We'll then need to linearize the residual, implying that
1130 * second derivatives of the potential energy must be computed. You might
1131 * want to compare this with the definition of `ADHelper` used int
1132 * previous function, where we used
1133 * `Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,double>`.
1134 *
1135 * @code
1136 *   using ADHelper = Differentiation::AD::EnergyFunctional<
1137 *   Differentiation::AD::NumberTypes::sacado_dfad_dfad,
1138 *   double>;
1139 *   using ADNumberType = typename ADHelper::ad_type;
1140 *  
1141 *   const FEValuesExtractors::Scalar u_fe(0);
1142 *  
1143 * @endcode
1144 *
1145 * Let us then again define the lambda function that does the integration on
1146 * a cell.
1147 *
1148
1149 *
1150 * To initialize an instance of the helper class, we now only require
1151 * that the number of independent variables (that is, the number
1152 * of degrees of freedom associated with the element solution vector)
1153 * are known up front. This is because the second-derivative matrix that
1154 * results from an energy functional is necessarily square (and also,
1155 * incidentally, symmetric).
1156 *
1157 * @code
1158 *   const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
1159 *   ScratchData & scratch_data,
1160 *   CopyData & copy_data) {
1161 *   const auto &fe_values = scratch_data.reinit(cell);
1162 *  
1163 *   FullMatrix<double> & cell_matrix = copy_data.matrices[0];
1164 *   Vector<double> & cell_rhs = copy_data.vectors[0];
1165 *   std::vector<types::global_dof_index> &local_dof_indices =
1166 *   copy_data.local_dof_indices[0];
1167 *   cell->get_dof_indices(local_dof_indices);
1168 *  
1169 *   const unsigned int n_independent_variables = local_dof_indices.size();
1170 *   ADHelper ad_helper(n_independent_variables);
1171 *  
1172 * @endcode
1173 *
1174 * Once more, we register all cell DoFs values with the helper, followed
1175 * by extracting the "sensitive" variant of these values that are to be
1176 * used in subsequent operations that must be differentiated -- one of
1177 * those being the calculation of the solution gradients.
1178 *
1179 * @code
1180 *   ad_helper.register_dof_values(current_solution, local_dof_indices);
1181 *  
1182 *   const std::vector<ADNumberType> &dof_values_ad =
1183 *   ad_helper.get_sensitive_dof_values();
1184 *  
1185 *   std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(
1186 *   fe_values.n_quadrature_points);
1187 *   fe_values[u_fe].get_function_gradients_from_local_dof_values(
1188 *   dof_values_ad, old_solution_gradients);
1189 *  
1190 * @endcode
1191 *
1192 * We next create a variable that stores the cell total energy.
1193 * Once more we emphasize that we explicitly zero-initialize this value,
1194 * thereby ensuring the integrity of the data for this starting value.
1195 *
1196
1197 *
1198 * The aim for our approach is then to compute the cell total
1199 * energy, which is the sum of the internal (due to right hand
1200 * side functions, typically linear in @f$U@f$) and external
1201 * energies. In this particular case, we have no external
1202 * energies (e.g., from source terms or Neumann boundary
1203 * conditions), so we'll focus on the internal energy part.
1204 *
1205
1206 *
1207 * In fact, computing @f$E(U)^K@f$ is almost trivial, requiring only
1208 * the following lines:
1209 *
1210 * @code
1211 *   ADNumberType energy_ad = ADNumberType(0.0);
1212 *   for (const unsigned int q : fe_values.quadrature_point_indices())
1213 *   {
1214 *   const ADNumberType psi = std::sqrt(1.0 + old_solution_gradients[q] *
1215 *   old_solution_gradients[q]);
1216 *  
1217 *   energy_ad += psi * fe_values.JxW(q);
1218 *   }
1219 *  
1220 * @endcode
1221 *
1222 * After we've computed the total energy on this cell, we'll
1223 * register it with the helper. Based on that, we may now
1224 * compute the desired quantities, namely the residual values
1225 * and their Jacobian at the evaluation point. As before, the
1226 * Newton right hand side needs to be the negative of the
1227 * residual:
1228 *
1229 * @code
1230 *   ad_helper.register_energy_functional(energy_ad);
1231 *  
1232 *   ad_helper.compute_residual(cell_rhs);
1233 *   cell_rhs *= -1.0;
1234 *  
1235 *   ad_helper.compute_linearization(cell_matrix);
1236 *   };
1237 *  
1238 * @endcode
1239 *
1240 * As in the previous two functions, the remainder of the function is as
1241 * before:
1242 *
1243 * @code
1244 *   const auto copier = [dofs_per_cell, this](const CopyData &copy_data) {
1245 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1246 *   const Vector<double> & cell_rhs = copy_data.vectors[0];
1247 *   const std::vector<types::global_dof_index> &local_dof_indices =
1248 *   copy_data.local_dof_indices[0];
1249 *  
1250 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1251 *   {
1252 *   for (unsigned int j = 0; j < dofs_per_cell; ++j)
1253 *   system_matrix.add(local_dof_indices[i],
1254 *   local_dof_indices[j],
1255 *   cell_matrix(i, j));
1256 *  
1257 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
1258 *   }
1259 *   };
1260 *  
1261 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
1262 *   cell_worker,
1263 *   copier,
1264 *   sample_scratch_data,
1265 *   sample_copy_data,
1267 *  
1268 *   hanging_node_constraints.condense(system_matrix);
1269 *   hanging_node_constraints.condense(system_rhs);
1270 *  
1271 *   std::map<types::global_dof_index, double> boundary_values;
1273 *   0,
1275 *   boundary_values);
1276 *   MatrixTools::apply_boundary_values(boundary_values,
1277 *   system_matrix,
1278 *   newton_update,
1279 *   system_rhs);
1280 *   }
1281 *  
1282 *  
1283 * @endcode
1284 *
1285 *
1286 * <a name="MinimalSurfaceProblemsolve"></a>
1287 * <h4>MinimalSurfaceProblem::solve</h4>
1288 *
1289
1290 *
1291 * The solve function is the same as is used in @ref step_15 "step-15".
1292 *
1293 * @code
1294 *   template <int dim>
1295 *   void MinimalSurfaceProblem<dim>::solve()
1296 *   {
1297 *   SolverControl solver_control(system_rhs.size(),
1298 *   system_rhs.l2_norm() * 1e-6);
1299 *   SolverCG<Vector<double>> solver(solver_control);
1300 *  
1301 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1302 *   preconditioner.initialize(system_matrix, 1.2);
1303 *  
1304 *   solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
1305 *  
1306 *   hanging_node_constraints.distribute(newton_update);
1307 *  
1308 *   const double alpha = determine_step_length();
1309 *   current_solution.add(alpha, newton_update);
1310 *   }
1311 *  
1312 *  
1313 * @endcode
1314 *
1315 *
1316 * <a name="MinimalSurfaceProblemrefine_mesh"></a>
1317 * <h4>MinimalSurfaceProblem::refine_mesh</h4>
1318 *
1319
1320 *
1321 * Nothing has changed since @ref step_15 "step-15" with respect to the mesh refinement
1322 * procedure and transfer of the solution between adapted meshes.
1323 *
1324 * @code
1325 *   template <int dim>
1326 *   void MinimalSurfaceProblem<dim>::refine_mesh()
1327 *   {
1328 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1329 *  
1331 *   dof_handler,
1332 *   QGauss<dim - 1>(fe.degree + 1),
1333 *   std::map<types::boundary_id, const Function<dim> *>(),
1334 *   current_solution,
1335 *   estimated_error_per_cell);
1336 *  
1338 *   estimated_error_per_cell,
1339 *   0.3,
1340 *   0.03);
1341 *  
1342 *   triangulation.prepare_coarsening_and_refinement();
1343 *   SolutionTransfer<dim> solution_transfer(dof_handler);
1344 *   solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
1345 *   triangulation.execute_coarsening_and_refinement();
1346 *  
1347 *   dof_handler.distribute_dofs(fe);
1348 *  
1349 *   Vector<double> tmp(dof_handler.n_dofs());
1350 *   solution_transfer.interpolate(current_solution, tmp);
1351 *   current_solution = tmp;
1352 *  
1353 *   hanging_node_constraints.clear();
1355 *   hanging_node_constraints);
1356 *   hanging_node_constraints.close();
1357 *  
1358 *   set_boundary_values();
1359 *  
1360 *   setup_system(false);
1361 *   }
1362 *  
1363 *  
1364 *  
1365 * @endcode
1366 *
1367 *
1368 * <a name="MinimalSurfaceProblemset_boundary_values"></a>
1369 * <h4>MinimalSurfaceProblem::set_boundary_values</h4>
1370 *
1371
1372 *
1373 * The choice of boundary conditions remains identical to @ref step_15 "step-15"...
1374 *
1375 * @code
1376 *   template <int dim>
1377 *   void MinimalSurfaceProblem<dim>::set_boundary_values()
1378 *   {
1379 *   std::map<types::global_dof_index, double> boundary_values;
1381 *   0,
1382 *   BoundaryValues<dim>(),
1383 *   boundary_values);
1384 *   for (auto &boundary_value : boundary_values)
1385 *   current_solution(boundary_value.first) = boundary_value.second;
1386 *  
1387 *   hanging_node_constraints.distribute(current_solution);
1388 *   }
1389 *  
1390 *  
1391 * @endcode
1392 *
1393 *
1394 * <a name="MinimalSurfaceProblemcompute_residual"></a>
1395 * <h4>MinimalSurfaceProblem::compute_residual</h4>
1396 *
1397
1398 *
1399 * ... as does the function used to compute the residual during the
1400 * solution iteration procedure. One could replace this by
1401 * differentiation of the energy functional if one really wanted,
1402 * but for simplicity we here simply copy what we already had in
1403 * @ref step_15 "step-15".
1404 *
1405 * @code
1406 *   template <int dim>
1407 *   double MinimalSurfaceProblem<dim>::compute_residual(const double alpha) const
1408 *   {
1409 *   Vector<double> residual(dof_handler.n_dofs());
1410 *  
1411 *   Vector<double> evaluation_point(dof_handler.n_dofs());
1412 *   evaluation_point = current_solution;
1413 *   evaluation_point.add(alpha, newton_update);
1414 *  
1415 *   FEValues<dim> fe_values(fe,
1416 *   quadrature_formula,
1418 *   update_JxW_values);
1419 *  
1420 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1421 *   const unsigned int n_q_points = quadrature_formula.size();
1422 *  
1423 *   Vector<double> cell_residual(dofs_per_cell);
1424 *   std::vector<Tensor<1, dim>> gradients(n_q_points);
1425 *  
1426 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1427 *  
1428 *   for (const auto &cell : dof_handler.active_cell_iterators())
1429 *   {
1430 *   cell_residual = 0;
1431 *   fe_values.reinit(cell);
1432 *  
1433 *   fe_values.get_function_gradients(evaluation_point, gradients);
1434 *  
1435 *   for (unsigned int q = 0; q < n_q_points; ++q)
1436 *   {
1437 *   const double coeff =
1438 *   1.0 / std::sqrt(1.0 + gradients[q] * gradients[q]);
1439 *  
1440 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1441 *   cell_residual(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
1442 *   * coeff // * a_n
1443 *   * gradients[q] // * \nabla u_n
1444 *   * fe_values.JxW(q)); // * dx
1445 *   }
1446 *  
1447 *   cell->get_dof_indices(local_dof_indices);
1448 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1449 *   residual(local_dof_indices[i]) += cell_residual(i);
1450 *   }
1451 *  
1452 *   hanging_node_constraints.condense(residual);
1453 *  
1454 *   for (const types::global_dof_index i :
1455 *   DoFTools::extract_boundary_dofs(dof_handler))
1456 *   residual(i) = 0;
1457 *  
1458 *   return residual.l2_norm();
1459 *   }
1460 *  
1461 *  
1462 *  
1463 * @endcode
1464 *
1465 *
1466 * <a name="MinimalSurfaceProblemdetermine_step_length"></a>
1467 * <h4>MinimalSurfaceProblem::determine_step_length</h4>
1468 *
1469
1470 *
1471 * The choice of step length (or, under-relaxation factor) for the nonlinear
1472 * iterations procedure remains fixed at the value chosen and discussed in
1473 * @ref step_15 "step-15".
1474 *
1475 * @code
1476 *   template <int dim>
1477 *   double MinimalSurfaceProblem<dim>::determine_step_length() const
1478 *   {
1479 *   return 0.1;
1480 *   }
1481 *  
1482 *  
1483 *  
1484 * @endcode
1485 *
1486 *
1487 * <a name="MinimalSurfaceProblemoutput_results"></a>
1488 * <h4>MinimalSurfaceProblem::output_results</h4>
1489 *
1490
1491 *
1492 * This last function to be called from `run()` outputs the current solution
1493 * (and the Newton update) in graphical form as a VTU file. It is entirely the
1494 * same as what has been used in previous tutorials.
1495 *
1496 * @code
1497 *   template <int dim>
1498 *   void MinimalSurfaceProblem<dim>::output_results(
1499 *   const unsigned int refinement_cycle) const
1500 *   {
1501 *   DataOut<dim> data_out;
1502 *  
1503 *   data_out.attach_dof_handler(dof_handler);
1504 *   data_out.add_data_vector(current_solution, "solution");
1505 *   data_out.add_data_vector(newton_update, "update");
1506 *   data_out.build_patches();
1507 *  
1508 *   const std::string filename =
1509 *   "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
1510 *   std::ofstream output(filename);
1511 *   data_out.write_vtu(output);
1512 *   }
1513 *  
1514 *  
1515 * @endcode
1516 *
1517 *
1518 * <a name="MinimalSurfaceProblemrun"></a>
1519 * <h4>MinimalSurfaceProblem::run</h4>
1520 *
1521
1522 *
1523 * In the run function, most remains the same as was first implemented
1524 * in @ref step_15 "step-15". The only observable changes are that we can now choose (via
1525 * the parameter file) what the final acceptable tolerance for the system
1526 * residual is, and that we can choose which method of assembly we wish to
1527 * utilize. To make the second choice clear, we output to the console some
1528 * message which indicates the selection. Since we're interested in comparing
1529 * the time taken to assemble for each of the three methods, we've also
1530 * added a timer that keeps a track of how much time is spent during assembly.
1531 * We also track the time taken to solve the linear system, so that we can
1532 * contrast those numbers to the part of the code which would normally take
1533 * the longest time to execute.
1534 *
1535 * @code
1536 *   template <int dim>
1537 *   void MinimalSurfaceProblem<dim>::run(const int formulation,
1538 *   const double tolerance)
1539 *   {
1540 *   std::cout << "******** Assembly approach ********" << std::endl;
1541 *   const std::array<std::string, 3> method_descriptions = {
1542 *   {"Unassisted implementation (full hand linearization).",
1543 *   "Automated linearization of the finite element residual.",
1544 *   "Automated computation of finite element residual and linearization using a variational formulation."}};
1545 *   AssertIndexRange(formulation, method_descriptions.size());
1546 *   std::cout << method_descriptions[formulation] << std::endl << std::endl;
1547 *  
1548 *  
1550 *  
1552 *   triangulation.refine_global(2);
1553 *  
1554 *   setup_system(/*first time=*/true);
1555 *   set_boundary_values();
1556 *  
1557 *   double last_residual_norm = std::numeric_limits<double>::max();
1558 *   unsigned int refinement_cycle = 0;
1559 *   do
1560 *   {
1561 *   std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
1562 *  
1563 *   if (refinement_cycle != 0)
1564 *   refine_mesh();
1565 *  
1566 *   std::cout << " Initial residual: " << compute_residual(0) << std::endl;
1567 *  
1568 *   for (unsigned int inner_iteration = 0; inner_iteration < 5;
1569 *   ++inner_iteration)
1570 *   {
1571 *   {
1572 *   TimerOutput::Scope t(timer, "Assemble");
1573 *  
1574 *   if (formulation == 0)
1575 *   assemble_system_unassisted();
1576 *   else if (formulation == 1)
1577 *   assemble_system_with_residual_linearization();
1578 *   else if (formulation == 2)
1579 *   assemble_system_using_energy_functional();
1580 *   else
1581 *   AssertThrow(false, ExcNotImplemented());
1582 *   }
1583 *  
1584 *   last_residual_norm = system_rhs.l2_norm();
1585 *  
1586 *   {
1587 *   TimerOutput::Scope t(timer, "Solve");
1588 *   solve();
1589 *   }
1590 *  
1591 *  
1592 *   std::cout << " Residual: " << compute_residual(0) << std::endl;
1593 *   }
1594 *  
1595 *   output_results(refinement_cycle);
1596 *  
1597 *   ++refinement_cycle;
1598 *   std::cout << std::endl;
1599 *   }
1600 *   while (last_residual_norm > tolerance);
1601 *   }
1602 *   } // namespace Step72
1603 *  
1604 * @endcode
1605 *
1606 *
1607 * <a name="Themainfunction"></a>
1608 * <h4>The main function</h4>
1609 *
1610
1611 *
1612 * Finally the main function. This follows the scheme of most other main
1613 * functions, with two obvious exceptions:
1614 * - We call Utilities::MPI::MPI_InitFinalize in order to set up (via a hidden
1615 * default parameter) the number of threads using the execution of
1616 * multithreaded tasks.
1617 * - We also have a few lines dedicates to reading in or initializing the
1618 * user-defined parameters that will be considered during the execution of the
1619 * program.
1620 *
1621 * @code
1622 *   int main(int argc, char *argv[])
1623 *   {
1624 *   try
1625 *   {
1626 *   using namespace Step72;
1627 *  
1628 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
1629 *  
1630 *   std::string prm_file;
1631 *   if (argc > 1)
1632 *   prm_file = argv[1];
1633 *   else
1634 *   prm_file = "parameters.prm";
1635 *  
1636 *   const MinimalSurfaceProblemParameters parameters;
1637 *   ParameterAcceptor::initialize(prm_file);
1638 *  
1639 *   MinimalSurfaceProblem<2> minimal_surface_problem_2d;
1640 *   minimal_surface_problem_2d.run(parameters.formulation,
1641 *   parameters.tolerance);
1642 *   }
1643 *   catch (std::exception &exc)
1644 *   {
1645 *   std::cerr << std::endl
1646 *   << std::endl
1647 *   << "----------------------------------------------------"
1648 *   << std::endl;
1649 *   std::cerr << "Exception on processing: " << std::endl
1650 *   << exc.what() << std::endl
1651 *   << "Aborting!" << std::endl
1652 *   << "----------------------------------------------------"
1653 *   << std::endl;
1654 *  
1655 *   return 1;
1656 *   }
1657 *   catch (...)
1658 *   {
1659 *   std::cerr << std::endl
1660 *   << std::endl
1661 *   << "----------------------------------------------------"
1662 *   << std::endl;
1663 *   std::cerr << "Unknown exception!" << std::endl
1664 *   << "Aborting!" << std::endl
1665 *   << "----------------------------------------------------"
1666 *   << std::endl;
1667 *   return 1;
1668 *   }
1669 *   return 0;
1670 *   }
1671 * @endcode
1672<a name="Results"></a><h1>Results</h1>
1673
1674
1675Since there was no change to the physics of the problem that has first been analyzed
1676in @ref step_15 "step-15", there is nothing to report about that. The only outwardly noticeable
1677difference between them is that, by default, this program will only run 9 mesh
1678refinement steps (as opposed to @ref step_15 "step-15", which executes 11 refinements).
1679This will be observable in the simulation status that appears between the
1680header text that prints which assembly method is being used, and the final
1681timings. (All timings reported below were obtained in release mode.)
1682
1683@code
1684Mesh refinement step 0
1685 Initial residual: 1.53143
1686 Residual: 1.08746
1687 Residual: 0.966748
1688 Residual: 0.859602
1689 Residual: 0.766462
1690 Residual: 0.685475
1691
1692...
1693
1694Mesh refinement step 9
1695 Initial residual: 0.00924594
1696 Residual: 0.00831928
1697 Residual: 0.0074859
1698 Residual: 0.0067363
1699 Residual: 0.00606197
1700 Residual: 0.00545529
1701@endcode
1702
1703So what is interesting for us to compare is how long the assembly process takes
1704for the three different implementations, and to put that into some greater context.
1705Below is the output for the hand linearization (as computed on a circa 2012
1706four core, eight thread laptop -- but we're only really interested in the
1707relative time between the different implementations):
1708@code
1709******** Assembly approach ********
1710Unassisted implementation (full hand linearization).
1711
1712...
1713
1714+---------------------------------------------+------------+------------+
1715| Total wallclock time elapsed since start | 35.1s | |
1716| | | |
1717| Section | no. calls | wall time | % of total |
1718+---------------------------------+-----------+------------+------------+
1719| Assemble | 50 | 1.56s | 4.5% |
1720| Solve | 50 | 30.8s | 88% |
1721+---------------------------------+-----------+------------+------------+
1722@endcode
1723And for the implementation that linearizes the residual in an automated
1724manner using the Sacado dynamic forward AD number type:
1725@code
1726******** Assembly approach ********
1727Automated linearization of the finite element residual.
1728
1729...
1730
1731+---------------------------------------------+------------+------------+
1732| Total wallclock time elapsed since start | 40.1s | |
1733| | | |
1734| Section | no. calls | wall time | % of total |
1735+---------------------------------+-----------+------------+------------+
1736| Assemble | 50 | 8.8s | 22% |
1737| Solve | 50 | 28.6s | 71% |
1738+---------------------------------+-----------+------------+------------+
1739@endcode
1740And, lastly, for the implementation that computes both the residual and
1741its linearization directly from an energy functional (using nested Sacado
1742dynamic forward AD numbers):
1743@code
1744******** Assembly approach ********
1745Automated computation of finite element residual and linearization using a variational formulation.
1746
1747...
1748
1749+---------------------------------------------+------------+------------+
1750| Total wallclock time elapsed since start | 48.8s | |
1751| | | |
1752| Section | no. calls | wall time | % of total |
1753+---------------------------------+-----------+------------+------------+
1754| Assemble | 50 | 16.7s | 34% |
1755| Solve | 50 | 29.3s | 60% |
1756+---------------------------------+-----------+------------+------------+
1757@endcode
1758
1759It's evident that the more work that is passed off to the automatic differentiation
1760framework to perform, the more time is spent during the assembly process. Accumulated
1761over all refinement steps, using one level of automatic differentiation resulted
1762in @f$5.65 \times@f$ more computational time being spent in the assembly stage when
1763compared to unassisted assembly, while assembling the discrete linear system took
1764@f$10.7 \times@f$ longer when deriving directly from the energy functional.
1765Unsurprisingly, the overall time spent solving the linear system remained unchanged.
1766This means that the proportion of time spent in the solve phase to the assembly phase
1767shifted significantly as the number of times automated differentiation was performed
1768at the finite element level. For many, this might mean that leveraging higher-order
1769differentiation (at the finite element level) in production code leads to an
1770unacceptable overhead, but it may still be useful during the prototyping phase.
1771A good compromise between the two may, therefore, be the automated linearization
1772of the finite element residual, which offers a lot of convenience at a measurable,
1773but perhaps not unacceptable, cost. Alternatively, one could consider
1774not re-building the Newton matrix in every step -- a topic that is
1775explored in substantial depth in @ref step_77 "step-77".
1776
1777Of course, in practice the actual overhead is very much dependent on the problem being evaluated
1778(e.g., how many components there are in the solution, what the nature of the function
1779being differentiated is, etc.). So the exact results presented here should be
1780interpreted within the context of this scalar problem alone, and when it comes to
1781other problems, some preliminary investigation by the user is certainly warranted.
1782
1783
1784<a name="Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1785
1786
1787Like @ref step_71 "step-71", there are a few items related to automatic differentiation that could
1788be evaluated further:
1789- The use of other AD frameworks should be investigated, with the outlook that
1790 alternative implementations may provide performance benefits.
1791- It is also worth evaluating AD number types other than those that have been
1792 hard-coded into this tutorial. With regard to twice differentiable types
1793 employed at the finite-element level, mixed differentiation modes ("RAD-FAD")
1794 should in principle be more computationally efficient than the single
1795 mode ("FAD-FAD") types employed here. The reason that the RAD-FAD type was not
1796 selected by default is that, at the time of writing, there remain some
1797 bugs in its implementation within the Sacado library that lead to memory leaks.
1798 This is documented in the @ref auto_symb_diff module.
1799- It might be the case that using reduced precision types (i.e., `float`) as the
1800 scalar types for the AD numbers could render a reduction in computational
1801 expense during assembly. Using `float` as the data type for the
1802 matrix and the residual is not unreasonable, given that the Newton
1803 update is only meant to get us closer to the solution, but not
1804 actually *to* the solution; as a consequence, it makes sense to
1805 consider using reduced-precision data types for computing these
1806 updates, and then accumulating these updates in a solution vector
1807 that uses the full `double` precision accuracy.
1808- One further method of possibly reducing resources during assembly is to frame
1809 the AD implementations as a constitutive model. This would be similar to the
1810 approach adopted in @ref step_71 "step-71", and pushes the starting point for the automatic
1811 differentiation one level higher up the chain of computations. This, in turn,
1812 means that less operations are tracked by the AD library, thereby reducing the
1813 cost of differentiating (even though one would perform the differentiation at
1814 each cell quadrature point).
1815- @ref step_77 "step-77" is yet another variation of @ref step_15 "step-15" that addresses a very
1816 different part of the problem: Line search and whether it is
1817 necessary to re-build the Newton matrix in every nonlinear
1818 iteration. Given that the results above show that using automatic
1819 differentiation comes at a cost, the techniques in @ref step_77 "step-77" have the
1820 potential to offset these costs to some degree. It would therefore
1821 be quite interesting to combine the current program with the
1822 techniques in @ref step_77 "step-77". For production codes, this would certainly be
1823 the way to go.
1824 *
1825 *
1826<a name="PlainProg"></a>
1827<h1> The plain program</h1>
1828@include "step-72.cc"
1829*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
const unsigned int dofs_per_cell
Definition fe_data.h:437
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, typename InputVector::value_type > * > &neumann_bc, const InputVector &solution, Vector< float > &error, const ComponentMask &component_mask=ComponentMask(), const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
@ wall_times
Definition timer.h:652
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
float depth
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
__global__ void reduction(Number *result, const Number *v, const size_type N)
__global__ void set(Number *val, const Number s, const size_type N)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:282
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void extract_boundary_dofs(const DoFHandler< dim, spacedim > &dof_handler, const ComponentMask &component_mask, std::vector< bool > &selected_dofs, const std::set< types::boundary_id > &boundary_ids={})
Definition dof_tools.cc:545
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
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
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Definition advection.h:131
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
constexpr const ReferenceCell Line
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:471
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
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)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
static constexpr double PI
Definition numbers.h:259
STL namespace.
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition types.h:33
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation