deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-72.h
Go to the documentation of this file.
1) const
489 *   {
490 *   return std::sin(2 * numbers::PI * (p[0] + p[1]));
491 *   }
492 *  
493 *  
494 * @endcode
495 *
496 *
497 * <a name="step_72-ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
499 *
500
501 *
502 *
503 * <a name="step_72-MinimalSurfaceProblemMinimalSurfaceProblem"></a>
505 *
506
507 *
509 *
510 * @code
511 *   template <int dim>
513 *   : dof_handler(triangulation)
514 *   , fe(2)
515 *   , quadrature_formula(fe.degree + 1)
516 *   {}
517 *  
518 *  
519 * @endcode
520 *
521 *
522 * <a name="step_72-MinimalSurfaceProblemsetup_system"></a>
524 *
525
526 *
527 * There have been no changes made to the function that sets up the class
529 * applied to the problem, and the linear system.
530 *
531 * @code
532 *   template <int dim>
534 *   {
535 *   dof_handler.distribute_dofs(fe);
536 *   current_solution.reinit(dof_handler.n_dofs());
537 *   newton_update.reinit(dof_handler.n_dofs());
538 *   system_rhs.reinit(dof_handler.n_dofs());
539 *  
542 *   0,
546 *   zero_constraints.close();
547 *  
550 *   0,
553 *   nonzero_constraints.close();
554 *  
555 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
557 *  
558 *   sparsity_pattern.copy_from(dsp);
559 *   system_matrix.reinit(sparsity_pattern);
560 *   }
561 *  
562 * @endcode
563 *
564 *
565 * <a name="step_72-Assemblingthelinearsystem"></a>
566 * <h4>Assembling the linear system</h4>
567 *
568
569 *
570 *
571 * <a name="step_72-Manualassembly"></a>
572 * <h5>Manual assembly</h5>
573 *
574
575 *
578 * assembly function as is detailed in @ref step_15 "step-15", but in this instance we
579 * use the MeshWorker::mesh_loop() function to multithread the assembly
580 * process. The reason for doing this is quite simple: When using
588 * (The MeshWorker::mesh_loop() function is first discussed in @ref step_12 "step-12" and
589 * @ref step_16 "step-16", if you'd like to read up on it.)
590 *
591
592 *
594 * three functions, so we'll use the assemble_system_unassisted() function
596 *
597 * @code
598 *   template <int dim>
600 *   {
601 *   system_matrix = 0;
602 *   system_rhs = 0;
603 *  
604 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
605 *  
606 * @endcode
607 *
609 * structures. The first, `ScratchData`, is to store all large data that
610 * is to be reused between threads. The `CopyData` will hold the
612 * independent matrix-vector pairs must be accumulated into the
613 * global linear system sequentially. Since we don't need anything
614 * on top of what the MeshWorker::ScratchData and MeshWorker::CopyData
617 * matrix, local right-hand side vector, and cell degree of freedom index
618 * vector -- the MeshWorker::CopyData therefore has `1` for all three
619 * of its template arguments.
620 *
621 * @code
622 *   using ScratchData = MeshWorker::ScratchData<dim>;
623 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
624 *  
625 * @endcode
626 *
627 * We also need to know what type of iterator we'll be working with
630 * be iterating over active cells owned by the @p dof_handler.
631 *
632 * @code
633 *   using CellIteratorType = decltype(dof_handler.begin_active());
634 *  
635 * @endcode
636 *
637 * Here we initialize the exemplar data structures. Since we know that
638 * we need to compute the shape function gradients, weighted Jacobian,
639 * and the position of the quadrate points in real space, we pass these
640 * flags into the class constructor.
641 *
642 * @code
643 *   const ScratchData sample_scratch_data(fe,
648 *   const CopyData sample_copy_data(dofs_per_cell);
649 *  
650 * @endcode
651 *
652 * Now we define a lambda function that will perform the assembly on
653 * a single cell. The three arguments are those that will be expected by
654 * MeshWorker::mesh_loop(), due to the arguments that we'll pass to that
655 * final call. We also capture the @p this pointer, which means that we'll
656 * have access to "this" (i.e., the current `MinimalSurfaceProblem<dim>`)
657 * class instance, and its private member data (since the lambda function is
658 * defined within a MinimalSurfaceProblem<dim> method).
659 *
660
661 *
662 * At the top of the function, we initialize the data structures
663 * that are dependent on the cell for which the work is being
665 * returns an instance to an FEValues object that is initialized
667 * `scratch_data` object.
668 *
669
670 *
671 * Similarly, we get aliases to the local matrix, local RHS
672 * vector, and local cell DoF indices from the `copy_data`
673 * instance that MeshWorker::mesh_loop() provides. We then
674 * initialize the cell DoF indices, knowing that the local matrix
675 * and vector are already correctly sized.
676 *
677 * @code
678 *   const auto cell_worker = [this](const CellIteratorType &cell,
679 *   ScratchData &scratch_data,
680 *   CopyData &copy_data) {
681 *   const auto &fe_values = scratch_data.reinit(cell);
682 *  
683 *   FullMatrix<double> &cell_matrix = copy_data.matrices[0];
684 *   Vector<double> &cell_rhs = copy_data.vectors[0];
685 *   std::vector<types::global_dof_index> &local_dof_indices =
686 *   copy_data.local_dof_indices[0];
687 *   cell->get_dof_indices(local_dof_indices);
688 *  
689 * @endcode
690 *
691 * For Newton's method, we require the gradient of the solution at the
692 * point about which the problem is being linearized.
693 *
694
695 *
696 * Once we have that, we can perform assembly for this cell in
697 * the usual way. One minor difference to @ref step_15 "step-15" is that we've
698 * used the (rather convenient) range-based loops to iterate
699 * over all quadrature points and degrees-of-freedom.
700 *
701 * @code
702 *   std::vector<Tensor<1, dim>> old_solution_gradients(
703 *   fe_values.n_quadrature_points);
704 *   fe_values.get_function_gradients(current_solution,
706 *  
707 *   for (const unsigned int q : fe_values.quadrature_point_indices())
708 *   {
709 *   const double coeff =
710 *   1.0 / std::sqrt(1.0 + old_solution_gradients[q] *
712 *  
713 *   for (const unsigned int i : fe_values.dof_indices())
714 *   {
715 *   for (const unsigned int j : fe_values.dof_indices())
716 *   cell_matrix(i, j) +=
717 *   (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
718 *   * coeff // * a_n
719 *   * fe_values.shape_grad(j, q)) // * \nabla \phi_j)
720 *   - // -
721 *   (fe_values.shape_grad(i, q) // (\nabla \phi_i
722 *   * coeff * coeff * coeff // * a_n^3
723 *   * (fe_values.shape_grad(j, q) // * (\nabla \phi_j
724 *   * old_solution_gradients[q]) // * \nabla u_n)
725 *   * old_solution_gradients[q])) // * \nabla u_n)))
726 *   * fe_values.JxW(q)); // * dx
727 *  
728 *   cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
729 *   * coeff // * a_n
730 *   * old_solution_gradients[q] // * \nabla u_n
731 *   * fe_values.JxW(q)); // * dx
732 *   }
733 *   }
734 *   };
735 *  
736 * @endcode
737 *
738 * The second lambda function that MeshWorker::mesh_loop() requires is
740 * in the global linear system. That is precisely what this one does,
743 * in the `copy_data` instance that is passed into this function. This
744 * `copy_data` has been filled with data during @a some call to the
745 * `cell_worker`.
746 *
747 * @code
748 *   const auto copier = [this](const CopyData &copy_data) {
749 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
750 *   const Vector<double> &cell_rhs = copy_data.vectors[0];
751 *   const std::vector<types::global_dof_index> &local_dof_indices =
752 *   copy_data.local_dof_indices[0];
753 *  
754 *   zero_constraints.distribute_local_to_global(
755 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
756 *   };
757 *  
758 * @endcode
759 *
762 * assembly. We pass a flag as the last parameter which states
764 * cells. Internally, MeshWorker::mesh_loop() then distributes the
765 * available work to different threads, making efficient use of
767 * offer.
768 *
769 * @code
770 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
771 *   cell_worker,
772 *   copier,
773 *   sample_scratch_data,
774 *   sample_copy_data,
775 *   MeshWorker::assemble_own_cells);
776 *   }
777 *  
778 * @endcode
779 *
780 *
781 * <a name="step_72-Assemblyviadifferentiationoftheresidualvector"></a>
782 * <h5>Assembly via differentiation of the residual vector</h5>
783 *
784
785 *
788 * from cell @f$K@f$ to the residual vector, and then let the
789 * AD machinery deal with how to compute the
790 * derivatives @f$J(U)_{ij}^K=\frac{\partial F(U)^K_i}{\partial U_j}@f$
791 * from it.
792 *
793
794 *
796 * @f[
799 * u|^{2}}} \nabla u \right] \, dV ,
800 * @f]
802 *
803
804 *
805 * Let us see how this is implemented in practice:
806 *
807 * @code
808 *   template <int dim>
810 *   {
811 *   system_matrix = 0;
812 *   system_rhs = 0;
813 *  
814 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
815 *  
816 *   using ScratchData = MeshWorker::ScratchData<dim>;
817 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
818 *   using CellIteratorType = decltype(dof_handler.begin_active());
819 *  
820 *   const ScratchData sample_scratch_data(fe,
825 *   const CopyData sample_copy_data(dofs_per_cell);
826 *  
827 * @endcode
828 *
829 * We'll define up front the AD data structures that we'll be using,
830 * utilizing the techniques shown in @ref step_71 "step-71".
831 * In this case, we choose the helper class that will automatically compute
832 * the linearization of the finite element residual using Sacado forward
833 * automatic differentiation types. These number types can be used to
835 * know that we'll only be linearizing the residual, which means that we
836 * only need to compute first-order derivatives. The return values from the
837 * calculations are to be of type `double`.
838 *
839
840 *
841 * We also need an extractor to retrieve some data related to the field
842 * solution to the problem.
843 *
844 * @code
845 *   using ADHelper = Differentiation::AD::ResidualLinearization<
846 *   Differentiation::AD::NumberTypes::sacado_dfad,
847 *   double>;
848 *   using ADNumberType = typename ADHelper::ad_type;
849 *  
850 *   const FEValuesExtractors::Scalar u_fe(0);
851 *  
852 * @endcode
853 *
854 * With this, let us define the lambda function that will be used
855 * to compute the cell contributions to the Jacobian matrix and
856 * the right hand side:
857 *
858 * @code
859 *   const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
860 *   ScratchData &scratch_data,
861 *   CopyData &copy_data) {
862 *   const auto &fe_values = scratch_data.reinit(cell);
863 *   const unsigned int dofs_per_cell = fe_values.get_fe().n_dofs_per_cell();
864 *  
865 *   FullMatrix<double> &cell_matrix = copy_data.matrices[0];
866 *   Vector<double> &cell_rhs = copy_data.vectors[0];
867 *   std::vector<types::global_dof_index> &local_dof_indices =
868 *   copy_data.local_dof_indices[0];
869 *   cell->get_dof_indices(local_dof_indices);
870 *  
871 * @endcode
872 *
873 * We'll now create and initialize an instance of the AD helper class.
876 * number of local degrees of freedom that our solution vector has,
877 * i.e., the number @f$j@f$ in the per-element representation of the
878 * discretized solution vector
880 * that indicates how many solution coefficients are associated with
881 * each finite element. In deal.II, this equals
883 * the number of entries in the local residual vector that we will be
886 * method](https://en.wikipedia.org/wiki/Galerkin_method)) the number of
887 * local solution coefficients matches the number of local residual
888 * equations.
889 *
890 * @code
891 *   const unsigned int n_independent_variables = local_dof_indices.size();
892 *   const unsigned int n_dependent_variables = dofs_per_cell;
893 *   ADHelper ad_helper(n_independent_variables, n_dependent_variables);
894 *  
895 * @endcode
896 *
897 * Next we inform the helper of the values of the solution, i.e., the
899 * wish to linearize. As this is done on each element individually, we
900 * have to extract the solution coefficients from the global solution
901 * vector. In other words, we define all of those coefficients @f$U_j@f$
902 * where @f$j@f$ is a local degree of freedom as the independent variables
903 * that enter the computation of the vector @f$F(U)^{K}@f$ (the dependent
904 * function).
905 *
906
907 *
908 * Then we get the complete set of degree of freedom values as
911 * from this point until the object goes out of scope. So it is
913 * compute derivatives of the residual entries.
914 *
915 * @code
916 *   ad_helper.register_dof_values(current_solution, local_dof_indices);
917 *  
918 *   const std::vector<ADNumberType> &dof_values_ad =
919 *   ad_helper.get_sensitive_dof_values();
920 *  
921 * @endcode
922 *
924 * compute all values, (spatial) gradients, and the like based on
925 * "sensitive" AD degree of freedom values. In this instance we are
926 * retrieving the solution gradients at each quadrature point. Observe
927 * that the solution gradients are now sensitive
929 * as the scalar type and the @p dof_values_ad vector provides the local
930 * DoF values.
931 *
932 * @code
934 *   fe_values.n_quadrature_points);
935 *   fe_values[u_fe].get_function_gradients_from_local_dof_values(
937 *  
938 * @endcode
939 *
940 * The next variable that we declare will store the cell residual vector
942 * <b>very important</b> detail:
943 * Note that each entry in the vector is hand-initialized with a value
945 * libraries appear not to safely initialize the internal data
948 * of this program mentions this out of, generally bad, experience). So
949 * out of an abundance of caution it's worthwhile zeroing the initial
950 * value explicitly. After that, apart from a sign change the residual
951 * assembly looks much the same as we saw for the cell RHS vector before:
952 * We loop over all quadrature points, ensure that the coefficient now
953 * encodes its dependence on the (sensitive) finite element DoF values by
954 * using the correct `ADNumberType`, and finally we assemble the
955 * components of the residual vector. For complete clarity, the finite
956 * element shape functions (and their gradients, etc.) as well as the
957 * "JxW" values remain scalar
958 * valued, but the @p coeff and the @p old_solution_gradients at each
959 * quadrature point are computed in terms of the independent
960 * variables.
961 *
962 * @code
963 *   std::vector<ADNumberType> residual_ad(n_dependent_variables,
964 *   ADNumberType(0.0));
965 *   for (const unsigned int q : fe_values.quadrature_point_indices())
966 *   {
967 *   const ADNumberType coeff =
968 *   1.0 / std::sqrt(1.0 + old_solution_gradients[q] *
969 *   old_solution_gradients[q]);
970 *  
971 *   for (const unsigned int i : fe_values.dof_indices())
972 *   {
973 *   residual_ad[i] += (fe_values.shape_grad(i, q) // \nabla \phi_i
974 *   * coeff // * a_n
975 *   * old_solution_gradients[q]) // * \nabla u_n
976 *   * fe_values.JxW(q); // * dx
977 *   }
978 *   }
979 *  
980 * @endcode
981 *
982 * Once we have the full cell residual vector computed, we can register
983 * it with the helper class.
984 *
985
986 *
987 * Thereafter, we compute the residual values (basically,
988 * extracting the real values from what we already computed) and
989 * their Jacobian (the linearization of each residual component
990 * with respect to all cell DoFs) at the evaluation point. For
991 * the purposes of assembly into the global linear system, we
992 * have to respect the sign difference between the residual and
993 * the RHS contribution: For Newton's method, the right hand
994 * side vector needs to be equal to the *negative* residual
995 * vector.
996 *
997 * @code
998 *   ad_helper.register_residual_vector(residual_ad);
999 *  
1000 *   ad_helper.compute_residual(cell_rhs);
1001 *   cell_rhs *= -1.0;
1002 *  
1003 *   ad_helper.compute_linearization(cell_matrix);
1004 *   };
1005 *  
1006 * @endcode
1007 *
1008 * The remainder of the function equals what we had previously:
1009 *
1010 * @code
1011 *   const auto copier = [this](const CopyData &copy_data) {
1012 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1013 *   const Vector<double> &cell_rhs = copy_data.vectors[0];
1014 *   const std::vector<types::global_dof_index> &local_dof_indices =
1015 *   copy_data.local_dof_indices[0];
1016 *  
1017 *   zero_constraints.distribute_local_to_global(
1018 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1019 *   };
1020 *  
1021 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
1022 *   cell_worker,
1023 *   copier,
1024 *   sample_scratch_data,
1025 *   sample_copy_data,
1027 *   }
1028 *  
1029 * @endcode
1030 *
1031 *
1032 * <a name="step_72-Assemblyviadifferentiationoftheenergyfunctional"></a>
1034 *
1035
1036 *
1037 * In this third approach, we compute residual and Jacobian as first
1038 * and second derivatives of the local energy functional
1039 * @f[
1040 * E\left( U \right)^K
1041 * \dealcoloneq \int\limits_{K} \Psi \left( u \right) \, dV
1042 * \approx \sum\limits_{q}^{n_{\textrm{q-points}}} \Psi \left( u \left(
1043 * \mathbf{X}_{q} \right) \right) \underbrace{\vert J_{q} \vert \times
1044 * W_{q}}_{\text{JxW(q)}}
1045 * @f]
1046 * with the energy density given by
1047 * @f[
1048 * \Psi \left( u \right) = \sqrt{1+|\nabla u|^{2}} .
1049 * @f]
1050 *
1051
1052 *
1053 * Let us again see how this is done:
1054 *
1055 * @code
1056 *   template <int dim>
1058 *   {
1059 *   system_matrix = 0;
1060 *   system_rhs = 0;
1061 *  
1062 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1063 *  
1064 *   using ScratchData = MeshWorker::ScratchData<dim>;
1065 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
1066 *   using CellIteratorType = decltype(dof_handler.begin_active());
1067 *  
1068 *   const ScratchData sample_scratch_data(fe,
1073 *   const CopyData sample_copy_data(dofs_per_cell);
1074 *  
1075 * @endcode
1076 *
1078 * class that will automatically compute both the residual and its
1079 * linearization from the cell contribution to an energy functional using
1080 * nested Sacado forward automatic differentiation types.
1081 * The selected number types can be used to compute both first and
1082 * second derivatives. We require this, as the residual defined as the
1084 * its gradient). We'll then need to linearize the residual, implying that
1085 * second derivatives of the potential energy must be computed. You might
1086 * want to compare this with the definition of `ADHelper` used int
1087 * previous function, where we used
1088 * `Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,double>`.
1089 *
1090 * @code
1091 *   using ADHelper = Differentiation::AD::EnergyFunctional<
1092 *   Differentiation::AD::NumberTypes::sacado_dfad_dfad,
1093 *   double>;
1094 *   using ADNumberType = typename ADHelper::ad_type;
1095 *  
1096 *   const FEValuesExtractors::Scalar u_fe(0);
1097 *  
1098 * @endcode
1099 *
1100 * Let us then again define the lambda function that does the integration on
1101 * a cell.
1102 *
1103
1104 *
1105 * To initialize an instance of the helper class, we now only require
1106 * that the number of independent variables (that is, the number
1107 * of degrees of freedom associated with the element solution vector)
1108 * are known up front. This is because the second-derivative matrix that
1109 * results from an energy functional is necessarily square (and also,
1110 * incidentally, symmetric).
1111 *
1112 * @code
1113 *   const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
1114 *   ScratchData &scratch_data,
1115 *   CopyData &copy_data) {
1116 *   const auto &fe_values = scratch_data.reinit(cell);
1117 *  
1118 *   FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1119 *   Vector<double> &cell_rhs = copy_data.vectors[0];
1120 *   std::vector<types::global_dof_index> &local_dof_indices =
1121 *   copy_data.local_dof_indices[0];
1122 *   cell->get_dof_indices(local_dof_indices);
1123 *  
1124 *   const unsigned int n_independent_variables = local_dof_indices.size();
1125 *   ADHelper ad_helper(n_independent_variables);
1126 *  
1127 * @endcode
1128 *
1129 * Once more, we register all cell DoFs values with the helper, followed
1130 * by extracting the "sensitive" variant of these values that are to be
1131 * used in subsequent operations that must be differentiated -- one of
1132 * those being the calculation of the solution gradients.
1133 *
1134 * @code
1135 *   ad_helper.register_dof_values(current_solution, local_dof_indices);
1136 *  
1137 *   const std::vector<ADNumberType> &dof_values_ad =
1138 *   ad_helper.get_sensitive_dof_values();
1139 *  
1140 *   std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(
1141 *   fe_values.n_quadrature_points);
1142 *   fe_values[u_fe].get_function_gradients_from_local_dof_values(
1143 *   dof_values_ad, old_solution_gradients);
1144 *  
1145 * @endcode
1146 *
1147 * We next create a variable that stores the cell total energy.
1148 * Once more we emphasize that we explicitly zero-initialize this value,
1149 * thereby ensuring the integrity of the data for this starting value.
1150 *
1151
1152 *
1153 * The aim for our approach is then to compute the cell total
1154 * energy, which is the sum of the internal (due to right hand
1155 * side functions, typically linear in @f$U@f$) and external
1156 * energies. In this particular case, we have no external
1157 * energies (e.g., from source terms or Neumann boundary
1158 * conditions), so we'll focus on the internal energy part.
1159 *
1160
1161 *
1163 * the following lines:
1164 *
1165 * @code
1167 *   for (const unsigned int q : fe_values.quadrature_point_indices())
1168 *   {
1171 *  
1172 *   energy_ad += psi * fe_values.JxW(q);
1173 *   }
1174 *  
1175 * @endcode
1176 *
1177 * After we've computed the total energy on this cell, we'll
1178 * register it with the helper. Based on that, we may now
1179 * compute the desired quantities, namely the residual values
1181 * Newton right hand side needs to be the negative of the
1182 * residual:
1183 *
1184 * @code
1185 *   ad_helper.register_energy_functional(energy_ad);
1186 *  
1187 *   ad_helper.compute_residual(cell_rhs);
1188 *   cell_rhs *= -1.0;
1189 *  
1190 *   ad_helper.compute_linearization(cell_matrix);
1191 *   };
1192 *  
1193 * @endcode
1194 *
1195 * As in the previous two functions, the remainder of the function is as
1196 * before:
1197 *
1198 * @code
1199 *   const auto copier = [this](const CopyData &copy_data) {
1200 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1201 *   const Vector<double> &cell_rhs = copy_data.vectors[0];
1202 *   const std::vector<types::global_dof_index> &local_dof_indices =
1203 *   copy_data.local_dof_indices[0];
1204 *  
1205 *   zero_constraints.distribute_local_to_global(
1206 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1207 *   };
1208 *  
1209 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
1210 *   cell_worker,
1211 *   copier,
1212 *   sample_scratch_data,
1213 *   sample_copy_data,
1215 *   }
1216 *  
1217 *  
1218 * @endcode
1219 *
1220 *
1221 * <a name="step_72-MinimalSurfaceProblemsolve"></a>
1222 * <h4>MinimalSurfaceProblem::solve</h4>
1223 *
1224
1225 *
1226 * The solve function is the same as is used in @ref step_15 "step-15".
1227 *
1228 * @code
1229 *   template <int dim>
1231 *   {
1232 *   SolverControl solver_control(system_rhs.size(),
1233 *   system_rhs.l2_norm() * 1e-6);
1234 *   SolverCG<Vector<double>> solver(solver_control);
1235 *  
1236 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1237 *   preconditioner.initialize(system_matrix, 1.2);
1238 *  
1239 *   solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
1240 *  
1241 *   zero_constraints.distribute(newton_update);
1242 *  
1243 *   const double alpha = determine_step_length();
1245 *   }
1246 *  
1247 *  
1248 * @endcode
1249 *
1250 *
1251 * <a name="step_72-MinimalSurfaceProblemrefine_mesh"></a>
1253 *
1254
1255 *
1256 * Nothing has changed since @ref step_15 "step-15" with respect to the mesh refinement
1257 * procedure and transfer of the solution between adapted meshes.
1258 *
1259 * @code
1260 *   template <int dim>
1262 *   {
1264 *  
1266 *   dof_handler,
1267 *   QGauss<dim - 1>(fe.degree + 1),
1268 *   std::map<types::boundary_id, const Function<dim> *>(),
1271 *  
1274 *   0.3,
1275 *   0.03);
1276 *  
1277 *   triangulation.prepare_coarsening_and_refinement();
1278 *  
1281 *   solution_transfer.prepare_for_coarsening_and_refinement(coarse_solution);
1282 *  
1283 *   triangulation.execute_coarsening_and_refinement();
1284 *  
1285 *   setup_system();
1286 *  
1287 *   solution_transfer.interpolate(current_solution);
1289 *   }
1290 *  
1291 *  
1292 *  
1293 * @endcode
1294 *
1295 *
1296 * <a name="step_72-MinimalSurfaceProblemcompute_residual"></a>
1297 * <h4>MinimalSurfaceProblem::compute_residual</h4>
1298 *
1299
1300 *
1301 * ... as does the function used to compute the residual during the
1302 * solution iteration procedure. One could replace this by
1304 * but for simplicity we here simply copy what we already had in
1305 * @ref step_15 "step-15".
1306 *
1307 * @code
1308 *   template <int dim>
1309 *   double MinimalSurfaceProblem<dim>::compute_residual(const double alpha) const
1310 *   {
1311 *   Vector<double> residual(dof_handler.n_dofs());
1312 *  
1313 *   Vector<double> evaluation_point(dof_handler.n_dofs());
1316 *  
1317 *   FEValues<dim> fe_values(fe,
1321 *  
1322 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1323 *   const unsigned int n_q_points = quadrature_formula.size();
1324 *  
1325 *   Vector<double> cell_residual(dofs_per_cell);
1326 *   std::vector<Tensor<1, dim>> gradients(n_q_points);
1327 *  
1328 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1329 *  
1330 *   for (const auto &cell : dof_handler.active_cell_iterators())
1331 *   {
1332 *   cell_residual = 0;
1333 *   fe_values.reinit(cell);
1334 *  
1335 *   fe_values.get_function_gradients(evaluation_point, gradients);
1336 *  
1337 *   for (unsigned int q = 0; q < n_q_points; ++q)
1338 *   {
1339 *   const double coeff =
1340 *   1.0 / std::sqrt(1.0 + gradients[q] * gradients[q]);
1341 *  
1342 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1343 *   cell_residual(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
1344 *   * coeff // * a_n
1345 *   * gradients[q] // * \nabla u_n
1346 *   * fe_values.JxW(q)); // * dx
1347 *   }
1348 *  
1349 *   cell->get_dof_indices(local_dof_indices);
1350 *   zero_constraints.distribute_local_to_global(cell_residual,
1351 *   local_dof_indices,
1352 *   residual);
1353 *   }
1354 *  
1355 *   return residual.l2_norm();
1356 *   }
1357 *  
1358 *  
1359 *  
1360 * @endcode
1361 *
1362 *
1363 * <a name="step_72-MinimalSurfaceProblemdetermine_step_length"></a>
1365 *
1366
1367 *
1368 * The choice of step length (or, under-relaxation factor) for the nonlinear
1370 * @ref step_15 "step-15".
1371 *
1372 * @code
1373 *   template <int dim>
1375 *   {
1376 *   return 0.1;
1377 *   }
1378 *  
1379 *  
1380 *  
1381 * @endcode
1382 *
1383 *
1384 * <a name="step_72-MinimalSurfaceProblemoutput_results"></a>
1386 *
1387
1388 *
1389 * This last function to be called from `run()` outputs the current solution
1390 * (and the Newton update) in graphical form as a VTU file. It is entirely the
1391 * same as what has been used in previous tutorials.
1392 *
1393 * @code
1396 *   const unsigned int refinement_cycle) const
1397 *   {
1398 *   DataOut<dim> data_out;
1399 *  
1400 *   data_out.attach_dof_handler(dof_handler);
1401 *   data_out.add_data_vector(current_solution, "solution");
1402 *   data_out.add_data_vector(newton_update, "update");
1403 *   data_out.build_patches();
1404 *  
1405 *   const std::string filename =
1406 *   "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
1407 *   std::ofstream output(filename);
1408 *   data_out.write_vtu(output);
1409 *   }
1410 *  
1411 *  
1412 * @endcode
1413 *
1414 *
1415 * <a name="step_72-MinimalSurfaceProblemrun"></a>
1416 * <h4>MinimalSurfaceProblem::run</h4>
1417 *
1418
1419 *
1421 * in @ref step_15 "step-15". The only observable changes are that we can now choose (via
1422 * the parameter file) what the final acceptable tolerance for the system
1423 * residual is, and that we can choose which method of assembly we wish to
1425 * message which indicates the selection. Since we're interested in comparing
1426 * the time taken to assemble for each of the three methods, we've also
1427 * added a timer that keeps a track of how much time is spent during assembly.
1428 * We also track the time taken to solve the linear system, so that we can
1430 * the longest time to execute.
1431 *
1432 * @code
1433 *   template <int dim>
1435 *   const double tolerance)
1436 *   {
1437 *   std::cout << "******** Assembly approach ********" << std::endl;
1438 *   const std::array<std::string, 3> method_descriptions = {
1439 *   {"Unassisted implementation (full hand linearization).",
1440 *   "Automated linearization of the finite element residual.",
1441 *   "Automated computation of finite element residual and linearization using a variational formulation."}};
1443 *   std::cout << method_descriptions[formulation] << std::endl << std::endl;
1444 *  
1445 *  
1447 *  
1449 *   triangulation.refine_global(2);
1450 *  
1451 *   setup_system();
1453 *  
1454 *   double last_residual_norm = std::numeric_limits<double>::max();
1455 *   unsigned int refinement_cycle = 0;
1456 *   do
1457 *   {
1458 *   std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
1459 *  
1460 *   if (refinement_cycle != 0)
1461 *   refine_mesh();
1462 *  
1463 *   std::cout << " Initial residual: " << compute_residual(0) << std::endl;
1464 *  
1465 *   for (unsigned int inner_iteration = 0; inner_iteration < 5;
1466 *   ++inner_iteration)
1467 *   {
1468 *   {
1469 *   TimerOutput::Scope t(timer, "Assemble");
1470 *  
1471 *   if (formulation == 0)
1473 *   else if (formulation == 1)
1475 *   else if (formulation == 2)
1477 *   else
1478 *   AssertThrow(false, ExcNotImplemented());
1479 *   }
1480 *  
1481 *   last_residual_norm = system_rhs.l2_norm();
1482 *  
1483 *   {
1484 *   TimerOutput::Scope t(timer, "Solve");
1485 *   solve();
1486 *   }
1487 *  
1488 *  
1489 *   std::cout << " Residual: " << compute_residual(0) << std::endl;
1490 *   }
1491 *  
1493 *  
1494 *   ++refinement_cycle;
1495 *   std::cout << std::endl;
1496 *   }
1497 *   while (last_residual_norm > tolerance);
1498 *   }
1499 *   } // namespace Step72
1500 *  
1501 * @endcode
1502 *
1503 *
1504 * <a name="step_72-Themainfunction"></a>
1505 * <h4>The main function</h4>
1506 *
1507
1508 *
1509 * Finally the main function. This follows the scheme of most other main
1511 * - We call Utilities::MPI::MPI_InitFinalize in order to set up (via a hidden
1512 * default parameter) the number of threads using the execution of
1513 * multithreaded tasks.
1516 * program.
1517 *
1518 * @code
1519 *   int main(int argc, char *argv[])
1520 *   {
1521 *   try
1522 *   {
1523 *   using namespace Step72;
1524 *  
1526 *  
1527 *   std::string prm_file;
1528 *   if (argc > 1)
1529 *   prm_file = argv[1];
1530 *   else
1531 *   prm_file = "parameters.prm";
1532 *  
1533 *   const MinimalSurfaceProblemParameters parameters;
1535 *  
1537 *   minimal_surface_problem_2d.run(parameters.formulation,
1538 *   parameters.tolerance);
1539 *   }
1540 *   catch (std::exception &exc)
1541 *   {
1542 *   std::cerr << std::endl
1543 *   << std::endl
1544 *   << "----------------------------------------------------"
1545 *   << std::endl;
1546 *   std::cerr << "Exception on processing: " << std::endl
1547 *   << exc.what() << std::endl
1548 *   << "Aborting!" << std::endl
1549 *   << "----------------------------------------------------"
1550 *   << std::endl;
1551 *  
1552 *   return 1;
1553 *   }
1554 *   catch (...)
1555 *   {
1556 *   std::cerr << std::endl
1557 *   << std::endl
1558 *   << "----------------------------------------------------"
1559 *   << std::endl;
1560 *   std::cerr << "Unknown exception!" << std::endl
1561 *   << "Aborting!" << std::endl
1562 *   << "----------------------------------------------------"
1563 *   << std::endl;
1564 *   return 1;
1565 *   }
1566 *   return 0;
1567 *   }
1568 * @endcode
1569<a name="step_72-Results"></a><h1>Results</h1>
1570
1571
1573in @ref step_15 "step-15", there is nothing to report about that. The only outwardly noticeable
1574difference between them is that, by default, this program will only run 9 mesh
1575refinement steps (as opposed to @ref step_15 "step-15", which executes 11 refinements).
1577header text that prints which assembly method is being used, and the final
1579
1580@code
1581Mesh refinement step 0
1582 Initial residual: 1.53143
1583 Residual: 1.08746
1584 Residual: 0.966748
1585 Residual: 0.859602
1586 Residual: 0.766462
1587 Residual: 0.685475
1588
1589...
1590
1591Mesh refinement step 9
1592 Initial residual: 0.00924594
1593 Residual: 0.00831928
1594 Residual: 0.0074859
1595 Residual: 0.0067363
1596 Residual: 0.00606197
1597 Residual: 0.00545529
1598@endcode
1599
1600So what is interesting for us to compare is how long the assembly process takes
1602Below is the output for the hand linearization (as computed on a circa 2012
1603four core, eight thread laptop -- but we're only really interested in the
1604relative time between the different implementations):
1605@code
1606******** Assembly approach ********
1607Unassisted implementation (full hand linearization).
1608
1609...
1610
1611+---------------------------------------------+------------+------------+
1612| Total wallclock time elapsed since start | 35.1s | |
1613| | | |
1614| Section | no. calls | wall time | % of total |
1615+---------------------------------+-----------+------------+------------+
1616| Assemble | 50 | 1.56s | 4.5% |
1617| Solve | 50 | 30.8s | 88% |
1618+---------------------------------+-----------+------------+------------+
1619@endcode
1620And for the implementation that linearizes the residual in an automated
1621manner using the Sacado dynamic forward AD number type:
1622@code
1623******** Assembly approach ********
1624Automated linearization of the finite element residual.
1625
1626...
1627
1628+---------------------------------------------+------------+------------+
1629| Total wallclock time elapsed since start | 40.1s | |
1630| | | |
1631| Section | no. calls | wall time | % of total |
1632+---------------------------------+-----------+------------+------------+
1633| Assemble | 50 | 8.8s | 22% |
1634| Solve | 50 | 28.6s | 71% |
1635+---------------------------------+-----------+------------+------------+
1636@endcode
1637And, lastly, for the implementation that computes both the residual and
1638its linearization directly from an energy functional (using nested Sacado
1639dynamic forward AD numbers):
1640@code
1641******** Assembly approach ********
1642Automated computation of finite element residual and linearization using a variational formulation.
1643
1644...
1645
1646+---------------------------------------------+------------+------------+
1647| Total wallclock time elapsed since start | 48.8s | |
1648| | | |
1649| Section | no. calls | wall time | % of total |
1650+---------------------------------+-----------+------------+------------+
1651| Assemble | 50 | 16.7s | 34% |
1652| Solve | 50 | 29.3s | 60% |
1653+---------------------------------+-----------+------------+------------+
1654@endcode
1655
1658over all refinement steps, using one level of automatic differentiation resulted
1662Unsurprisingly, the overall time spent solving the linear system remained unchanged.
1665at the finite element level. For many, this might mean that leveraging higher-order
1669of the finite element residual, which offers a lot of convenience at a measurable,
1671not re-building the Newton matrix in every step -- a topic that is
1672explored in substantial depth in @ref step_77 "step-77".
1673
1675(e.g., how many components there are in the solution, what the nature of the function
1679
1680
1681<a name="step_72-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1682
1683
1684Like @ref step_71 "step-71", there are a few items related to automatic differentiation that could
1690 employed at the finite-element level, mixed differentiation modes ("RAD-FAD")
1692 mode ("FAD-FAD") types employed here. The reason that the RAD-FAD type was not
1693 selected by default is that, at the time of writing, there remain some
1695 This is documented in the @ref auto_symb_diff topic.
1696- It might be the case that using reduced precision types (i.e., `float`) as the
1698 expense during assembly. Using `float` as the data type for the
1699 matrix and the residual is not unreasonable, given that the Newton
1700 update is only meant to get us closer to the solution, but not
1701 actually *to* the solution; as a consequence, it makes sense to
1702 consider using reduced-precision data types for computing these
1703 updates, and then accumulating these updates in a solution vector
1704 that uses the full `double` precision accuracy.
1707 approach adopted in @ref step_71 "step-71", and pushes the starting point for the automatic
1711 each cell quadrature point).
1712- @ref step_77 "step-77" is yet another variation of @ref step_15 "step-15" that addresses a very
1714 necessary to re-build the Newton matrix in every nonlinear
1715 iteration. Given that the results above show that using automatic
1716 differentiation comes at a cost, the techniques in @ref step_77 "step-77" have the
1717 potential to offset these costs to some degree. It would therefore
1719 techniques in @ref step_77 "step-77". For production codes, this would certainly be
1720 the way to go.
1721 *
1722 *
1723<a name="step_72-PlainProg"></a>
1724<h1> The plain program</h1>
1725@include "step-72.cc"
1726*/
const unsigned int dofs_per_cell
Definition fe_data.h:436
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, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, 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)
Tensor< rank, dim, Number > sum(const Tensor< rank, dim, Number > &local, const MPI_Comm mpi_communicator)
constexpr void clear()
friend class Tensor
Definition tensor.h:865
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
Definition tensor.h:851
@ wall_times
Definition timer.h:651
Point< 2 > second
Definition grid_out.cc:4630
Point< 2 > first
Definition grid_out.cc:4629
unsigned int level
Definition grid_out.cc:4632
static ::ExceptionBase & ExcNotImplemented()
#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:281
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={}, 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.
std::vector< index_type > data
Definition mpi.cc:740
std::size_t size
Definition mpi.cc:739
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > &center={}, 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.
constexpr char U
constexpr types::blas_int zero
constexpr types::blas_int one
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:74
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:130
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:193
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)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:474
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={})
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)
constexpr double PI
Definition numbers.h:241
STL namespace.
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition types.h:32
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation