450 * <a name=
"ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
451 * <h3>The <code>MinimalSurfaceProblem</code>
class implementation</h3>
456 * <a name=
"Constructorandsetupfunctions"></a>
461 * The following few
functions are also essentially copies of what
462 * @ref step_15
"step-15" already does, and so there is little to discuss.
466 * MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
475 *
void MinimalSurfaceProblem<dim>::setup_system(
const bool initial_step)
481 * dof_handler.distribute_dofs(fe);
482 * current_solution.reinit(dof_handler.n_dofs());
484 * hanging_node_constraints.clear();
486 * hanging_node_constraints);
487 * hanging_node_constraints.close();
493 * hanging_node_constraints.condense(dsp);
495 * sparsity_pattern.copy_from(dsp);
496 * jacobian_matrix.reinit(sparsity_pattern);
497 * jacobian_matrix_factorization.reset();
505 * <a name=
"AssemblingandfactorizingtheJacobianmatrix"></a>
506 * <h4>Assembling and factorizing the Jacobian
matrix</h4>
510 * The following function is then responsible
for assembling and factorizing
511 * the Jacobian
matrix. The
first half of the function is in essence the
512 * `assemble_system()` function of @ref step_15
"step-15", except that it does not deal with
513 * also forming a right hand side vector (i.e., the residual) since we
do not
514 *
always have to
do these operations at the same time.
518 * We put the whole assembly functionality into a code block enclosed by curly
520 * time is spent in
this code block, excluding everything that happens in
this
521 * function after the matching closing brace `}`.
525 *
void MinimalSurfaceProblem<dim>::compute_and_factorize_jacobian(
531 * std::cout <<
" Computing Jacobian matrix" << std::endl;
533 *
const QGauss<dim> quadrature_formula(fe.degree + 1);
535 * jacobian_matrix = 0;
538 * quadrature_formula,
542 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
543 *
const unsigned int n_q_points = quadrature_formula.size();
547 * std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
549 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
551 *
for (
const auto &cell : dof_handler.active_cell_iterators())
555 * fe_values.reinit(cell);
557 * fe_values.get_function_gradients(evaluation_point,
558 * evaluation_point_gradients);
560 *
for (
unsigned int q = 0; q < n_q_points; ++q)
562 *
const double coeff =
563 * 1.0 /
std::sqrt(1 + evaluation_point_gradients[q] *
564 * evaluation_point_gradients[q]);
566 *
for (
unsigned int i = 0; i < dofs_per_cell; ++i)
568 *
for (
unsigned int j = 0; j < dofs_per_cell; ++j)
570 * (((fe_values.shape_grad(i, q)
572 * * fe_values.shape_grad(j, q))
574 * (fe_values.shape_grad(i, q)
575 * * coeff * coeff * coeff
577 * (fe_values.shape_grad(j, q)
578 * * evaluation_point_gradients[q])
579 * * evaluation_point_gradients[q]))
580 * * fe_values.JxW(q));
584 * cell->get_dof_indices(local_dof_indices);
585 * hanging_node_constraints.distribute_local_to_global(cell_matrix,
590 * std::map<types::global_dof_index, double> boundary_values;
605 * The
second half of the function then deals with factorizing the
607 *
object and by assigning it to the member variable
608 * `jacobian_matrix_factorization`, we also destroy whatever
object that
609 * pointer previously pointed to (
if any). Then we tell the
object to
610 * factorize the Jacobian.
614 * As above, we enclose
this block of code into curly braces and use a timer
615 * to assess how
long this part of the program takes.
619 * (Strictly speaking, we don
't actually need the matrix any more after we
620 * are done here, and could throw the matrix object away. A code intended to
621 * be memory efficient would do this, and only create the matrix object in
622 * this function, rather than as a member variable of the surrounding class.
623 * We omit this step here because using the same coding style as in previous
624 * tutorial programs breeds familiarity with the common style and helps make
625 * these tutorial programs easier to read.)
629 * TimerOutput::Scope t(computing_timer, "factorizing the Jacobian");
631 * std::cout << " Factorizing Jacobian matrix" << std::endl;
633 * jacobian_matrix_factorization = std::make_unique<SparseDirectUMFPACK>();
634 * jacobian_matrix_factorization->factorize(jacobian_matrix);
643 * <a name="Computingtheresidualvector"></a>
644 * <h4>Computing the residual vector</h4>
648 * The second part of what `assemble_system()` used to do in @ref step_15 "step-15" is
649 * computing the residual vector, i.e., the right hand side vector of the
650 * Newton linear systems. We have broken this out of the previous function,
651 * but the following function will be easy to understand if you understood
652 * what `assemble_system()` in @ref step_15 "step-15" did. Importantly, however, we need to
653 * compute the residual not linearized around the current solution vector, but
654 * whatever we get from KINSOL. This is necessary for operations such as line
655 * search where we want to know what the residual @f$F(U^k + \alpha_k \delta
656 * U^K)@f$ is for different values of @f$\alpha_k@f$; KINSOL in those cases simply
657 * gives us the argument to the function @f$F@f$ and we then compute the residual
658 * @f$F(\cdot)@f$ at this point.
662 * The function prints the norm of the so-computed residual at the end as a
663 * way for us to follow along the progress of the program.
667 * void MinimalSurfaceProblem<dim>::compute_residual(
668 * const Vector<double> &evaluation_point,
669 * Vector<double> & residual)
671 * TimerOutput::Scope t(computing_timer, "assembling the residual");
673 * std::cout << " Computing residual vector..." << std::flush;
675 * const QGauss<dim> quadrature_formula(fe.degree + 1);
676 * FEValues<dim> fe_values(fe,
677 * quadrature_formula,
678 * update_gradients | update_quadrature_points |
679 * update_JxW_values);
681 * const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
682 * const unsigned int n_q_points = quadrature_formula.size();
684 * Vector<double> cell_residual(dofs_per_cell);
685 * std::vector<Tensor<1, dim>> evaluation_point_gradients(n_q_points);
687 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
689 * for (const auto &cell : dof_handler.active_cell_iterators())
692 * fe_values.reinit(cell);
694 * fe_values.get_function_gradients(evaluation_point,
695 * evaluation_point_gradients);
698 * for (unsigned int q = 0; q < n_q_points; ++q)
700 * const double coeff =
701 * 1.0 / std::sqrt(1 + evaluation_point_gradients[q] *
702 * evaluation_point_gradients[q]);
704 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
706 * (fe_values.shape_grad(i, q) // \nabla \phi_i
708 * * evaluation_point_gradients[q] // * \nabla u_n
709 * * fe_values.JxW(q)); // * dx
712 * cell->get_dof_indices(local_dof_indices);
713 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
714 * residual(local_dof_indices[i]) += cell_residual(i);
717 * hanging_node_constraints.condense(residual);
719 * for (const types::global_dof_index i :
720 * DoFTools::extract_boundary_dofs(dof_handler))
723 * for (const types::global_dof_index i :
724 * DoFTools::extract_hanging_node_dofs(dof_handler))
727 * std::cout << " norm=" << residual.l2_norm() << std::endl;
735 * <a name="SolvinglinearsystemswiththeJacobianmatrix"></a>
736 * <h4>Solving linear systems with the Jacobian matrix</h4>
740 * Next up is the function that implements the solution of a linear system
741 * with the Jacobian matrix. Since we have already factored the matrix when we
742 * built the matrix, solving a linear system comes down to applying the
743 * inverse matrix to the given right hand side vector: This is what the
744 * SparseDirectUMFPACK::vmult() function does that we use here. Following
745 * this, we have to make sure that we also address the values of hanging nodes
746 * in the solution vector, and this is done using
747 * AffineConstraints::distribute().
751 * The function takes an additional, but unused, argument `tolerance` that
752 * indicates how accurately we have to solve the linear system. The meaning of
753 * this argument is discussed in the introduction in the context of the
754 * "Eisenstat Walker trick", but since we are using a direct rather than an
755 * iterative solver, we are not using this opportunity to solve linear systems
760 * void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
761 * Vector<double> & solution,
762 * const double /*tolerance*/)
764 * TimerOutput::Scope t(computing_timer, "linear system solve");
766 * std::cout << " Solving linear system" << std::endl;
768 * jacobian_matrix_factorization->vmult(solution, rhs);
770 * hanging_node_constraints.distribute(solution);
778 * <a name="Refiningthemeshsettingboundaryvaluesandgeneratinggraphicaloutput"></a>
779 * <h4>Refining the mesh, setting boundary values, and generating graphical output</h4>
783 * The following three functions are again simply copies of the ones in
784 * @ref step_15 "step-15":
788 * void MinimalSurfaceProblem<dim>::refine_mesh()
790 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
792 * KellyErrorEstimator<dim>::estimate(
794 * QGauss<dim - 1>(fe.degree + 1),
795 * std::map<types::boundary_id, const Function<dim> *>(),
797 * estimated_error_per_cell);
799 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
800 * estimated_error_per_cell,
804 * triangulation.prepare_coarsening_and_refinement();
806 * SolutionTransfer<dim> solution_transfer(dof_handler);
807 * solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
809 * triangulation.execute_coarsening_and_refinement();
811 * dof_handler.distribute_dofs(fe);
813 * Vector<double> tmp(dof_handler.n_dofs());
814 * solution_transfer.interpolate(current_solution, tmp);
815 * current_solution = std::move(tmp);
817 * hanging_node_constraints.clear();
819 * DoFTools::make_hanging_node_constraints(dof_handler,
820 * hanging_node_constraints);
821 * hanging_node_constraints.close();
823 * hanging_node_constraints.distribute(current_solution);
825 * set_boundary_values();
827 * setup_system(/*initial_step=*/false);
833 * void MinimalSurfaceProblem<dim>::set_boundary_values()
835 * std::map<types::global_dof_index, double> boundary_values;
836 * VectorTools::interpolate_boundary_values(dof_handler,
838 * BoundaryValues<dim>(),
840 * for (const auto &boundary_value : boundary_values)
841 * current_solution(boundary_value.first) = boundary_value.second;
843 * hanging_node_constraints.distribute(current_solution);
849 * void MinimalSurfaceProblem<dim>::output_results(
850 * const unsigned int refinement_cycle)
852 * TimerOutput::Scope t(computing_timer, "graphical output");
854 * DataOut<dim> data_out;
856 * data_out.attach_dof_handler(dof_handler);
857 * data_out.add_data_vector(current_solution, "solution");
858 * data_out.build_patches();
860 * const std::string filename =
861 * "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
862 * std::ofstream output(filename);
863 * data_out.write_vtu(output);
871 * <a name="Therunfunctionandtheoveralllogicoftheprogram"></a>
872 * <h4>The run() function and the overall logic of the program</h4>
876 * The only function that *really* is interesting in this program is the one
877 * that drives the overall algorithm of starting on a coarse mesh, doing some
878 * mesh refinement cycles, and on each mesh using KINSOL to find the solution
879 * of the nonlinear algebraic equation we obtain from discretization on this
880 * mesh. The `refine_mesh()` function above makes sure that the solution on
881 * one mesh is used as the starting guess on the next mesh. We also use a
882 * TimerOutput object to measure how much time every operation on each mesh
883 * costs, and reset the timer at the beginning of each cycle.
887 * As discussed in the introduction, it is not necessary to solve problems on
888 * coarse meshes particularly accurately since these will only solve as
889 * starting guesses for the next mesh. As a consequence, we will use a target
891 * @f$\tau=10^{-3} \frac{1}{10^k}@f$ for the @f$k@f$th mesh refinement cycle.
895 * All of this is encoded in the first part of this function:
899 * void MinimalSurfaceProblem<dim>::run()
901 * GridGenerator::hyper_ball(triangulation);
902 * triangulation.refine_global(2);
904 * setup_system(/*initial_step=*/true);
905 * set_boundary_values();
907 * for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
908 * ++refinement_cycle)
910 * computing_timer.reset();
911 * std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
913 * if (refinement_cycle != 0)
916 * const double target_tolerance = 1e-3 * std::pow(0.1, refinement_cycle);
917 * std::cout << " Target_tolerance: " << target_tolerance << std::endl
922 * This is where the fun starts. At the top we create the KINSOL solver
923 * object and feed it with an object that encodes a number of additional
924 * specifics (of which we only change the nonlinear tolerance we want to
925 * reach; but you might want to look into what other members the
926 * SUNDIALS::KINSOL::AdditionalData class has and play with them).
930 * typename SUNDIALS::KINSOL<Vector<double>>::AdditionalData
932 * additional_data.function_tolerance = target_tolerance;
934 * SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
938 * Then we have to describe the operations that were already mentioned
939 * in the introduction. In essence, we have to teach KINSOL how to (i)
940 * resize a vector to the correct size, (ii) compute the residual
941 * vector, (iii) compute the Jacobian matrix (during which we also
942 * compute its factorization), and (iv) solve a linear system with the
947 * All four of these operations are represented by member variables of
948 * the SUNDIALS::KINSOL class that are of type `std::function`, i.e.,
949 * they are objects to which we can assign a pointer to a function or,
950 * as we do here, a "lambda function" that takes the appropriate
951 * arguments and returns the appropriate information. By convention,
952 * KINSOL wants that functions doing something nontrivial return an
953 * integer where zero indicates success. It turns out that we can do
954 * all of this in just 25 lines of code.
958 * (If you're not familiar what
"lambda functions" are, take
959 * a look at @ref step_12
"step-12" or at the
960 * [wikipedia page](https:
962 * wants to define a function with a certain
set of
963 * arguments, but (i) not make it a named
functions because,
964 * typically, the function is used in only one place and it
965 * seems unnecessary to give it a global name; and (ii) that
966 * the function has access to some of the variables that
967 * exist at the place where it is defined, including member
969 * ultimately quite useful.)
973 * At the very
end of the code block we then tell KINSOL to go to work
974 * and solve our problem. The member
functions called from the
975 *
'residual',
'setup_jacobian', and
'solve_jacobian_system' functions
976 * will then print output to screen that allows us to follow along
977 * with the progress of the program.
981 * x.reinit(dof_handler.n_dofs());
984 * nonlinear_solver.residual =
987 * compute_residual(evaluation_point, residual);
992 * nonlinear_solver.setup_jacobian =
995 * compute_and_factorize_jacobian(current_u);
1000 * nonlinear_solver.solve_with_jacobian = [&](
const Vector<double> &rhs,
1002 *
const double tolerance) {
1003 * this->solve(rhs, dst, tolerance);
1008 * nonlinear_solver.solve(current_solution);
1013 * The rest is then just house-keeping: Writing data to a file
for
1014 * visualizing, and showing a summary of the timing collected so that we
1015 * can interpret how
long each operation has taken, how often it was
1019 * output_results(refinement_cycle);
1021 * computing_timer.print_summary();
1023 * std::cout << std::endl;
1033 *
using namespace Step77;
1035 * MinimalSurfaceProblem<2> laplace_problem_2d;
1036 * laplace_problem_2d.run();
1038 *
catch (std::exception &exc)
1040 * std::cerr << std::endl
1042 * <<
"----------------------------------------------------"
1044 * std::cerr <<
"Exception on processing: " << std::endl
1045 * << exc.what() << std::endl
1046 * <<
"Aborting!" << std::endl
1047 * <<
"----------------------------------------------------"
1054 * std::cerr << std::endl
1056 * <<
"----------------------------------------------------"
1058 * std::cerr <<
"Unknown exception!" << std::endl
1059 * <<
"Aborting!" << std::endl
1060 * <<
"----------------------------------------------------"
1067<a name=
"Results"></a><h1>Results</h1>
1070When running the program, you get output that looks like
this:
1072Mesh refinement step 0
1073 Target_tolerance: 0.001
1075 Computing residual vector...
norm=0.231202
1076 Computing Jacobian
matrix
1077 Factorizing Jacobian
matrix
1078 Solving linear system
1079 Computing residual vector...
norm=0.231202
1080 Computing residual vector...
norm=0.171585
1081 Solving linear system
1082 Computing residual vector...
norm=0.171585
1083 Computing residual vector...
norm=0.127245
1084 Computing residual vector...
norm=0.0796471
1085 Solving linear system
1086 Computing residual vector...
norm=0.0796471
1087 Computing residual vector...
norm=0.0625301
1088 Solving linear system
1089 Computing residual vector...
norm=0.0625301
1090 Computing residual vector...
norm=0.0498864
1091 Solving linear system
1092 Computing residual vector...
norm=0.0498864
1093 Computing residual vector...
norm=0.0407765
1094 Solving linear system
1095 Computing residual vector...
norm=0.0407765
1096 Computing residual vector...
norm=0.0341589
1097 Solving linear system
1098 Computing residual vector...
norm=0.0341589
1099 Computing residual vector...
norm=0.0292867
1100 Solving linear system
1101 Computing residual vector...
norm=0.0292867
1102 Computing residual vector...
norm=0.0256309
1103 Computing residual vector...
norm=0.0223448
1104 Solving linear system
1105 Computing residual vector...
norm=0.0223448
1106 Computing residual vector...
norm=0.0202797
1107 Computing residual vector...
norm=0.0183817
1108 Solving linear system
1109 Computing residual vector...
norm=0.0183817
1110 Computing residual vector...
norm=0.0170464
1111 Computing residual vector...
norm=0.0157967
1112 Computing Jacobian
matrix
1113 Factorizing Jacobian
matrix
1114 Solving linear system
1115 Computing residual vector...
norm=0.0157967
1116 Computing residual vector...
norm=0.0141572
1117 Computing residual vector...
norm=0.012657
1118 Solving linear system
1119 Computing residual vector...
norm=0.012657
1120 Computing residual vector...
norm=0.0116863
1121 Computing residual vector...
norm=0.0107696
1122 Solving linear system
1123 Computing residual vector...
norm=0.0107696
1124 Computing residual vector...
norm=0.0100986
1125 Computing residual vector...
norm=0.00944829
1126 Computing residual vector...
norm=0.00822576
1127 Solving linear system
1128 Computing residual vector...
norm=0.00822576
1129 Computing residual vector...
norm=0.00781983
1130 Computing residual vector...
norm=0.00741619
1131 Computing residual vector...
norm=0.00661792
1132 Solving linear system
1133 Computing residual vector...
norm=0.00661792
1134 Computing residual vector...
norm=0.00630571
1135 Computing residual vector...
norm=0.00599457
1136 Computing residual vector...
norm=0.00537663
1137 Solving linear system
1138 Computing residual vector...
norm=0.00537663
1139 Computing residual vector...
norm=0.00512813
1140 Computing residual vector...
norm=0.00488033
1141 Computing residual vector...
norm=0.00438751
1142 Computing residual vector...
norm=0.00342052
1143 Solving linear system
1144 Computing residual vector...
norm=0.00342052
1145 Computing residual vector...
norm=0.00326581
1146 Computing residual vector...
norm=0.00311176
1147 Computing residual vector...
norm=0.00280617
1148 Computing residual vector...
norm=0.00220992
1149 Solving linear system
1150 Computing residual vector...
norm=0.00220992
1151 Computing residual vector...
norm=0.00209976
1152 Computing residual vector...
norm=0.00199943
1153 Solving linear system
1154 Computing residual vector...
norm=0.00199942
1155 Computing residual vector...
norm=0.00190953
1156 Computing residual vector...
norm=0.00182005
1157 Computing residual vector...
norm=0.00164259
1158 Computing residual vector...
norm=0.00129652
1161+---------------------------------------------+------------+------------+
1162| Total wallclock time elapsed since start | 0.192s | |
1164| Section | no. calls | wall time | % of total |
1165+---------------------------------+-----------+------------+------------+
1166| assembling the Jacobian | 2 | 0.0141s | 7.4% |
1167| assembling the residual | 61 | 0.168s | 88% |
1168| factorizing the Jacobian | 2 | 0.0016s | 0.83% |
1169| graphical output | 1 | 0.00385s | 2% |
1170| linear system solve | 19 | 0.0013s | 0.68% |
1171+---------------------------------+-----------+------------+------------+
1174Mesh refinement step 1
1175 Target_tolerance: 0.0001
1177 Computing residual vector...
norm=0.0883422
1178 Computing Jacobian
matrix
1179 Factorizing Jacobian
matrix
1180 Solving linear system
1181 Computing residual vector...
norm=0.0883422
1182 Computing residual vector...
norm=0.0607066
1183 Solving linear system
1184 Computing residual vector...
norm=0.0607066
1185 Computing residual vector...
norm=0.0437266
1186 Solving linear system
1187 Computing residual vector...
norm=0.0437266
1188 Computing residual vector...
norm=0.0327999
1189 Solving linear system
1190 Computing residual vector...
norm=0.0327999
1191 Computing residual vector...
norm=0.0255418
1192 Solving linear system
1193 Computing residual vector...
norm=0.0255417
1194 Computing residual vector...
norm=0.0206042
1195 Solving linear system
1196 Computing residual vector...
norm=0.0206042
1197 Computing residual vector...
norm=0.0171602
1198 Solving linear system
1199 Computing residual vector...
norm=0.0171602
1200 Computing residual vector...
norm=0.014689
1201 Solving linear system
1206The way
this should be interpreted is most easily explained by looking at
1207the
first few lines of the output on the
first mesh:
1209Mesh refinement step 0
1210Mesh refinement step 0
1211 Target_tolerance: 0.001
1213 Computing residual vector...
norm=0.231202
1214 Computing Jacobian
matrix
1215 Factorizing Jacobian
matrix
1216 Solving linear system
1217 Computing residual vector...
norm=0.231202
1218 Computing residual vector...
norm=0.171585
1219 Solving linear system
1220 Computing residual vector...
norm=0.171585
1221 Computing residual vector...
norm=0.127245
1222 Computing residual vector...
norm=0.0796471
1223 Solving linear system
1224 Computing residual vector...
norm=0.0796471
1227What is happening is
this:
1228- In the
first residual computation, KINSOL computes the residual to see whether
1229 the desired tolerance has been reached. The answer is no, so it requests the
1230 user program to compute the Jacobian
matrix (and the function then also
1232- KINSOL then instructs us to solve a linear system of the form
1233 @f$J_k \, \delta U_k = -F_k@f$ with
this matrix and the previously computed
1235- It is then time to determine how far we want to go in
this direction,
1236 i.e.,
do line search. To
this end, KINSOL
requires us to compute the
1237 residual vector @f$F(U_k + \alpha_k \delta U_k)@f$
for different step lengths
1238 @f$\alpha_k@f$. For the
first step above, it finds an acceptable @f$\alpha_k@f$
1239 after two tries, the
second time around it takes three tries.
1240- Having found a suitable updated solution @f$U_{k+1}@f$, the process is
1241 repeated except now KINSOL is happy with the current Jacobian
matrix
1242 and does not instruct us to re-build the
matrix and its factorization,
1243 and instead asks us to solve a linear system with that same
matrix.
1245The program also writes the solution to a VTU file at the
end
1246of each mesh refinement cycle, and it looks as follows:
1247<table width=
"60%" align=
"center">
1250 <img src=
"https://www.dealii.org/images/steps/developer/step-77.solution.png" alt=
"">
1256The key takeaway messages of
this program are the following:
1258- The solution is the same as the one we computed in @ref step_15
"step-15", i.e., the
1259 interfaces to %
SUNDIALS' KINSOL package really did what they were supposed
1260 to do. This should not come as a surprise, but the important point is that
1261 we don't have to spend the time implementing the complex algorithms that
1262 underlie advanced nonlinear solvers ourselves.
1264- KINSOL is able to avoid all sorts of operations such as rebuilding the
1265 Jacobian
matrix when that is not actually necessary. Comparing the
1266 number of linear solves in the output above with the number of times
1267 we rebuild the Jacobian and compute its factorization should make it
1268 clear that
this leads to very substantial savings in terms of compute
1269 times, without us having to implement the intricacies of algorithms
1270 that determine when we need to rebuild
this information.
1272<a name=
"extensions"></a>
1273<a name=
"Possibilitiesforextensions"></a><h3> Possibilities
for extensions </h3>
1276For all but the small problems we consider here, a sparse direct solver
1277requires too much time and memory -- we need an iterative solver like
1278we use in many other programs. The trade-off between constructing an
1279expensive preconditioner (say, a geometric or algebraic multigrid method)
1280is different in the current
case, however: Since we can re-use the same
1281matrix for numerous linear solves, we can
do the same
for the preconditioner
1282and putting more work into building a good preconditioner can more easily
1283be justified than
if we used it only
for a single linear solve as one
1284does
for many other situations.
1286But iterative solvers also afford other opportunities. For example (and as
1287discussed briefly in the introduction), we may not need to solve to
1288very high accuracy (small tolerances) in early nonlinear iterations as
long
1289as we are still far away from the actual solution. This was the basis of the
1290Eisenstat-Walker trick mentioned there.
1292KINSOL provides the function that does the linear solution with a target
1293tolerance that needs to be reached. We ignore it in the program above
1294because the direct solver we use does not need a tolerance and instead
1295solves the linear system exactly (up to round-off, of course), but iterative
1296solvers could make use of
this kind of information -- and, in fact, should.
1299<a name=
"PlainProg"></a>
1300<h1> The plain program</h1>
1301@include
"step-77.cc"
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
__global__ void set(Number *val, const Number s, const size_type N)
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ 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.)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
VectorType::value_type * end(VectorType &V)
int(&) functions(const void *v1, const void *v2)
static constexpr double PI
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation