]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finish step-15.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 3 Jun 2012 17:45:29 +0000 (17:45 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 3 Jun 2012 17:45:29 +0000 (17:45 +0000)
git-svn-id: https://svn.dealii.org/trunk@25599 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/doxygen/tutorial/steps.cmapx
deal.II/doc/doxygen/tutorial/steps.png
deal.II/doc/doxygen/tutorial/toc-list.html
deal.II/doc/doxygen/tutorial/toc-topics.html
deal.II/doc/news/news.html
deal.II/examples/step-15/doc/builds-on
deal.II/examples/step-15/doc/intro.dox
deal.II/examples/step-15/doc/results.dox
deal.II/examples/step-15/doc/tooltip
deal.II/examples/step-15/step-15.cc
deal.II/examples/step-41/doc/builds-on

index 66bc3f06b4cc166da3244bed75c8bede3fb3e1dd..31b62cfbd6a6015216ef57c649431a09981694f8 100644 (file)
@@ -1,47 +1,47 @@
 <map id="StepsMap" name="StepsMap">
-<area shape="poly" id="node1" href="../deal.II/step_1.html" title="Creating a grid. Refining it. Writing it to a file" alt="" coords="548,60,548,28,525,5,493,5,471,28,471,60,493,83,525,83"/>
-<area shape="poly" id="node2" href="../deal.II/step_2.html" title="Assigning degrees of freedom to a grid." alt="" coords="548,185,548,153,525,131,493,131,471,153,471,185,493,208,525,208"/>
-<area shape="poly" id="node3" href="../deal.II/step_3.html" title="Solving Poisson&#39;s equation." alt="" coords="548,311,548,279,525,256,493,256,471,279,471,311,493,333,525,333"/>
-<area shape="poly" id="node4" href="../deal.II/step_4.html" title="Dimension independent programming. Boundary conditions." alt="" coords="548,436,548,404,525,381,493,381,471,404,471,436,493,459,525,459"/>
-<area shape="poly" id="node5" href="../deal.II/step_5.html" title="Reading a grid from disk. Computations on successively refined grids." alt="" coords="437,561,437,529,415,507,383,507,360,529,360,561,383,584,415,584"/>
-<area shape="rect" id="node10" href="../deal.II/step_10.html" title="Higher order mappings." alt="" coords="239,529,275,562"/>
-<area shape="rect" id="node15" href="../deal.II/step_15.html" title="1d problems. A nonlinear problem." alt="" coords="299,529,335,562"/>
-<area shape="rect" id="node20" href="../deal.II/step_20.html" title="Mixed finite elements for the mixed Laplacian. Block solvers." alt="" coords="462,533,497,558"/>
-<area shape="rect" id="node23" href="../deal.II/step_23.html" title="Time dependent problems. The wave equation." alt="" coords="146,533,181,558"/>
-<area shape="rect" id="node28" href="../deal.II/step_29.html" title="A complex&#45;valued Helmholtz equation. Sparse direct solvers." alt="" coords="522,529,557,562"/>
-<area shape="rect" id="node33" href="../deal.II/step_34.html" title="Boundary element methods for potential flow." alt="" coords="582,533,617,558"/>
-<area shape="rect" id="node35" href="../deal.II/step_36.html" title="Finding eigenvalues of the Schr&ouml;dinger equation." alt="" coords="642,529,677,562"/>
-<area shape="rect" id="node37" href="../deal.II/step_38.html" title="Solve the Laplace Beltrami operator on a Half Sphere." alt="" coords="625,654,660,687"/>
-<area shape="rect" id="node40" href="../deal.II/step_41.html" title="Solving the obstacle problem (a variational inequality)" alt="" coords="754,533,789,558"/>
-<area shape="poly" id="node6" href="../deal.II/step_6.html" title="Adaptive local refinement. Higher order elements" alt="" coords="408,687,408,655,385,632,353,632,331,655,331,687,353,709,385,709"/>
-<area shape="rect" id="node7" href="../deal.II/step_7.html" title="Helmholtz equation. Computing errors. Boundary integrals." alt="" coords="153,758,185,791"/>
-<area shape="rect" id="node8" href="../deal.II/step_8.html" title="Systems of PDE. Elasticity." alt="" coords="442,758,475,791"/>
-<area shape="rect" id="node9" href="../deal.II/step_9.html" title="Advection equation. Multithreading. Refinement criteria." alt="" coords="499,758,532,791"/>
-<area shape="rect" id="node13" href="../deal.II/step_13.html" title="Modularity. Software design." alt="" coords="41,758,76,791"/>
-<area shape="rect" id="node16" href="../deal.II/step_16.html" title="Multigrid on adaptive meshes." alt="" coords="101,841,136,874"/>
-<area shape="rect" id="node22" href="../deal.II/step_22.html" title="The Stokes equation on adaptive meshes." alt="" coords="609,762,644,787"/>
-<area shape="rect" id="node26" href="../deal.II/step_27.html" title="hp&#45;adaptive finite element methods." alt="" coords="322,758,357,791"/>
-<area shape="rect" id="node27" href="../deal.II/step_28.html" title="Handling multiple meshes at the same time. Neutron transport." alt="" coords="210,758,245,791"/>
-<area shape="rect" id="node38" href="../deal.II/step_39.html" title="Interior Penalty for the Laplace equation. Adaptive refinement. Multigrid." alt="" coords="265,923,300,957"/>
-<area shape="rect" id="node39" href="../deal.II/step_40.html" title="Solving the Laplace equation on adaptive meshes on thousands of processors." alt="" coords="553,923,588,957"/>
-<area shape="rect" id="node43" href="../deal.II/step_45.html" title="Periodic boundary conditions" alt="" coords="382,758,417,791"/>
-<area shape="rect" id="node12" href="../deal.II/step_12.html" title="Discontinuous Galerkin for linear advection." alt="" coords="205,841,240,874"/>
-<area shape="rect" id="node17" href="../deal.II/step_17.html" title="Parallel computing using MPI. Using PETSc." alt="" coords="497,845,532,870"/>
-<area shape="rect" id="node42" href="../deal.II/step_44.html" title="Quasi&#45;static finite&#45;strain elasticity" alt="" coords="378,1010,413,1035"/>
-<area shape="rect" id="node44" href="../deal.II/step_46.html" title="Coupling different physical models (flow, elasticity) in different parts of the domain" alt="" coords="377,841,412,874"/>
-<area shape="rect" id="node11" href="../deal.II/step_11.html" title="Higher order mappings. Dealing with constraints." alt="" coords="239,654,275,687"/>
-<area shape="rect" id="node29" href="../deal.II/step_30.html" title="Anisotropic refinement for DG methods." alt="" coords="205,923,240,957"/>
-<area shape="rect" id="node32" href="../deal.II/step_33.html" title="Hyperbolic conservation laws: the Euler equations of gas dynamics." alt="" coords="145,927,180,953"/>
-<area shape="rect" id="node14" href="../deal.II/step_14.html" title="Duality based error estimates. Adaptivity." alt="" coords="41,841,76,874"/>
-<area shape="rect" id="node36" href="../deal.II/step_37.html" title="Matrix&#45;free methods. Multigrid. Cell&#45;based finite element operator." alt="" coords="85,923,120,957"/>
-<area shape="rect" id="node18" href="../deal.II/step_18.html" title="Quasistatic elasticity. More parallel computing." alt="" coords="437,927,472,953"/>
-<area shape="rect" id="node31" href="../deal.II/step_32.html" title="A parallel Boussinesq flow solver for thermal convection in the earth mantle." alt="" coords="553,1010,588,1035"/>
-<area shape="rect" id="node19" href="../deal.II/step_19.html" title="Handling input parameter files. Converting output formats." alt="" coords="437,841,472,874"/>
-<area shape="rect" id="node21" href="../deal.II/step_21.html" title="Two&#45;phase flow in porous media." alt="" coords="565,658,600,683"/>
-<area shape="rect" id="node41" href="../deal.II/step_43.html" title="Efficient ways to solve two&#45;phase flow problems on adaptive meshes in 2d and 3d." alt="" coords="697,927,732,953"/>
-<area shape="rect" id="node30" href="../deal.II/step_31.html" title="Boussinesq flow for thermal convection." alt="" coords="669,845,704,870"/>
-<area shape="rect" id="node34" href="../deal.II/step_35.html" title="A projection solver for the Navier&#45;Stokes equations." alt="" coords="609,845,644,870"/>
-<area shape="rect" id="node24" href="../deal.II/step_24.html" title="The wave equation with absorbing boundary conditions. Extracting point values." alt="" coords="53,658,88,683"/>
-<area shape="rect" id="node25" href="../deal.II/step_25.html" title="The nonlinear sine&#45;Gordon soliton equation" alt="" coords="5,927,40,953"/>
-<area shape="rect" id="node45" href="../deal.II/step_48.html" title="Parallelization via MPI. The wave equation, in linear and nonlinear variants. Mass lumping. Cell&#45;based finite element operator." alt="" coords="85,1006,120,1039"/>
+<area shape="poly" id="node1" href="../deal.II/step_1.html" title="Creating a grid. Refining it. Writing it to a file" alt="" coords="572,60 572,28 549,5 517,5 495,28 495,60 517,83 549,83"/>
+<area shape="poly" id="node2" href="../deal.II/step_2.html" title="Assigning degrees of freedom to a grid." alt="" coords="572,185 572,153 549,131 517,131 495,153 495,185 517,208 549,208"/>
+<area shape="poly" id="node3" href="../deal.II/step_3.html" title="Solving Poisson&#39;s equation." alt="" coords="572,311 572,279 549,256 517,256 495,279 495,311 517,333 549,333"/>
+<area shape="poly" id="node4" href="../deal.II/step_4.html" title="Dimension independent programming. Boundary conditions." alt="" coords="572,436 572,404 549,381 517,381 495,404 495,436 517,459 549,459"/>
+<area shape="poly" id="node5" href="../deal.II/step_5.html" title="Reading a grid from disk. Computations on successively refined grids." alt="" coords="404,561 404,529 381,507 349,507 327,529 327,561 349,584 381,584"/>
+<area shape="rect" id="node10" href="../deal.II/step_10.html" title="Higher order mappings." alt="" coords="487,529,521,562"/>
+<area shape="rect" id="node20" href="../deal.II/step_20.html" title="Mixed finite elements for the mixed Laplacian. Block solvers." alt="" coords="164,531,199,560"/>
+<area shape="rect" id="node23" href="../deal.II/step_23.html" title="Time dependent problems. The wave equation." alt="" coords="428,531,463,560"/>
+<area shape="rect" id="node28" href="../deal.II/step_29.html" title="A complex&#45;valued Helmholtz equation. Sparse direct solvers." alt="" coords="545,529,580,562"/>
+<area shape="rect" id="node33" href="../deal.II/step_34.html" title="Boundary element methods for potential flow." alt="" coords="655,531,689,560"/>
+<area shape="rect" id="node35" href="../deal.II/step_36.html" title="Finding eigenvalues of the Schr&ouml;dinger equation." alt="" coords="713,529,748,562"/>
+<area shape="rect" id="node37" href="../deal.II/step_38.html" title="Solve the Laplace Beltrami operator on a Half Sphere." alt="" coords="655,654,689,687"/>
+<area shape="poly" id="node6" href="../deal.II/step_6.html" title="Adaptive local refinement. Higher order elements" alt="" coords="404,687 404,655 381,632 349,632 327,655 327,687 349,709 381,709"/>
+<area shape="rect" id="node7" href="../deal.II/step_7.html" title="Helmholtz equation. Computing errors. Boundary integrals." alt="" coords="549,758,581,791"/>
+<area shape="rect" id="node8" href="../deal.II/step_8.html" title="Systems of PDE. Elasticity." alt="" coords="151,758,183,791"/>
+<area shape="rect" id="node9" href="../deal.II/step_9.html" title="Advection equation. Multithreading. Refinement criteria." alt="" coords="207,758,239,791"/>
+<area shape="rect" id="node13" href="../deal.II/step_13.html" title="Modularity. Software design." alt="" coords="431,758,465,791"/>
+<area shape="rect" id="node15" href="../deal.II/step_15.html" title="A nonlinear elliptic problem. Newton&#39;s method." alt="" coords="489,758,524,791"/>
+<area shape="rect" id="node16" href="../deal.II/step_16.html" title="Multigrid on adaptive meshes." alt="" coords="657,841,692,874"/>
+<area shape="rect" id="node22" href="../deal.II/step_22.html" title="The Stokes equation on adaptive meshes." alt="" coords="33,760,68,789"/>
+<area shape="rect" id="node26" href="../deal.II/step_27.html" title="hp&#45;adaptive finite element methods." alt="" coords="92,758,127,791"/>
+<area shape="rect" id="node27" href="../deal.II/step_28.html" title="Handling multiple meshes at the same time. Neutron transport." alt="" coords="263,758,297,791"/>
+<area shape="rect" id="node38" href="../deal.II/step_39.html" title="Interior Penalty for the Laplace equation. Adaptive refinement. Multigrid." alt="" coords="597,923,632,957"/>
+<area shape="rect" id="node39" href="../deal.II/step_40.html" title="Solving the Laplace equation on adaptive meshes on thousands of processors." alt="" coords="376,923,411,957"/>
+<area shape="rect" id="node43" href="../deal.II/step_45.html" title="Periodic boundary conditions" alt="" coords="321,758,356,791"/>
+<area shape="rect" id="node12" href="../deal.II/step_12.html" title="Discontinuous Galerkin for linear advection." alt="" coords="548,841,583,874"/>
+<area shape="rect" id="node17" href="../deal.II/step_17.html" title="Parallel computing using MPI. Using PETSc." alt="" coords="267,843,301,872"/>
+<area shape="rect" id="node42" href="../deal.II/step_44.html" title="Quasi&#45;static finite&#45;strain elasticity" alt="" coords="296,1005,331,1035"/>
+<area shape="rect" id="node44" href="../deal.II/step_46.html" title="Coupling different physical models (flow, elasticity) in different parts of the domain" alt="" coords="149,841,184,874"/>
+<area shape="rect" id="node11" href="../deal.II/step_11.html" title="Higher order mappings. Dealing with constraints." alt="" coords="487,654,521,687"/>
+<area shape="rect" id="node29" href="../deal.II/step_30.html" title="Anisotropic refinement for DG methods." alt="" coords="480,923,515,957"/>
+<area shape="rect" id="node32" href="../deal.II/step_33.html" title="Hyperbolic conservation laws: the Euler equations of gas dynamics." alt="" coords="539,925,573,955"/>
+<area shape="rect" id="node14" href="../deal.II/step_14.html" title="Duality based error estimates. Adaptivity." alt="" coords="431,841,465,874"/>
+<area shape="rect" id="node40" href="../deal.II/step_41.html" title="Solving the obstacle problem (a variational inequality)" alt="" coords="489,843,524,872"/>
+<area shape="rect" id="node36" href="../deal.II/step_37.html" title="Matrix&#45;free methods. Multigrid. Cell&#45;based finite element operator." alt="" coords="656,923,691,957"/>
+<area shape="rect" id="node18" href="../deal.II/step_18.html" title="Quasistatic elasticity. More parallel computing." alt="" coords="267,925,301,955"/>
+<area shape="rect" id="node31" href="../deal.II/step_32.html" title="A parallel Boussinesq flow solver for thermal convection in the earth mantle." alt="" coords="212,1005,247,1035"/>
+<area shape="rect" id="node19" href="../deal.II/step_19.html" title="Handling input parameter files. Converting output formats." alt="" coords="208,841,243,874"/>
+<area shape="rect" id="node21" href="../deal.II/step_21.html" title="Two&#45;phase flow in porous media." alt="" coords="33,656,68,685"/>
+<area shape="rect" id="node41" href="../deal.II/step_43.html" title="Efficient ways to solve two&#45;phase flow problems on adaptive meshes in 2d and 3d." alt="" coords="5,925,40,955"/>
+<area shape="rect" id="node30" href="../deal.II/step_31.html" title="Boussinesq flow for thermal convection." alt="" coords="32,843,67,872"/>
+<area shape="rect" id="node34" href="../deal.II/step_35.html" title="A projection solver for the Navier&#45;Stokes equations." alt="" coords="91,843,125,872"/>
+<area shape="rect" id="node24" href="../deal.II/step_24.html" title="The wave equation with absorbing boundary conditions. Extracting point values." alt="" coords="428,656,463,685"/>
+<area shape="rect" id="node25" href="../deal.II/step_25.html" title="The nonlinear sine&#45;Gordon soliton equation" alt="" coords="715,925,749,955"/>
+<area shape="rect" id="node45" href="../deal.II/step_48.html" title="Parallelization via MPI. The wave equation, in linear and nonlinear variants. Mass lumping. Cell&#45;based finite element operator." alt="" coords="656,1005,691,1035"/>
 </map>
index 583cbaa611cd26dfd6f5f4b0ff077c6d728aba09..eab5fdc358f41d9651f2eb8e2c122ee735bdeec4 100644 (file)
Binary files a/deal.II/doc/doxygen/tutorial/steps.png and b/deal.II/doc/doxygen/tutorial/steps.png differ
index 09dc3c5097cf1327926d9a7f1a665901905d3708..4a56b86be2b5cfbccb4e19b7499809643430c88d 100644 (file)
 
   <tr valign="top">
       <td><a href="../../doxygen/deal.II/step_15.html">Step-15</a></td>
-      <td> 1d problems, nonlinear solvers,
-      transferring a solution across mesh refinement.
+      <td> A nonlinear elliptic problem: The minimal surface equation.
+       Newton's method. Transferring a solution across mesh refinement.
       </td></tr>
 
   <tr valign="top">
       <td><a href="../../doxygen/deal.II/step_43.html">Step-43</a></td>
       <td> Advanced techniques for the simulation of porous media flow.
       </td></tr>
-      
+
   <tr valign="top">
       <td><a href="../../doxygen/deal.II/step_44.html">Step-44</a></td>
       <td> Finite strain hyperelasticity based on a three-field formulation.
index 77646263e82e21cdfa6a2e3c4984144f21c6ebdf..e2fafe353704068c1c3ef8a4caafe7bbe050a92f 100644 (file)
     <td><a href="../../doxygen/deal.II/step_7.html">Step-7</a></td>
     <td> Evaluating errors
     </td>
+
+  <tr valign="top">
+    <td></td>
+    <td><a href="../../doxygen/deal.II/step_15.html">Step-15</a></td>
+    <td> Nonlinear problems, Newton's method
+    </td>
   </tr>
 
 
     </td>
   </tr>
 
+  <tr valign="top">
+    <td></td>
+    <td>
+      <a href="../../doxygen/deal.II/step_15.html">Step-15</a>
+    </td>
+    <td> Minimal surface equation
+    </td>
+  </tr>
+
   <tr valign="top">
     <td></td>
     <td>
index 07b1ca7909dbcdda4be3c5622fd633ebeb347d78..0642179bd55a50bbcd65be05aec75316a11734c7 100644 (file)
     <a href="../doxygen/deal.II/changes_after_7_1.html">here.</a></p>
 
     <dl>
+      <dt>
+        <strong>2012/06/03: New step-15</strong>
+      </dt>
+      <dd>
+        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.
+      </dd>
+
+
       <dt>
         <strong>2012/05/05: Matrix-free implementation framework</strong>
       </dt>
index 48a0f738761717ac5d7b461e0b04951a515849bc..17402734c787ca2681ca4585d5db3399ad0e9c49 100644 (file)
@@ -1 +1 @@
-step-4
+step-6
index b9cabbb93115f2f7b6206a88738b39e308622d09..7c3e096e115a038ec0778d9d067ef5994b7f9b9a 100644 (file)
@@ -234,7 +234,10 @@ some more in the <a href="#Results">results section</a>.
 
 <h3> Summary of the algorithm and testcase </h3>
 
-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 <code>run()</code> function is different and works as
+follows:
 <ol>
 <li>
   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.
 </li>
 </ol>
+
+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.
index b85402b9aec1f30aa97e8f926a4d76104214eefc..8cdb8f3b4c2c22a241a7e311d8887d8c0c431629 100644 (file)
@@ -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:
 
-
 <TABLE WIDTH="100%">
 <tr>
 <td>
@@ -103,6 +92,93 @@ and mesh are shown here:
 <a name="extensions"></a>
 <h3>Possibilities for extensions</h3>
 
+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.
+
+
 <h4> Step length control </h4>
 
-<h4> Integrating mesh refinement and nonlinear solvers </h4>
+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 <i>step length</i> 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 $0<r<1$. It is obvious where the term
+"backtracking" comes from: we try a long step, but if that doesn't work we try
+a shorter step, and ever shorter step, etc. The function
+<code>determine_step_length()</code> 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 &mdash; substantially fewer than we need with the current
+version of the program.
+
+
+<h4> Integrating mesh refinement and nonlinear and linear solvers </h4>
+
+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.
index 4e056430e0133073ca0c8ce7370d7b12339a4900..70b5c4b8bb72523818ce86f4373daa4d23e66024 100644 (file)
@@ -1 +1 @@
-1d problems. A nonlinear problem.
+A nonlinear elliptic problem. Newton's method.
index 41d36fa8b86667ed177a011a41d8e1b3b31c8799..a3375ddf0e85939d9aac853bd526bd3037acaead 100644 (file)
 #include <fstream>
 #include <iostream>
 
-                                 // 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 <deal.II/numerics/solution_transfer.h>
 
-                                 // 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 <code>MinimalSurfaceProblem</code> 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 <code>setup_system</code> 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
+                                  //   <code>refine_mesh</code>.
+                                  // - We then also need new functions:
+                                  //   <code>set_boundary_values()</code>
+                                  //   takes care of setting the boundary
+                                  //   values on the solution vector
+                                  //   correctly, as discussed at the end of
+                                  //   the
+                                  //   introduction. <code>compute_residual()</code>
+                                  //   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,
+                                  //   <code>determine_step_length()</code>
+                                  //   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 <int dim>
   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<double>       present_solution;
       Vector<double>       newton_update;
       Vector<double>       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 <int dim>
   class BoundaryValues : public Function<dim>
@@ -146,6 +165,7 @@ namespace Step15
                            const unsigned int  component = 0) const;
   };
 
+
   template <int dim>
   double BoundaryValues<dim>::value (const Point<dim> &p,
                                     const unsigned int /*component*/) const
@@ -170,7 +190,6 @@ namespace Step15
 
 
 
-                                  //
   template <int dim>
   MinimalSurfaceProblem<dim>::~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
+                                  // <code>refine_mesh()</code>).
 
   template <int dim>
-  void MinimalSurfaceProblem<dim>::setup_system ()
+  void MinimalSurfaceProblem<dim>::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
-                                    // <code>refine_grid</code>
-                                    // 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
-                                        // <code>refine_grid</code>
-                                        // 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 <int dim>
   void MinimalSurfaceProblem<dim>::assemble_system ()
   {
@@ -265,16 +277,19 @@ namespace Step15
     system_rhs = 0;
 
     FEValues<dim> 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<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-    Vector<double>       cell_rhs (dofs_per_cell);
+    FullMatrix<double>           cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>               cell_rhs (dofs_per_cell);
 
-    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    std::vector<Tensor<1, dim> > old_solution_gradients(n_q_points);
+
+    std::vector<unsigned int>    local_dof_indices (dofs_per_cell);
 
     typename DoFHandler<dim>::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<Tensor<1, dim> > 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<dofs_per_cell; ++i)
              {
-               for (unsigned int j = 0; j < dofs_per_cell; ++j)
+               for (unsigned int j=0; j<dofs_per_cell; ++j)
                  {
                    cell_matrix(i, j) += (fe_values.shape_grad(i, q_point)
                                          * coeff
                                          * (fe_values.shape_grad(j, q_point)
-                                            - coeff * coeff
+                                            -
+                                            coeff * coeff
                                             * (fe_values.shape_grad(j, q_point)
-                                               * gradients[q_point])
-                                            * gradients[q_point])
+                                               *
+                                               old_solution_gradients[q_point])
+                                            * old_solution_gradients[q_point]
+                                         )
                                          * fe_values.JxW(q_point));
                  }
 
-               cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff
-                               * gradients[q_point] * fe_values.JxW(q_point));
+               cell_rhs(i) -= (fe_values.shape_grad(i, q_point)
+                               * coeff
+                               * old_solution_gradients[q_point]
+                               * fe_values.JxW(q_point));
              }
          }
 
@@ -351,16 +377,20 @@ namespace Step15
            system_rhs(local_dof_indices[i]) += cell_rhs(i);
          }
       }
+
+                                    // Finally, we remove hanging nodes from
+                                    // the system and apply zero boundary
+                                    // values to the linear system that
+                                    // defines the Newton updates $\delta
+                                    // u^n$:
     hanging_node_constraints.condense (system_matrix);
     hanging_node_constraints.condense (system_rhs);
-    std::map<unsigned int,double> boundary_values;
-
 
+    std::map<unsigned int,double> boundary_values;
     VectorTools::interpolate_boundary_values (dof_handler,
                                              0,
                                              ZeroFunction<dim>(),
                                              boundary_values);
-
     MatrixTools::apply_boundary_values (boundary_values,
                                        system_matrix,
                                        newton_update,
@@ -368,103 +398,14 @@ namespace Step15
   }
 
 
-  template <int dim>
-  double MinimalSurfaceProblem<dim>::compute_residual (const double alpha) const
-  {
-    const QGauss<dim>  quadrature_formula(3);
-
-    Vector<double> residual (dof_handler.n_dofs());
-
-    Vector<double> linearization_point (dof_handler.n_dofs());
-    linearization_point = present_solution;
-    linearization_point.add (alpha, newton_update);
-
-    FEValues<dim> 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<double>       cell_rhs (dofs_per_cell);
-
-    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
-    typename DoFHandler<dim>::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<Tensor<1, dim> > 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<dofs_per_cell; ++i)
-         residual(local_dof_indices[i]) += cell_rhs(i);
-      }
-    hanging_node_constraints.condense (residual);
-
-    std::map<unsigned int,double> boundary_values;
-    VectorTools::interpolate_boundary_values (dof_handler,
-                                             0,
-                                             ZeroFunction<dim>(),
-                                             boundary_values);
-    for (std::map<unsigned int,double>::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 <int dim>
   void MinimalSurfaceProblem<dim>::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 <int dim>
-  double MinimalSurfaceProblem<dim>::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 <int dim>
-  void MinimalSurfaceProblem<dim>::refine_grid ()
+  void MinimalSurfaceProblem<dim>::refine_mesh ()
   {
     Vector<float> 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<dim> 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<double> 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
-                                    // <code>setup_system</code>
-                                    // function in the
-                                    // <code>run</code> 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
+                                    // <code>setup_dofs()</code> 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 <int dim>
   void MinimalSurfaceProblem<dim>::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<unsigned int, double> boundary_values2;
-    VectorTools::interpolate_boundary_values(dof_handler, 0,
-                                            BoundaryValues<dim>(), boundary_values2);
-    for (std::map<unsigned int, double>::const_iterator p =
-          boundary_values2.begin(); p != boundary_values2.end(); ++p)
+    std::map<unsigned int, double> boundary_values;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             BoundaryValues<dim>(),
+                                             boundary_values);
+    for (std::map<unsigned int, double>::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<F(u^n),\varphi_i\right>$ 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<F(u^n+\alpha^n\;\delta u^n),\varphi_i\right>$
+                                  // 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 <int dim>
-  void MinimalSurfaceProblem<dim>::run ()
+  double MinimalSurfaceProblem<dim>::compute_residual (const double alpha) const
   {
+    Vector<double> 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<double> evaluation_point (dof_handler.n_dofs());
+    evaluation_point = present_solution;
+    evaluation_point.add (alpha, newton_update);
 
-    GridGenerator::hyper_ball (triangulation);
-    static const HyperBallBoundary<dim> boundary;
-    triangulation.set_boundary (0, boundary);
-    triangulation.refine_global(2);
+    const QGauss<dim>  quadrature_formula(3);
+    FEValues<dim> 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<double>               cell_rhs (dofs_per_cell);
+    std::vector<Tensor<1, dim> > gradients(n_q_points);
+
+    std::vector<unsigned int>    local_dof_indices (dofs_per_cell);
+
+    typename DoFHandler<dim>::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
+                                        // <code>assemble_system()</code>. 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_point<n_q_points; ++q_point)
          {
-           refine_grid();
+           const double coeff = 1/std::sqrt(1 +
+                                            gradients[q_point] *
+                                            gradients[q_point]);
 
-           std::cout<<"********mesh-refinement:"<<refinement+1<<" ********"<<std::endl;
-           refinement++;
+           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<dofs_per_cell; ++i)
+         residual(local_dof_indices[i]) += cell_rhs(i);
+      }
+
+                                    // At the end of this function we also
+                                    // have to deal with the hanging node
+                                    // constraints and with the issue of
+                                    // boundary values. With regard to the
+                                    // latter, we have to set to zero the
+                                    // elements of the residual vector for
+                                    // all entries that correspond to degrees
+                                    // of freedom that sit at the
+                                    // boundary. The reason is that because
+                                    // the value of the solution there is
+                                    // fixed, they are of course no "real"
+                                    // degrees of freedom and so, strictly
+                                    // speaking, we shouldn't have assembled
+                                    // entries in the residual vector for
+                                    // them. However, as we always do, we
+                                    // want to do exactly the same thing on
+                                    // every cell and so we didn't not want
+                                    // to deal with the question of whether a
+                                    // particular degree of freedom sits at
+                                    // the boundary in the integration
+                                    // above. Rather, we will simply set to
+                                    // zero these entries after the fact. To
+                                    // this end, we first need to determine
+                                    // which degrees of freedom do in fact
+                                    // belong to the boundary and then loop
+                                    // over all of those and set the residual
+                                    // entry to zero. This happens in the
+                                    // following lines:
+    hanging_node_constraints.condense (residual);
+
+    std::vector<bool> boundary_dofs (dof_handler.n_dofs());
+    DoFTools::extract_boundary_dofs (dof_handler,
+                                    std::vector<bool>(1,true),
+                                    boundary_dofs);
+    for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+      if (boundary_dofs[i] == true)
+       residual(i) = 0;
+
+                                    // At the end of the function, we return
+                                    // the norm of the residual:
+    return residual.l2_norm();
+  }
+
 
-                                        // First thing to do after
-                                        // refining the mesh, is to
-                                        // setup the vectors,
-                                        // matrices, etc., which is
-                                        // done in the
-                                        // <code>setup_system</code>
-                                        // 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 <int dim>
+  double MinimalSurfaceProblem<dim>::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 <int dim>
+  void MinimalSurfaceProblem<dim>::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<dim> 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:"<<compute_residual(0)<<std::endl;
+                                    // The Newton iteration starts
+                                    // next. During the first step we do not
+                                    // have information about the residual
+                                    // prior to this step and so we continue
+                                    // the Newton iteration until we have
+                                    // reached at least one iteration and
+                                    // until residual is less than $10^{-3}$.
+                                    //
+                                    // At the beginning of the loop, we do a
+                                    // bit of setup work. In the first go
+                                    // around, we compute the solution on the
+                                    // twice globally refined mesh after
+                                    // setting up the basic data
+                                    // structures. In all following mesh
+                                    // refinement loops, the mesh will be
+                                    // refined adaptively.
+    double previous_res = 0;
+    while (first_step || (previous_res>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:"<<compute_residual(0)<<std::endl;
-         }
 
-                                        // The fifth solution, as
-                                        // well as the Newton update,
-                                        // on every mesh will be
-                                        // written in a vtk-file, in
-                                        // order to show the
-                                        // convergence of the
-                                        // solution.
-
-       Assert (refinement < 100, ExcNotImplemented());
+           first_step = false;
+           std::cout << "  Residual: "
+                     << compute_residual(0)
+                     << std::endl;
+         }
 
+                                        // Every fifth iteration, i.e., just
+                                        // before we refine the mesh again,
+                                        // we output the solution as well as
+                                        // the Newton update. This happens as
+                                        // in all programs before:
        DataOut<dim> 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 ()
index 48a0f738761717ac5d7b461e0b04951a515849bc..78300ce0a286a73c3d8761491d233e1c5e803437 100644 (file)
@@ -1 +1 @@
-step-4
+step-15

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.