From: Wolfgang Bangerth Date: Sun, 3 Jun 2012 17:45:29 +0000 (+0000) Subject: Finish step-15. X-Git-Tag: v8.0.0~2531 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fa9c1b2dda7b40255321ef7e0b851bd8e2ddcb27;p=dealii.git Finish step-15. git-svn-id: https://svn.dealii.org/trunk@25599 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/doxygen/tutorial/steps.cmapx b/deal.II/doc/doxygen/tutorial/steps.cmapx index 66bc3f06b4..31b62cfbd6 100644 --- a/deal.II/doc/doxygen/tutorial/steps.cmapx +++ b/deal.II/doc/doxygen/tutorial/steps.cmapx @@ -1,47 +1,47 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/deal.II/doc/doxygen/tutorial/steps.png b/deal.II/doc/doxygen/tutorial/steps.png index 583cbaa611..eab5fdc358 100644 Binary files a/deal.II/doc/doxygen/tutorial/steps.png and b/deal.II/doc/doxygen/tutorial/steps.png differ diff --git a/deal.II/doc/doxygen/tutorial/toc-list.html b/deal.II/doc/doxygen/tutorial/toc-list.html index 09dc3c5097..4a56b86be2 100644 --- a/deal.II/doc/doxygen/tutorial/toc-list.html +++ b/deal.II/doc/doxygen/tutorial/toc-list.html @@ -120,8 +120,8 @@ Step-15 - 1d problems, nonlinear solvers, - transferring a solution across mesh refinement. + A nonlinear elliptic problem: The minimal surface equation. + Newton's method. Transferring a solution across mesh refinement. @@ -280,7 +280,7 @@ Step-43 Advanced techniques for the simulation of porous media flow. - + Step-44 Finite strain hyperelasticity based on a three-field formulation. diff --git a/deal.II/doc/doxygen/tutorial/toc-topics.html b/deal.II/doc/doxygen/tutorial/toc-topics.html index 77646263e8..e2fafe3537 100644 --- a/deal.II/doc/doxygen/tutorial/toc-topics.html +++ b/deal.II/doc/doxygen/tutorial/toc-topics.html @@ -64,6 +64,12 @@ Step-7 Evaluating errors + + + + Step-15 + Nonlinear problems, Newton's method + @@ -387,6 +393,15 @@ + + + + Step-15 + + Minimal surface equation + + + diff --git a/deal.II/doc/news/news.html b/deal.II/doc/news/news.html index 07b1ca7909..0642179bd5 100644 --- a/deal.II/doc/news/news.html +++ b/deal.II/doc/news/news.html @@ -32,6 +32,20 @@ here.

+
+ 2012/06/03: New step-15 +
+
+ Step-15 used to be a program that demonstrated techniques for solving + a rather quirky, nonlinear, 1d problem. It was an interesting problem, + but it served little purpose in teaching deal.II or in demonstrating + nonlinear solution techniques. It has now been replaced by a more + illustrative example, contributed in large part by Sven Wetterauer, in + which we solve the nonlinear elliptic minimal surface equation using a + Newton iteration. +
+ +
2012/05/05: Matrix-free implementation framework
diff --git a/deal.II/examples/step-15/doc/builds-on b/deal.II/examples/step-15/doc/builds-on index 48a0f73876..17402734c7 100644 --- a/deal.II/examples/step-15/doc/builds-on +++ b/deal.II/examples/step-15/doc/builds-on @@ -1 +1 @@ -step-4 +step-6 diff --git a/deal.II/examples/step-15/doc/intro.dox b/deal.II/examples/step-15/doc/intro.dox index b9cabbb931..7c3e096e11 100644 --- a/deal.II/examples/step-15/doc/intro.dox +++ b/deal.II/examples/step-15/doc/intro.dox @@ -234,7 +234,10 @@ some more in the results section.

Summary of the algorithm and testcase

-Overall, the algorithm we use in this program works as follows: +Overall, the program we have here is not unlike step-6 in many regards. The +layout of the main class is essentially the same. On the other hand, the +driving algorithm in the run() function is different and works as +follows:
  1. Start with the function $u^{0}\equiv 0$ and modify it in such a way @@ -282,3 +285,8 @@ Overall, the algorithm we use in this program works as follows: Set $n\leftarrow n+1$ and go to step 2.
+ +The testcase we solve is chosen as follows: We seek to find the solution of +minimal surface over the unit disk $\Omega=\{\mathbf x: \|\mathbf +x\|<1\}\subset \R^2$ where the surface attains the values +$u(x,y)|{\partial\Omega} = g(x,y):=\sin(2 \pi (x+y))$ along the boundary. diff --git a/deal.II/examples/step-15/doc/results.dox b/deal.II/examples/step-15/doc/results.dox index b85402b9ae..8cdb8f3b4c 100644 --- a/deal.II/examples/step-15/doc/results.dox +++ b/deal.II/examples/step-15/doc/results.dox @@ -3,46 +3,35 @@ The output of the program looks as follows: @code -residual:17.1036 -residual:0.0493566 -residual:0.0437197 -residual:0.0388877 -residual:0.0346918 -mesh-refinement:1 -residual:0.0479423 -residual:0.0427299 -residual:0.0381861 -residual:0.03418 -residual:0.0306245 -mesh-refinement:2 +* ******** Initial mesh ******** + Initial residual: 1.53143 + Residual: 1.08746 + Residual: 0.966748 + Residual: 0.859602 + Residual: 0.766462 + Residual: 0.685475 +* ******** Refined mesh 1 ******** + Initial residual: 0.865774 + Residual: 0.759295 + Residual: 0.675281 + Residual: 0.603523 + Residual: 0.540744 + Residual: 0.485238 +* ******** Refined mesh 2 ******** + Initial residual: 0.425581 + Residual: 0.382042 + Residual: 0.343307 + Residual: 0.308718 .... -residual:0.00354605 -residual:0.00318818 -residual:0.00286721 -residual:0.00257911 -residual:0.00232034 -mesh-refinement:7 -residual:0.00211166 -residual:0.00189985 -residual:0.00170948 -residual:0.00153832 -residual:0.00138438 -mesh-refinement:8 -residual:0.00143114 -residual:0.00128729 -residual:0.00115805 -residual:0.00104189 -residual:0.000937454 -mesh-refinement:9 - @endcode +Obviously, the scheme converges, if not very fast. We will come back to +strategies for accelerating the method below. -We can visualize the solution after each set of five Newton +One can visualize the solution after each set of five Newton iterations, i.e., on each of the meshes on which we approximate the solution. This yields the following set of images: -
@@ -103,6 +92,93 @@ and mesh are shown here:

Possibilities for extensions

+The program shows the basic structure of a solver for a nonlinear, stationary +problem. However, it does not converge particularly fast, for good reasons: + +- The program always takes a step size of 0.1. This precludes the rapid, + quadratic convergence for which Newton's method is typically chosen. +- It does not connect the nonlinear iteration with the mesh refinement + iteration. + +Obviously, a better program would have to address these two points. +We will discuss them in the following. + +

Step length control

-

Integrating mesh refinement and nonlinear solvers

+Newton's method has two well known properties: +- It does not converge from arbitrarily chosen starting points. Rather, a + starting point has to be close enough to the solution to guarantee + convergence. However, we can enlarge the area from which Newton's method + converges by damping the iteration using a step length 0<$\alpha^n\le + 1$. +- It exhibits rapid convergence of quadratic order if (i) the step length is + chosen as $\alpha^n=1$, and (ii) it does in fact converge with this choice + of step length. + +A consequence of these two observations is that a successful strategy is to +choose $\alpha^n<1$ for the initial iterations until the iterate has come +close enough to allow for convergence with full step length, at which point we +want to switch to $\alpha^n=1$. The question is how to choose $\alpha^n$ in an +automatic fashion that satisfies these criteria. + +We do not want to review the literature on this topic here, but only briefly +mention that there are two fundamental approaches to the problem: backtracking +line search and trust region methods. The former is more widely used for +partial differential equations and essentially does the following: +- Compute a search direction +- See if the resulting residual of $u^n + \alpha^n\;\delta u^n$ with + $\alpha^n=1$ is "substantially smaller" than that of $u^n$ alone. +- If so, then take $\alpha^n=1$. +- If not, try whether the residual is "substantially smaller" with + $\alpha^n=2/3$. +- If so, then take $\alpha^n=2/3$. +- If not, try whether the residual is "substantially smaller" with + $\alpha^n=(2/3)^2$. +- Etc. +One can of course choose other factors $r, r^2, \ldots$ than the $2/3, +(2/3)^2, \ldots$ chosen above, for $0determine_step_length() is written the way it is to support +exactly this kind of use case. + +Whether we accept a particular step length $\alpha^n$ depends on how we define +"substantially smaller". There are a number of ways to do so, but without +going into detail let us just mention that the most common ones are to use the +Wolfe and Armijo-Goldstein conditions. For these, one can show the following: +- There is always a step length $\alpha^n$ for which the conditions are + satisfied, i.e., the iteration never gets stuck as long as the problem is + convex. +- If we are close enough to the solution, then the conditions allow for + $\alpha^n$, thereby enabling quadratic convergence. + +We will not dwell on this here any further but leave the implementation of +such algorithms as an exercise. We note, however, that when implemented +correctly then it is a common observation that most reasonably nonlinear +problems can be solved in anywhere between 5 and 15 Newton iterations to +engineering accuracy — substantially fewer than we need with the current +version of the program. + + +

Integrating mesh refinement and nonlinear and linear solvers

+ +We currently do exactly 5 iterations on each mesh. But is this optimal? One +could ask the following questions: +- Maybe it is worthwhile doing more iterations on the initial meshes since + there, computations are cheap. +- On the other hand, we do not want to do too many iterations on every mesh: + yes, we could drive the residual to zero on every mesh, but that would only + mean that the nonlinear iteration error is far smaller than the + discretization error. +- Should we use solve the linear systems in each Newton step with higher or + lower accuracy? + +Ultimately, what this boils down to is that we somehow need to couple the +discretization error on the current mesh with the nonlinear residual we want +to achieve with the Newton iterations on a given mesh, and to the linear +iteration we want to achieve with the CG method within each Newton +iterations. + +How to do this is, again, not entirely trivial, and we again leave it as a +future exercise. diff --git a/deal.II/examples/step-15/doc/tooltip b/deal.II/examples/step-15/doc/tooltip index 4e056430e0..70b5c4b8bb 100644 --- a/deal.II/examples/step-15/doc/tooltip +++ b/deal.II/examples/step-15/doc/tooltip @@ -1 +1 @@ -1d problems. A nonlinear problem. +A nonlinear elliptic problem. Newton's method. diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index 41d36fa8b8..a3375ddf0e 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -52,15 +52,20 @@ #include #include - // We will use adaptive mesh refinement between Newton - // interations. To do so, we need to be able to work - // with a solution on the new mesh, although it was - // computed on the old one. The SolutionTransfer - // class transfers the solution to the new mesh. + // We will use adaptive mesh refinement + // between Newton interations. To do so, we + // need to be able to work with a solution on + // the new mesh, although it was computed on + // the old one. The SolutionTransfer class + // transfers the solution from the old to the + // new mesh: #include - // As in previous programs: + // We then open a namepsace for this program + // and import everything from the dealii + // namespace into it, as in previous + // programs: namespace Step15 { using namespace dealii; @@ -68,18 +73,47 @@ namespace Step15 // @sect3{The MinimalSurfaceProblem class template} - // The class template is basically - // the same as in step 6. Four - // additions are made: There are - // two solution vectors, one for - // the Newton update, and one for - // the solution of the original - // pde. Also we need a double for - // the residual of the Newton - // method, an integer, which counts - // the mesh refinements and a bool - // for the boundary condition in - // the first Newton step. + // The class template is basically the same + // as in step-6. Four additions are made: + // - There are two solution vectors, one for + // the Newton update $\delta u^n$, and one + // for the current iterate $u^n$. + // - The setup_system function + // takes an argument that denotes whether + // this is the first time it is called or + // not. The difference is that the first + // time around we need to distributed + // degrees of freedom and set the + // solution vector for $u^n$ to the + // correct size. The following times, the + // function is called after we have + // already done these steps as part of + // refining the mesh in + // refine_mesh. + // - We then also need new functions: + // set_boundary_values() + // takes care of setting the boundary + // values on the solution vector + // correctly, as discussed at the end of + // the + // introduction. compute_residual() + // is a function that computes the norm + // of the nonlinear (discrete) + // residual. We use this function to + // monitor convergence of the Newton + // iteration. The function takes a step + // length $\alpha^n$ as argument to + // compute the residual of $u^n + + // \alpha^n \; \delta u^n$. This is + // something one typically needs for step + // length control, although we will not + // use this feature here. Finally, + // determine_step_length() + // computes the step length $\alpha^n$ in + // each Newton iteration. As discussed in + // the introduction, we here use a fixed + // step length and leave implementing a + // better strategy as an exercise. template class MinimalSurfaceProblem @@ -91,10 +125,10 @@ namespace Step15 void run (); private: - void setup_system (); + void setup_system (const bool initial_step); void assemble_system (); void solve (); - void refine_grid (); + void refine_mesh (); void set_boundary_values (); double compute_residual (const double alpha) const; double determine_step_length () const; @@ -112,29 +146,14 @@ namespace Step15 Vector present_solution; Vector newton_update; Vector system_rhs; - - - - unsigned int refinement; - - // As described in the - // Introduction, the first - // Newton iteration is special, - // because of the boundary - // condition. To implement - // these correctly, there is a - // bool, which is true in the - // first step and false ever - // after. - bool first_step; }; // @sect3{Boundary condition} // The boundary condition is - // implemented just like in step 4. - // It was chosen as $g(x,y)=sin(2 - // \pi (x+y))$ in this example. + // implemented just like in step-4. + // It is chosen as $g(x,y)=\sin(2 + // \pi (x+y))$: template class BoundaryValues : public Function @@ -146,6 +165,7 @@ namespace Step15 const unsigned int component = 0) const; }; + template double BoundaryValues::value (const Point &p, const unsigned int /*component*/) const @@ -170,7 +190,6 @@ namespace Step15 - // template MinimalSurfaceProblem::~MinimalSurfaceProblem () { @@ -179,50 +198,32 @@ namespace Step15 // @sect4{MinimalSurfaceProblem::setup_system} - // As always in the setup-system - // function, we setup the variables - // of the finite element - // method. There are same - // differences to step 6, because - // we don't have to solve one pde - // over all, but one in every - // Newton step. Also the starting - // function has to be setup in the - // first step. + // As always in the setup-system function, + // we setup the variables of the finite + // element method. There are same + // differences to step-6, because there we + // start solving the PDE from scratch in + // every refinement cycle whereas here we + // need to take the solution from the + // previous mesh onto the current + // mesh. Consequently, we can't just reset + // solution vectors. The argument passed to + // this function thus indicates whether we + // can distributed degrees of freedom (plus + // compute constraints) and set the + // solution vector to zero or whether this + // has happened elsewhere already + // (specifically, in + // refine_mesh()). template - void MinimalSurfaceProblem::setup_system () + void MinimalSurfaceProblem::setup_system (const bool initial_step) { - - // This function will be called, - // every time we refine the mesh - // to resize the system matrix, - // Newton update - and right hand - // side vector and to set the - // right values of hanging nodes - // to get a continuous solution. - // But only the first time, the - // starting solution has to be - // initialized. Also the vector - // of the solution will be - // resized in the - // refine_grid - // function, while the vector is - // transferred to the new mesh. - - if (first_step) + if (initial_step) { dof_handler.distribute_dofs (fe); present_solution.reinit (dof_handler.n_dofs()); - // The constraint matrix, - // holding a list of the - // hanging nodes, will be - // setup in the - // refine_grid - // function after refining - // the mesh. - hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints); @@ -232,7 +233,7 @@ namespace Step15 // The remaining parts of the // function are the same as in - // step 6. + // step-6. newton_update.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); @@ -248,14 +249,25 @@ namespace Step15 // @sect4{MinimalSurfaceProblem::assemble_system} - // This function does the same as - // in the previous tutorials. The - // only additional step is the - // correct implementation of the - // boundary condition and the usage - // of the gradients of the old - // solution. - + // This function does the same as in the + // previous tutorials except that now, of + // course, the matrix and right hand side + // functions depend on the previous + // iteration's solution. As discussed in + // the introduction, we need to use zero + // boundary values for the Newton updates; + // we compute them at the end of this + // function. + // + // The top of the function contains the + // usual boilerplate code, setting up the + // objects that allow us to evaluate shape + // functions at quadrature points and + // temporary storage locations for the + // local matrices and vectors, as well as + // for the gradients of the previous + // solution at the quadrature points. We + // then start the loop over all cells: template void MinimalSurfaceProblem::assemble_system () { @@ -265,16 +277,19 @@ namespace Step15 system_rhs = 0; FEValues fe_values (fe, quadrature_formula, - update_gradients | - update_quadrature_points | update_JxW_values); + update_gradients | + update_quadrature_points | + update_JxW_values); - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - Vector cell_rhs (dofs_per_cell); + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); + std::vector > old_solution_gradients(n_q_points); + + std::vector local_dof_indices (dofs_per_cell); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -286,57 +301,68 @@ namespace Step15 fe_values.reinit (cell); - + // For the assembly of the linear + // system, we have to obtain the + // values of the previous solution's + // gradients at the quadrature + // points. There is a standard way of + // doing this: the + // FEValues::get_function function + // takes a vector that represents a + // finite element field defined on a + // DoFHandler, and evaluates the + // gradients of this field at the + // quadrature points of the cell with + // which the FEValues object has last + // been reinitialized. The values of + // the gradients at all quadrature + // points are then written into the + // second argument: + fe_values.get_function_gradients(present_solution, + old_solution_gradients); + + // With this, we can then do the + // integration loop over all + // quadrature points and shape + // functions. Having just computed + // the gradients of the old solution + // in the quadrature points, we are + // able to compute the coefficients + // $a_{n}$ in these points. The + // assembly of the system itself then + // looks similar to what we always do + // with the exception of the + // nonlinear terms, as does copying + // the results from the local objects + // into the global ones: for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { + const double coeff + = 1.0 / std::sqrt(1 + + old_solution_gradients[q_point] * + old_solution_gradients[q_point]); - // To setup up the linear - // system, the gradient of - // the old solution in the - // quadrature points is - // needed. For this purpose - // there is is a function, - // which will write these - // gradients in a vector, - // where every component of - // the vector is a vector - // itself: - - std::vector > gradients(n_q_points); - fe_values.get_function_gradients(present_solution, gradients); - - // Having the gradients of - // the old solution in the - // quadrature points, we - // are able to compute the - // coefficients $a_{n}$ in - // these points. - - const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]); - - // The assembly of the - // system then is the same - // as always, except of the - // damping parameter of the - // Newton method, which we - // set on 0.1 in this case. - - for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int i=0; i boundary_values; - + std::map boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, ZeroFunction(), boundary_values); - MatrixTools::apply_boundary_values (boundary_values, system_matrix, newton_update, @@ -368,103 +398,14 @@ namespace Step15 } - template - double MinimalSurfaceProblem::compute_residual (const double alpha) const - { - const QGauss quadrature_formula(3); - - Vector residual (dof_handler.n_dofs()); - - Vector linearization_point (dof_handler.n_dofs()); - linearization_point = present_solution; - linearization_point.add (alpha, newton_update); - - FEValues fe_values (fe, quadrature_formula, - update_gradients | - update_quadrature_points | update_JxW_values); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - Vector cell_rhs (dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - cell_rhs = 0; - - fe_values.reinit (cell); - - - for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - { - - // To setup up the linear - // system, the gradient of - // the old solution in the - // quadrature points is - // needed. For this purpose - // there is is a function, - // which will write these - // gradients in a vector, - // where every component of - // the vector is a vector - // itself: - - std::vector > gradients(n_q_points); - fe_values.get_function_gradients(linearization_point, gradients); - - // Having the gradients of - // the old solution in the - // quadrature points, we - // are able to compute the - // coefficients $a_{n}$ in - // these points. - - const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]); - - // The assembly of the - // system then is the same - // as always, except of the - // damping parameter of the - // Newton method, which we - // set on 0.1 in this case. - - for (unsigned int i = 0; i < dofs_per_cell; ++i) - cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff - * gradients[q_point] * fe_values.JxW(q_point)); - } - - cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(), - boundary_values); - for (std::map::const_iterator p = boundary_values.begin(); - p != boundary_values.end(); ++p) - residual(p->first) = 0; - - return residual.l2_norm(); - } // @sect4{MinimalSurfaceProblem::solve} - // The solve function is the same - // as always, we just have to - // implement the minimal residual - // method as a solver and apply the - // Newton update to the solution. - + // The solve function is the same as + // always. At the end of the solution + // process we update the current solution + // by setting $u^{n+1}=u^n+\alpha^n\;\delta + // u^n$. template void MinimalSurfaceProblem::solve () { @@ -480,32 +421,23 @@ namespace Step15 hanging_node_constraints.distribute (newton_update); - // In this step, the old solution - // is updated to the new one: const double alpha = determine_step_length(); - std::cout << " step length alpha=" << alpha << std::endl; present_solution.add (alpha, newton_update); } - template - double MinimalSurfaceProblem::determine_step_length() const - { - return 0.1; - } - // @sect4{MinimalSurfaceProblem::refine_grid} - - // The first part of this function - // is the same as in step 6. But - // after refining the mesh we have - // to transfer the old solution to - // the new one, which is done with - // the help of the SolutionTransfer - // class. - + // @sect4{MinimalSurfaceProblem::refine_mesh} + // The first part of this function is the + // same as in step-6... However, after + // refining the mesh we have to transfer + // the old solution to the new one which we + // do with the help of the SolutionTransfer + // class. The process is slightly + // convoluted, so let us describe it in + // detail: template - void MinimalSurfaceProblem::refine_grid () + void MinimalSurfaceProblem::refine_mesh () { Vector estimated_error_per_cell (triangulation.n_active_cells()); @@ -519,248 +451,416 @@ namespace Step15 estimated_error_per_cell, 0.3, 0.03); - // Then we need an additional - // step: if, for example, you - // flag a cell that is once more - // refined than its neighbor, and - // that neighbor is not flagged - // for refinement, we would end - // up with a jump of two - // refinement levels across a - // cell interface. To avoid - // these situations, the library - // will silently also have to - // refine the neighbor cell - // once. It does so by calling - // the + // Then we need an additional step: if, + // for example, you flag a cell that is + // once more refined than its neighbor, + // and that neighbor is not flagged for + // refinement, we would end up with a + // jump of two refinement levels across a + // cell interface. To avoid these + // situations, the library will silently + // also have to refine the neighbor cell + // once. It does so by calling the // Triangulation::prepare_coarsening_and_refinement - // function before actually doing - // the refinement and coarsening. - // This function flags a set of - // additional cells for - // refinement or coarsening, to + // function before actually doing the + // refinement and coarsening. This + // function flags a set of additional + // cells for refinement or coarsening, to // enforce rules like the - // one-hanging-node rule. The - // cells that are flagged for - // refinement and coarsening - // after calling this function + // one-hanging-node rule. The cells that + // are flagged for refinement and + // coarsening after calling this function // are exactly the ones that will // actually be refined or - // coarsened. Since the - // SolutionTransfer class needs - // this information in order to - // store the data from the old - // mesh and transfer to the new - // one. - + // coarsened. Usually, you don't have to + // do this by hand + // (Triangulation::execute_coarsening_and_refinement + // does this for you). However, we need + // to initialize the SolutionTransfer + // class and it needs to know the final + // set of cells that will be coarsened or + // refined in order to store the data + // from the old mesh and transfer to the + // new one. Thus, we call the function by + // hand: triangulation.prepare_coarsening_and_refinement (); // With this out of the way, we - // initialize a SolutionTransfer - // object with the present - // DoFHandler and attach the - // solution vector to it: - + // initialize a SolutionTransfer object + // with the present DoFHandler and attach + // the solution vector to it, followed by + // doing the actual refinement and + // distribution of degrees of freedom on + // the new mesh SolutionTransfer solution_transfer(dof_handler); solution_transfer.prepare_for_coarsening_and_refinement(present_solution); - // Then we do the actual - // refinement, and distribute - // degrees of freedom on the new - // mesh: - triangulation.execute_coarsening_and_refinement(); - dof_handler.distribute_dofs(fe); - // Finally, we retrieve the old - // solution interpolated to the - // new mesh. Since the - // SolutionTransfer function does - // not actually store the values - // of the old solution, but - // rather indices, we need to - // preserve the old solution - // vector until we have gotten - // the new interpolated - // values. Thus, we have the new - // values written into a - // temporary vector, and only - // afterwards write them into the - // solution vector object: + dof_handler.distribute_dofs(fe); + // Finally, we retrieve the old solution + // interpolated to the new mesh. Since + // the SolutionTransfer function does not + // actually store the values of the old + // solution, but rather indices, we need + // to preserve the old solution vector + // until we have gotten the new + // interpolated values. Thus, we have the + // new values written into a temporary + // vector, and only afterwards write them + // into the solution vector object. Once + // we have this solution we have to make + // sure that the $u^n$ we now have + // actually has the correct boundary + // values. As explained at the end of the + // introduction, this is not + // automatically the case even if the + // solution before refinement had the + // correct boundary values, and so we + // have to explicitly make sure that it + // now has: Vector tmp(dof_handler.n_dofs()); - solution_transfer.interpolate(present_solution,tmp); - present_solution=tmp; + solution_transfer.interpolate(present_solution, tmp); + present_solution = tmp; set_boundary_values (); - // On the new mesh, there are - // different hanging nodes, which - // shall be enlisted in a matrix - // like before. To ensure there - // are no hanging nodes of the - // old mesh in the matrix, it's - // first cleared: + // On the new mesh, there are different + // hanging nodes, which we have to + // compute again. To ensure there are no + // hanging nodes of the old mesh in the + // object, it's first cleared. To be on + // the safe side, we then also make sure + // that the current solution's vector + // entries satisfy the hanging node + // constraints: + hanging_node_constraints.clear(); - // After doing so, the hanging - // nodes of the new mesh can be - // enlisted in the matrix, like - // before. Calling the - // setup_system - // function in the - // run function - // again after this, the hanging - // nodes don't have to be - // enlisted there once more. - - DoFTools::make_hanging_node_constraints(dof_handler, hanging_node_constraints); + DoFTools::make_hanging_node_constraints(dof_handler, + hanging_node_constraints); hanging_node_constraints.close(); - hanging_node_constraints.distribute(present_solution); + + hanging_node_constraints.distribute (present_solution); + + // We end the function by updating all + // the remaining data structures, + // indicating to + // setup_dofs() that this is + // not the first go-around and that it + // needs to preserve the content of the + // solution vector: + setup_system (false); } + + // @sect4{MinimalSurfaceProblem::set_boundary_values} + + // The next function ensures that the + // solution vector's entries respect the + // boundary values for our problem. Having + // refined the mesh (or just started + // computations), there might be new nodal + // points on the boundary. These have + // values that are simply interpolated from + // the previous mesh (or are just zero), + // instead of the correct boundary + // values. This is fixed up by setting all + // boundary nodes explicit to the right + // value: template void MinimalSurfaceProblem::set_boundary_values () { - // Having refined the mesh, there - // might be new nodal points on - // the boundary. These have just - // interpolated values, but not - // the right boundary - // values. This is fixed up, by - // setting all boundary nodals - // explicit to the right value: - - std::map boundary_values2; - VectorTools::interpolate_boundary_values(dof_handler, 0, - BoundaryValues(), boundary_values2); - for (std::map::const_iterator p = - boundary_values2.begin(); p != boundary_values2.end(); ++p) + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + boundary_values); + for (std::map::const_iterator + p = boundary_values.begin(); + p != boundary_values.end(); ++p) present_solution(p->first) = p->second; } - // @sect4{MinimalSurfaceProblem::run} - // In the run function, the first - // grid is build. Also in this - // function, the Newton iteration - // is implemented. + // @sect4{MinimalSurfaceProblem::compute_residual} + + // In order to monitor convergence, we need + // a way to compute the norm of the + // (discrete) residual, i.e., the norm of + // the vector + // $\left$ with + // $F(u)=-\nabla \cdot \left( + // \frac{1}{\sqrt{1+|\nabla u|^{2}}}\nabla + // u \right)$ as discussed in the + // introduction. It turns out that + // (although we don't use this feature in + // the current version of the program) one + // needs to compute the residual + // $\left$ + // when determining optimal step lengths, + // and so this is what we implement here: + // the function takes the step length + // $\alpha^n$ as an argument. The original + // functionality is of course obtained by + // passing a zero as argument. + // + // In the function below, we first set up a + // vector for the residual, and then a + // vector for the evaluation point + // $u^n+\alpha^n\;\delta u^n$. This is + // followed by the same boilerplate code we + // use for all integration operations: template - void MinimalSurfaceProblem::run () + double MinimalSurfaceProblem::compute_residual (const double alpha) const { + Vector residual (dof_handler.n_dofs()); - // The integer refinement counts - // the mesh - // refinements. Obviously - // starting the program, it - // should be zero. - refinement=0; - first_step=true; - - // As described in the - // introduction, the domain is a - // unitball around the - // origin. The Mesh is globally - // refined two times, not to - // start on the coarse mesh, - // which consists only of five - // cells. + Vector evaluation_point (dof_handler.n_dofs()); + evaluation_point = present_solution; + evaluation_point.add (alpha, newton_update); - GridGenerator::hyper_ball (triangulation); - static const HyperBallBoundary boundary; - triangulation.set_boundary (0, boundary); - triangulation.refine_global(2); + const QGauss quadrature_formula(3); + FEValues fe_values (fe, quadrature_formula, + update_gradients | + update_quadrature_points | + update_JxW_values); - // The Newton iteration starts - // here. During the first step, - // there is no residual computed, - // so the bool is needed here to - // enter the iteration - // scheme. Later the Newton - // method will continue until the - // residual is less than - // $10^{-3}$. + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); - double previous_res = 0; - while (first_step || (previous_res>1e-3)) + Vector cell_rhs (dofs_per_cell); + std::vector > gradients(n_q_points); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) { + cell_rhs = 0; + fe_values.reinit (cell); + + // The actual computation is much as + // in + // assemble_system(). We + // first evaluate the gradients of + // $u^n+\alpha^n\,\delta u^n$ at the + // quadrature points, then compute + // the coefficient $a_n$, and then + // plug it all into the formula for + // the residual: + fe_values.get_function_gradients (evaluation_point, + gradients); - // In the first step, we - // compute the solution on - // the two times globally - // refined mesh. After that - // the mesh will be refined - // adaptively, in order to - // not get too many - // cells. The refinement is - // the first thing done every - // time we restart the - // process in the while-loop. - if (!first_step) + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_dofs (dof_handler.n_dofs()); + DoFTools::extract_boundary_dofs (dof_handler, + std::vector(1,true), + boundary_dofs); + for (unsigned int i=0; isetup_system - // function. - setup_system(); + // @sect4{MinimalSurfaceProblem::determine_step_length} - if (first_step) - set_boundary_values (); + // As discussed in the introduction, + // Newton's method frequently does not + // converge if we always take full steps, + // i.e., compute $u^{n+1}=u^n+\delta + // u^n$. Rather, one needs a damping + // parameter (step length) $\alpha^n$ and + // set $u^{n+1}=u^n+\alpha^n\; delta + // u^n$. This function is the one called to + // compute $\alpha^n$. + // + // Here, we simply always return 0.1. This + // is of course a sub-optimal choice: + // ideally, what one wants is that the step + // size goes to one as we get closer to the + // solution, so that we get to enjoy the + // rapid quadratic convergence of Newton's + // method. We will discuss better + // strategies below in the results section. + template + double MinimalSurfaceProblem::determine_step_length() const + { + return 0.1; + } + + + + // @sect4{MinimalSurfaceProblem::run} + + // In the run function, we build the first + // grid and then have the top-level logic + // for the Newton iteration. The function + // has two variables, one that indicates + // whether this is the first time we solve + // for a Newton update and one that + // indicates the refinement level of the + // mesh: + template + void MinimalSurfaceProblem::run () + { + unsigned int refinement = 0; + bool first_step = true; + + // As described in the introduction, the + // domain is the unit disk around the + // origin, created in the same way as + // shown in step-6. The mesh is globally + // refined twice followed later on by + // several adaptive cycles: + GridGenerator::hyper_ball (triangulation); + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + triangulation.refine_global(2); - // On every mesh there are - // done five Newton steps, in - // order to get a better - // solution, before the mesh - // gets too fine and the - // computations take more - // time. - std::cout<<"initial residual:"<1e-3)) + { + if (first_step == true) + { + std::cout << "******** Initial mesh " + << " ********" + << std::endl; - for(unsigned int i=0; i<5;++i) + setup_system (true); + set_boundary_values (); + } + else { + ++refinement; + std::cout << "******** Refined mesh " << refinement + << " ********" + << std::endl; - // In every Newton step - // the system matrix and - // the right hand side - // have to be computed. + refine_mesh(); + } + // On every mesh we do exactly five + // Newton steps. We print the initial + // residual here and then start the + // iterations on this mesh. + // + // In every Newton step the system + // matrix and the right hand side + // have to be computed first, after + // which we store the norm of the + // right hand side as the residual to + // check against when deciding + // whether to stop the iterations. We + // then solve the linear system (the + // function also updates + // $u^{n+1}=u^n+\alpha^n\;\delta + // u^n$) and output the residual at + // the end of this Newton step: + std::cout << " Initial residual: " + << compute_residual(0) + << std::endl; + + for (unsigned int inner_iteration=0; inner_iteration<5; ++inner_iteration) + { assemble_system (); previous_res = system_rhs.l2_norm(); solve (); - first_step=false; - std::cout<<"residual:"< data_out; data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (newton_update, "update"); data_out.add_data_vector (present_solution, "solution"); + data_out.add_data_vector (newton_update, "update"); data_out.build_patches (); - const std::string filename = "solution-" + Utilities::int_to_string (refinement, 2) + ".vtk"; + const std::string filename = "solution-" + + Utilities::int_to_string (refinement, 2) + + ".vtk"; std::ofstream output (filename.c_str()); data_out.write_vtk (output); @@ -770,7 +870,7 @@ namespace Step15 // @sect4{The main function} - // Finally the main function, this + // Finally the main function. This // follows the scheme of all other // main functions: int main () diff --git a/deal.II/examples/step-41/doc/builds-on b/deal.II/examples/step-41/doc/builds-on index 48a0f73876..78300ce0a2 100644 --- a/deal.II/examples/step-41/doc/builds-on +++ b/deal.II/examples/step-41/doc/builds-on @@ -1 +1 @@ -step-4 +step-15