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)
705 * cell_residual(i) = (fe_values.shape_grad(i, q) // \nabla \phi_i
707 * * evaluation_point_gradients[q] // * u_n
708 * * fe_values.JxW(q)); // * dx
711 * cell->get_dof_indices(local_dof_indices);
712 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
713 * residual(local_dof_indices[i]) += cell_residual(i);
716 * hanging_node_constraints.condense(residual);
718 * for (const types::global_dof_index i :
719 * DoFTools::extract_boundary_dofs(dof_handler))
722 * for (const types::global_dof_index i :
723 * DoFTools::extract_hanging_node_dofs(dof_handler))
726 * std::cout << " norm=" << residual.l2_norm() << std::endl;
734 * <a name="SolvinglinearsystemswiththeJacobianmatrix"></a>
735 * <h4>Solving linear systems with the Jacobian matrix</h4>
739 * Next up is the function that implements the solution of a linear system
740 * with the Jacobian matrix. Since we have already factored the matrix when we
741 * built the matrix, solving a linear system comes down to applying the
742 * inverse matrix to the given right hand side vector: This is what the
743 * SparseDirectUMFPACK::vmult() function does that we use here. Following
744 * this, we have to make sure that we also address the values of hanging nodes
745 * in the solution vector, and this is done using
746 * AffineConstraints::distribute().
750 * The function takes an additional, but unused, argument `tolerance` that
751 * indicates how accurately we have to solve the linear system. The meaning of
752 * this argument is discussed in the introduction in the context of the
753 * "Eisenstat Walker trick", but since we are using a direct rather than an
754 * iterative solver, we are not using this opportunity to solve linear systems
759 * void MinimalSurfaceProblem<dim>::solve(const Vector<double> &rhs,
760 * Vector<double> & solution,
761 * const double /*tolerance*/)
763 * TimerOutput::Scope t(computing_timer, "linear system solve");
765 * std::cout << " Solving linear system" << std::endl;
767 * jacobian_matrix_factorization->vmult(solution, rhs);
769 * hanging_node_constraints.distribute(solution);
777 * <a name="Refiningthemeshsettingboundaryvaluesandgeneratinggraphicaloutput"></a>
778 * <h4>Refining the mesh, setting boundary values, and generating graphical output</h4>
782 * The following three functions are again simply copies of the ones in
783 * @ref step_15 "step-15":
787 * void MinimalSurfaceProblem<dim>::refine_mesh()
789 * Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
791 * KellyErrorEstimator<dim>::estimate(
793 * QGauss<dim - 1>(fe.degree + 1),
794 * std::map<types::boundary_id, const Function<dim> *>(),
796 * estimated_error_per_cell);
798 * GridRefinement::refine_and_coarsen_fixed_number(triangulation,
799 * estimated_error_per_cell,
803 * triangulation.prepare_coarsening_and_refinement();
805 * SolutionTransfer<dim> solution_transfer(dof_handler);
806 * solution_transfer.prepare_for_coarsening_and_refinement(current_solution);
808 * triangulation.execute_coarsening_and_refinement();
810 * dof_handler.distribute_dofs(fe);
812 * Vector<double> tmp(dof_handler.n_dofs());
813 * solution_transfer.interpolate(current_solution, tmp);
814 * current_solution = std::move(tmp);
816 * hanging_node_constraints.clear();
818 * DoFTools::make_hanging_node_constraints(dof_handler,
819 * hanging_node_constraints);
820 * hanging_node_constraints.close();
822 * hanging_node_constraints.distribute(current_solution);
824 * set_boundary_values();
826 * setup_system(/*initial_step=*/false);
832 * void MinimalSurfaceProblem<dim>::set_boundary_values()
834 * std::map<types::global_dof_index, double> boundary_values;
835 * VectorTools::interpolate_boundary_values(dof_handler,
837 * BoundaryValues<dim>(),
839 * for (const auto &boundary_value : boundary_values)
840 * current_solution(boundary_value.first) = boundary_value.second;
842 * hanging_node_constraints.distribute(current_solution);
848 * void MinimalSurfaceProblem<dim>::output_results(
849 * const unsigned int refinement_cycle)
851 * TimerOutput::Scope t(computing_timer, "graphical output");
853 * DataOut<dim> data_out;
855 * data_out.attach_dof_handler(dof_handler);
856 * data_out.add_data_vector(current_solution, "solution");
857 * data_out.build_patches();
859 * const std::string filename =
860 * "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
861 * std::ofstream output(filename);
862 * data_out.write_vtu(output);
870 * <a name="Therunfunctionandtheoveralllogicoftheprogram"></a>
871 * <h4>The run() function and the overall logic of the program</h4>
875 * The only function that *really* is interesting in this program is the one
876 * that drives the overall algorithm of starting on a coarse mesh, doing some
877 * mesh refinement cycles, and on each mesh using KINSOL to find the solution
878 * of the nonlinear algebraic equation we obtain from discretization on this
879 * mesh. The `refine_mesh()` function above makes sure that the solution on
880 * one mesh is used as the starting guess on the next mesh. We also use a
881 * TimerOutput object to measure how much time every operation on each mesh
882 * costs, and reset the timer at the beginning of each cycle.
886 * As discussed in the introduction, it is not necessary to solve problems on
887 * coarse meshes particularly accurately since these will only solve as
888 * starting guesses for the next mesh. As a consequence, we will use a target
890 * @f$\tau=10^{-3} \frac{1}{10^k}@f$ for the @f$k@f$th mesh refinement cycle.
894 * All of this is encoded in the first part of this function:
898 * void MinimalSurfaceProblem<dim>::run()
900 * GridGenerator::hyper_ball(triangulation);
901 * triangulation.refine_global(2);
903 * setup_system(/*initial_step=*/true);
904 * set_boundary_values();
906 * for (unsigned int refinement_cycle = 0; refinement_cycle < 6;
907 * ++refinement_cycle)
909 * computing_timer.reset();
910 * std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
912 * if (refinement_cycle != 0)
915 * const double target_tolerance = 1e-3 * std::pow(0.1, refinement_cycle);
916 * std::cout << " Target_tolerance: " << target_tolerance << std::endl
921 * This is where the fun starts. At the top we create the KINSOL solver
922 * object and feed it with an object that encodes a number of additional
923 * specifics (of which we only change the nonlinear tolerance we want to
924 * reach; but you might want to look into what other members the
925 * SUNDIALS::KINSOL::AdditionalData class has and play with them).
929 * typename SUNDIALS::KINSOL<Vector<double>>::AdditionalData
931 * additional_data.function_tolerance = target_tolerance;
933 * SUNDIALS::KINSOL<Vector<double>> nonlinear_solver(additional_data);
937 * Then we have to describe the operations that were already mentioned
938 * in the introduction. In essence, we have to teach KINSOL how to (i)
939 * resize a vector to the correct size, (ii) compute the residual
940 * vector, (iii) compute the Jacobian matrix (during which we also
941 * compute its factorization), and (iv) solve a linear system with the
946 * All four of these operations are represented by member variables of
947 * the SUNDIALS::KINSOL class that are of type `std::function`, i.e.,
948 * they are objects to which we can assign a pointer to a function or,
949 * as we do here, a "lambda function" that takes the appropriate
950 * arguments and returns the appropriate information. By convention,
951 * KINSOL wants that functions doing something nontrivial return an
952 * integer where zero indicates success. It turns out that we can do
953 * all of this in just 25 lines of code.
957 * (If you're not familiar what
"lambda functions" are, take
958 * a look at @ref step_12
"step-12" or at the
959 * [wikipedia page](https:
961 * wants to define a function with a certain
set of
962 * arguments, but (i) not make it a named
functions because,
963 * typically, the function is used in only
one place and it
964 * seems unnecessary to give it a global name; and (ii) that
965 * the function has access to some of the variables that
966 * exist at the place where it is defined, including member
968 * ultimately quite useful.)
972 * At the very
end of the code block we then tell KINSOL to go to work
973 * and solve our problem. The member
functions called from the
974 *
'residual',
'setup_jacobian', and
'solve_jacobian_system' functions
975 * will then print output to screen that allows us to follow along
976 * with the progress of the program.
980 * x.reinit(dof_handler.n_dofs());
983 * nonlinear_solver.residual =
986 * compute_residual(evaluation_point, residual);
991 * nonlinear_solver.setup_jacobian =
994 * compute_and_factorize_jacobian(current_u);
999 * nonlinear_solver.solve_with_jacobian = [&](
const Vector<double> &rhs,
1001 *
const double tolerance) {
1002 * this->solve(rhs, dst, tolerance);
1007 * nonlinear_solver.solve(current_solution);
1012 * The rest is then just house-keeping: Writing data to a file
for
1013 * visualizing, and showing a summary of the timing collected so that we
1014 * can interpret how
long each operation has taken, how often it was
1018 * output_results(refinement_cycle);
1020 * computing_timer.print_summary();
1022 * std::cout << std::endl;
1032 *
using namespace Step77;
1034 * MinimalSurfaceProblem<2> laplace_problem_2d;
1035 * laplace_problem_2d.run();
1037 *
catch (std::exception &exc)
1039 * std::cerr << std::endl
1041 * <<
"----------------------------------------------------"
1043 * std::cerr <<
"Exception on processing: " << std::endl
1044 * << exc.what() << std::endl
1045 * <<
"Aborting!" << std::endl
1046 * <<
"----------------------------------------------------"
1053 * std::cerr << std::endl
1055 * <<
"----------------------------------------------------"
1057 * std::cerr <<
"Unknown exception!" << std::endl
1058 * <<
"Aborting!" << std::endl
1059 * <<
"----------------------------------------------------"
1066<a name=
"Results"></a><h1>Results</h1>
1069When running the program, you get output that looks like
this:
1071Mesh refinement step 0
1072 Target_tolerance: 0.001
1074 Computing residual vector...
norm=0.231202
1075 Computing Jacobian
matrix
1076 Factorizing Jacobian
matrix
1077 Solving linear system
1078 Computing residual vector...
norm=0.231202
1079 Computing residual vector...
norm=0.171585
1080 Solving linear system
1081 Computing residual vector...
norm=0.171585
1082 Computing residual vector...
norm=0.127245
1083 Computing residual vector...
norm=0.0796471
1084 Solving linear system
1085 Computing residual vector...
norm=0.0796471
1086 Computing residual vector...
norm=0.0625301
1087 Solving linear system
1088 Computing residual vector...
norm=0.0625301
1089 Computing residual vector...
norm=0.0498864
1090 Solving linear system
1091 Computing residual vector...
norm=0.0498864
1092 Computing residual vector...
norm=0.0407765
1093 Solving linear system
1094 Computing residual vector...
norm=0.0407765
1095 Computing residual vector...
norm=0.0341589
1096 Solving linear system
1097 Computing residual vector...
norm=0.0341589
1098 Computing residual vector...
norm=0.0292867
1099 Solving linear system
1100 Computing residual vector...
norm=0.0292867
1101 Computing residual vector...
norm=0.0256309
1102 Computing residual vector...
norm=0.0223448
1103 Solving linear system
1104 Computing residual vector...
norm=0.0223448
1105 Computing residual vector...
norm=0.0202797
1106 Computing residual vector...
norm=0.0183817
1107 Solving linear system
1108 Computing residual vector...
norm=0.0183817
1109 Computing residual vector...
norm=0.0170464
1110 Computing residual vector...
norm=0.0157967
1111 Computing Jacobian
matrix
1112 Factorizing Jacobian
matrix
1113 Solving linear system
1114 Computing residual vector...
norm=0.0157967
1115 Computing residual vector...
norm=0.0141572
1116 Computing residual vector...
norm=0.012657
1117 Solving linear system
1118 Computing residual vector...
norm=0.012657
1119 Computing residual vector...
norm=0.0116863
1120 Computing residual vector...
norm=0.0107696
1121 Solving linear system
1122 Computing residual vector...
norm=0.0107696
1123 Computing residual vector...
norm=0.0100986
1124 Computing residual vector...
norm=0.00944829
1125 Computing residual vector...
norm=0.00822576
1126 Solving linear system
1127 Computing residual vector...
norm=0.00822576
1128 Computing residual vector...
norm=0.00781983
1129 Computing residual vector...
norm=0.00741619
1130 Computing residual vector...
norm=0.00661792
1131 Solving linear system
1132 Computing residual vector...
norm=0.00661792
1133 Computing residual vector...
norm=0.00630571
1134 Computing residual vector...
norm=0.00599457
1135 Computing residual vector...
norm=0.00537663
1136 Solving linear system
1137 Computing residual vector...
norm=0.00537663
1138 Computing residual vector...
norm=0.00512813
1139 Computing residual vector...
norm=0.00488033
1140 Computing residual vector...
norm=0.00438751
1141 Computing residual vector...
norm=0.00342052
1142 Solving linear system
1143 Computing residual vector...
norm=0.00342052
1144 Computing residual vector...
norm=0.00326581
1145 Computing residual vector...
norm=0.00311176
1146 Computing residual vector...
norm=0.00280617
1147 Computing residual vector...
norm=0.00220992
1148 Solving linear system
1149 Computing residual vector...
norm=0.00220992
1150 Computing residual vector...
norm=0.00209976
1151 Computing residual vector...
norm=0.00199943
1152 Solving linear system
1153 Computing residual vector...
norm=0.00199942
1154 Computing residual vector...
norm=0.00190953
1155 Computing residual vector...
norm=0.00182005
1156 Computing residual vector...
norm=0.00164259
1157 Computing residual vector...
norm=0.00129652
1160+---------------------------------------------+------------+------------+
1161| Total wallclock time elapsed since start | 0.192s | |
1163| Section | no. calls | wall time | % of total |
1164+---------------------------------+-----------+------------+------------+
1165| assembling the Jacobian | 2 | 0.0141s | 7.4% |
1166| assembling the residual | 61 | 0.168s | 88% |
1167| factorizing the Jacobian | 2 | 0.0016s | 0.83% |
1168| graphical output | 1 | 0.00385s | 2% |
1169| linear system solve | 19 | 0.0013s | 0.68% |
1170+---------------------------------+-----------+------------+------------+
1173Mesh refinement step 1
1174 Target_tolerance: 0.0001
1176 Computing residual vector...
norm=0.0883422
1177 Computing Jacobian
matrix
1178 Factorizing Jacobian
matrix
1179 Solving linear system
1180 Computing residual vector...
norm=0.0883422
1181 Computing residual vector...
norm=0.0607066
1182 Solving linear system
1183 Computing residual vector...
norm=0.0607066
1184 Computing residual vector...
norm=0.0437266
1185 Solving linear system
1186 Computing residual vector...
norm=0.0437266
1187 Computing residual vector...
norm=0.0327999
1188 Solving linear system
1189 Computing residual vector...
norm=0.0327999
1190 Computing residual vector...
norm=0.0255418
1191 Solving linear system
1192 Computing residual vector...
norm=0.0255417
1193 Computing residual vector...
norm=0.0206042
1194 Solving linear system
1195 Computing residual vector...
norm=0.0206042
1196 Computing residual vector...
norm=0.0171602
1197 Solving linear system
1198 Computing residual vector...
norm=0.0171602
1199 Computing residual vector...
norm=0.014689
1200 Solving linear system
1205The way
this should be interpreted is most easily explained by looking at
1206the
first few lines of the output on the
first mesh:
1208Mesh refinement step 0
1209Mesh refinement step 0
1210 Target_tolerance: 0.001
1212 Computing residual vector...
norm=0.231202
1213 Computing Jacobian
matrix
1214 Factorizing Jacobian
matrix
1215 Solving linear system
1216 Computing residual vector...
norm=0.231202
1217 Computing residual vector...
norm=0.171585
1218 Solving linear system
1219 Computing residual vector...
norm=0.171585
1220 Computing residual vector...
norm=0.127245
1221 Computing residual vector...
norm=0.0796471
1222 Solving linear system
1223 Computing residual vector...
norm=0.0796471
1226What is happening is
this:
1227- In the
first residual computation, KINSOL computes the residual to see whether
1228 the desired tolerance has been reached. The answer is no, so it requests the
1229 user program to compute the Jacobian
matrix (and the function then also
1231- KINSOL then instructs us to solve a linear system of the form
1232 @f$J_k \, \delta U_k = -F_k@f$ with
this matrix and the previously computed
1234- It is then time to determine how far we want to go in
this direction,
1235 i.e.,
do line search. To
this end, KINSOL
requires us to compute the
1236 residual vector @f$F(U_k + \alpha_k \delta U_k)@f$
for different step lengths
1237 @f$\alpha_k@f$. For the
first step above, it finds an acceptable @f$\alpha_k@f$
1238 after two tries, the
second time around it takes three tries.
1239- Having found a suitable updated solution @f$U_{k+1}@f$, the process is
1240 repeated except now KINSOL is happy with the current Jacobian
matrix
1241 and does not instruct us to re-build the
matrix and its factorization,
1242 and instead asks us to solve a linear system with that same
matrix.
1244The program also writes the solution to a VTU file at the
end
1245of each mesh refinement cycle, and it looks as follows:
1246<table width=
"60%" align=
"center">
1249 <img src=
"https://www.dealii.org/images/steps/developer/step-77.solution.png" alt=
"">
1255The key takeaway messages of
this program are the following:
1257- The solution is the same as the
one we computed in @ref step_15
"step-15", i.e., the
1258 interfaces to %
SUNDIALS' KINSOL package really did what they were supposed
1259 to do. This should not come as a surprise, but the important point is that
1260 we don't have to spend the time implementing the complex algorithms that
1261 underlie advanced nonlinear solvers ourselves.
1263- KINSOL is able to avoid all sorts of operations such as rebuilding the
1264 Jacobian
matrix when that is not actually necessary. Comparing the
1265 number of linear solves in the output above with the number of times
1266 we rebuild the Jacobian and compute its factorization should make it
1267 clear that
this leads to very substantial savings in terms of compute
1268 times, without us having to implement the intricacies of algorithms
1269 that determine when we need to rebuild
this information.
1271<a name=
"extensions"></a>
1272<a name=
"Possibilitiesforextensions"></a><h3> Possibilities
for extensions </h3>
1275For all but the small problems we consider here, a sparse direct solver
1276requires too much time and memory -- we need an iterative solver like
1277we use in many other programs. The trade-off between constructing an
1278expensive preconditioner (say, a geometric or algebraic multigrid method)
1279is different in the current
case, however: Since we can re-use the same
1280matrix for numerous linear solves, we can
do the same
for the preconditioner
1281and putting more work into building a good preconditioner can more easily
1282be justified than
if we used it only
for a single linear solve as
one
1283does
for many other situations.
1285But iterative solvers also afford other opportunities. For example (and as
1286discussed briefly in the introduction), we may not need to solve to
1287very high accuracy (small tolerances) in early nonlinear iterations as
long
1288as we are still far away from the actual solution. This was the basis of the
1289Eisenstat-Walker trick mentioned there.
1291KINSOL provides the function that does the linear solution with a target
1292tolerance that needs to be reached. We ignore it in the program above
1293because the direct solver we use does not need a tolerance and instead
1294solves the linear system exactly (up to round-off, of course), but iterative
1295solvers could make use of
this kind of information -- and, in fact, should.
1298<a name=
"PlainProg"></a>
1299<h1> The plain program</h1>
1300@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.
static const 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.)
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