From 97a4b05ada3303b1c577fef4568c087021d43848 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 14 Feb 2014 01:44:45 +0000 Subject: [PATCH] improve step-7: switch to vtk; less smooth solution; finer coarse mesh; more cycles when adaptive refining; possible extensions git-svn-id: https://svn.dealii.org/trunk@32483 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-7/doc/intro.dox | 2 +- deal.II/examples/step-7/doc/results.dox | 293 ++++++++++++------------ deal.II/examples/step-7/step-7.cc | 45 ++-- 3 files changed, 174 insertions(+), 166 deletions(-) diff --git a/deal.II/examples/step-7/doc/intro.dox b/deal.II/examples/step-7/doc/intro.dox index 2beb63deee..26dcd7634d 100644 --- a/deal.II/examples/step-7/doc/intro.dox +++ b/deal.II/examples/step-7/doc/intro.dox @@ -110,7 +110,7 @@ where the centers $x_i$ of the exponentials are $x_1=(-\frac 12,\frac 12)$, $x_2=(-\frac 12,-\frac 12)$, and $x_3=(\frac 12,-\frac 12)$. -The half width is set to $\sigma=\frac 13$. +The half width is set to $\sigma=\frac {1}{8}$. We further choose $\Gamma_1=\Gamma \cap\{\{x=1\} \cup \{y=1\}\}$, and there set $g_1$ such that it resembles the exact values of $u$. Likewise, we choose diff --git a/deal.II/examples/step-7/doc/results.dox b/deal.II/examples/step-7/doc/results.dox index b881ada781..1ba7d38af0 100644 --- a/deal.II/examples/step-7/doc/results.dox +++ b/deal.II/examples/step-7/doc/results.dox @@ -2,9 +2,9 @@ The program generates two kinds of output. The first are the output -files solution-adaptive-q1.gmv, -solution-global-q1.gmv, and -solution-global-q2.gmv. We show the latter in a 3d view +files solution-adaptive-q1.vtk, +solution-global-q1.vtk, and +solution-global-q2.vtk. We show the latter in a 3d view here: @@ -19,126 +19,155 @@ screen while running: @code examples/\step-7> make run -============================ Running \step-7 Solving with Q1 elements, adaptive refinement ============================================= Cycle 0: - Number of active cells: 4 - Number of degrees of freedom: 9 + Number of active cells: 64 + Number of degrees of freedom: 81 Cycle 1: - Number of active cells: 13 - Number of degrees of freedom: 22 + Number of active cells: 124 + Number of degrees of freedom: 157 Cycle 2: - Number of active cells: 31 - Number of degrees of freedom: 46 + Number of active cells: 280 + Number of degrees of freedom: 341 Cycle 3: - Number of active cells: 64 - Number of degrees of freedom: 87 + Number of active cells: 577 + Number of degrees of freedom: 690 Cycle 4: - Number of active cells: 127 - Number of degrees of freedom: 160 + Number of active cells: 1099 + Number of degrees of freedom: 1264 Cycle 5: - Number of active cells: 244 - Number of degrees of freedom: 297 + Number of active cells: 2191 + Number of degrees of freedom: 2452 Cycle 6: - Number of active cells: 466 - Number of degrees of freedom: 543 - -cycle cells dofs L2 H1 Linfty - 0 4 9 1.198e+00 2.732e+00 1.383e+00 - 1 13 22 8.795e-02 1.193e+00 1.816e-01 - 2 31 46 8.147e-02 1.167e+00 1.654e-01 - 3 64 87 7.702e-02 1.077e+00 1.310e-01 - 4 127 160 4.643e-02 7.988e-01 6.745e-02 - 5 244 297 2.470e-02 5.568e-01 3.668e-02 - 6 466 543 1.622e-02 4.107e-01 2.966e-02 + Number of active cells: 4165 + Number of degrees of freedom: 4510 +Cycle 7: + Number of active cells: 7915 + Number of degrees of freedom: 8440 +Cycle 8: + Number of active cells: 15196 + Number of degrees of freedom: 15912 + +cycle cells dofs L2 H1 Linfty + 0 64 81 1.576e-01 1.418e+00 2.707e-01 + 1 124 157 4.285e-02 1.285e+00 1.469e-01 + 2 280 341 1.593e-02 7.909e-01 8.034e-02 + 3 577 690 9.359e-03 5.096e-01 2.784e-02 + 4 1099 1264 2.865e-03 3.038e-01 9.822e-03 + 5 2191 2452 1.480e-03 2.106e-01 5.679e-03 + 6 4165 4510 6.907e-04 1.462e-01 2.338e-03 + 7 7915 8440 4.743e-04 1.055e-01 1.442e-03 + 8 15196 15912 1.920e-04 7.468e-02 7.259e-04 Solving with Q1 elements, global refinement =========================================== Cycle 0: - Number of active cells: 4 - Number of degrees of freedom: 9 -Cycle 1: - Number of active cells: 16 - Number of degrees of freedom: 25 -Cycle 2: Number of active cells: 64 Number of degrees of freedom: 81 -Cycle 3: +Cycle 1: Number of active cells: 256 Number of degrees of freedom: 289 -Cycle 4: +Cycle 2: Number of active cells: 1024 Number of degrees of freedom: 1089 -Cycle 5: +Cycle 3: Number of active cells: 4096 Number of degrees of freedom: 4225 -Cycle 6: +Cycle 4: Number of active cells: 16384 Number of degrees of freedom: 16641 -cycle cells dofs L2 H1 Linfty - 0 4 9 1.198e+00 2.732e+00 1.383e+00 - 1 16 25 8.281e-02 1.190e+00 1.808e-01 - 2 64 81 8.142e-02 1.129e+00 1.294e-01 - 3 256 289 2.113e-02 5.828e-01 4.917e-02 - 4 1024 1089 5.319e-03 2.934e-01 1.359e-02 - 5 4096 4225 1.332e-03 1.469e-01 3.482e-03 - 6 16384 16641 3.332e-04 7.350e-02 8.758e-04 - -n cells H1 L2 - 0 4 2.732e+00 - 1.198e+00 - - - 1 16 1.190e+00 1.20 8.281e-02 14.47 3.86 - 2 64 1.129e+00 0.08 8.142e-02 1.02 0.02 - 3 256 5.828e-01 0.95 2.113e-02 3.85 1.95 - 4 1024 2.934e-01 0.99 5.319e-03 3.97 1.99 - 5 4096 1.469e-01 1.00 1.332e-03 3.99 2.00 - 6 16384 7.350e-02 1.00 3.332e-04 4.00 2.00 +cycle cells dofs L2 H1 Linfty + 0 64 81 1.576e-01 1.418e+00 2.707e-01 + 1 256 289 4.280e-02 1.285e+00 1.444e-01 + 2 1024 1089 1.352e-02 7.556e-01 7.772e-02 + 3 4096 4225 3.423e-03 3.822e-01 2.332e-02 + 4 16384 16641 8.586e-04 1.917e-01 6.097e-03 + +n cells H1 L2 +0 64 1.418e+00 - - 1.576e-01 - - +1 256 1.285e+00 1.10 0.14 4.280e-02 3.68 1.88 +2 1024 7.556e-01 1.70 0.77 1.352e-02 3.17 1.66 +3 4096 3.822e-01 1.98 0.98 3.423e-03 3.95 1.98 +4 16384 1.917e-01 1.99 1.00 8.586e-04 3.99 2.00 Solving with Q2 elements, global refinement =========================================== Cycle 0: - Number of active cells: 4 - Number of degrees of freedom: 25 -Cycle 1: - Number of active cells: 16 - Number of degrees of freedom: 81 -Cycle 2: Number of active cells: 64 Number of degrees of freedom: 289 -Cycle 3: +Cycle 1: Number of active cells: 256 Number of degrees of freedom: 1089 -Cycle 4: +Cycle 2: Number of active cells: 1024 Number of degrees of freedom: 4225 -Cycle 5: +Cycle 3: Number of active cells: 4096 Number of degrees of freedom: 16641 -Cycle 6: +Cycle 4: Number of active cells: 16384 Number of degrees of freedom: 66049 -cycle cells dofs L2 H1 Linfty - 0 4 25 1.433e+00 2.445e+00 1.286e+00 - 1 16 81 7.912e-02 1.168e+00 1.728e-01 - 2 64 289 7.755e-03 2.511e-01 1.991e-02 - 3 256 1089 9.969e-04 6.235e-02 2.764e-03 - 4 1024 4225 1.265e-04 1.571e-02 3.527e-04 - 5 4096 16641 1.587e-05 3.937e-03 4.343e-05 - 6 16384 66049 1.986e-06 9.847e-04 5.402e-06 - -n cells H1 L2 - 0 4 2.445e+00 - 1.433e+00 - - - 1 16 1.168e+00 1.07 7.912e-02 18.11 4.18 - 2 64 2.511e-01 2.22 7.755e-03 10.20 3.35 - 3 256 6.235e-02 2.01 9.969e-04 7.78 2.96 - 4 1024 1.571e-02 1.99 1.265e-04 7.88 2.98 - 5 4096 3.937e-03 2.00 1.587e-05 7.97 2.99 - 6 16384 9.847e-04 2.00 1.986e-06 7.99 3.00 +cycle cells dofs L2 H1 Linfty + 0 64 289 1.606e-01 1.278e+00 3.029e-01 + 1 256 1089 7.638e-03 5.248e-01 4.816e-02 + 2 1024 4225 8.601e-04 1.086e-01 4.827e-03 + 3 4096 16641 1.107e-04 2.756e-02 7.802e-04 + 4 16384 66049 1.393e-05 6.915e-03 9.971e-05 + +n cells H1 L2 +0 64 1.278e+00 - - 1.606e-01 - - +1 256 5.248e-01 2.43 1.28 7.638e-03 21.03 4.39 +2 1024 1.086e-01 4.83 2.27 8.601e-04 8.88 3.15 +3 4096 2.756e-02 3.94 1.98 1.107e-04 7.77 2.96 +4 16384 6.915e-03 3.99 1.99 1.393e-05 7.94 2.99 + +Solving with Q2 elements, adaptive refinement +=========================================== + +Cycle 0: + Number of active cells: 64 + Number of degrees of freedom: 289 +Cycle 1: + Number of active cells: 124 + Number of degrees of freedom: 577 +Cycle 2: + Number of active cells: 289 + Number of degrees of freedom: 1353 +Cycle 3: + Number of active cells: 547 + Number of degrees of freedom: 2531 +Cycle 4: + Number of active cells: 1057 + Number of degrees of freedom: 4919 +Cycle 5: + Number of active cells: 2059 + Number of degrees of freedom: 9223 +Cycle 6: + Number of active cells: 3913 + Number of degrees of freedom: 17887 +Cycle 7: + Number of active cells: 7441 + Number of degrees of freedom: 33807 +Cycle 8: + Number of active cells: 14212 + Number of degrees of freedom: 64731 + +cycle cells dofs L2 H1 Linfty + 0 64 289 1.606e-01 1.278e+00 3.029e-01 + 1 124 577 7.891e-03 5.256e-01 4.852e-02 + 2 289 1353 1.070e-03 1.155e-01 4.868e-03 + 3 547 2531 5.962e-04 5.101e-02 1.876e-03 + 4 1057 4919 1.977e-04 3.094e-02 7.923e-04 + 5 2059 9223 7.738e-05 1.974e-02 7.270e-04 + 6 3913 17887 2.925e-05 8.772e-03 1.463e-04 + 7 7441 33807 1.024e-05 4.121e-03 8.567e-05 + 8 14212 64731 3.761e-06 2.108e-03 2.167e-05 @endcode @@ -151,73 +180,37 @@ are the quadratic and cubic rates in the $L_2$ norm. -Finally, the program generated various LaTeX tables. We show here -the convergence table of the Q2 element with global refinement, after -converting the format to HTML: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-n cellsH1-errorL2-error
042.445e+00-1.433e+00--
1161.168e+001.077.912e-0218.114.18
2642.511e-012.227.755e-0310.203.35
32566.235e-022.019.969e-047.782.96
410241.571e-021.991.265e-047.882.98
540963.937e-032.001.587e-057.972.99
6163849.847e-042.001.986e-067.993.00
+Finally, the program also generated LaTeX versions of the tables (not shown +here). + + +

Possible extensions

+ +

Higher Order Elements

+ +Go ahead and run the program with higher order elements (Q3, Q4, ...). You +will notice that assertions in several parts of the code will trigger (for +example in the generation of the filename for the data output). After fixing +these you will not see the correct convergence orders that the theory +predicts. This is because the orders for the quadrature formulas are +hard-coded in this program and this order is not enough for higher order +discretizations. What is a good way to pick the orders dynamically? + +

Convergence Comparison

+ +Is Q1 or Q2 better? What about adaptive versus global refinement? A (somewhat +unfair but typical) metric to compare them, is to look at the error as a +function of the number of unknowns. + +To see this, create a plot in log-log style with the number of unknowns on the +x axis and the L2 error on the y axis. You can add reference lines for +$h^2=N^{-1}$ and $h^3=N^{-3/2}$ and check that global and adaptive refinement +follow those. + +Note that changing the half width of the peaks influences if adaptive or +global refinement is more efficient (if the solution is very smooth, local +refinement does not give any advantage over global refinement). Verify this. +Finally, a more fair comparison would be to plot runtime (switch to release +mode first!) instead of number of unknowns on the x axis. Picking a better +linear solver might be appropriate though. diff --git a/deal.II/examples/step-7/step-7.cc b/deal.II/examples/step-7/step-7.cc index dee41ef7e9..0f12d98b4b 100644 --- a/deal.II/examples/step-7/step-7.cc +++ b/deal.II/examples/step-7/step-7.cc @@ -152,7 +152,7 @@ namespace Step7 // concrete instantiation by substituting dim with a concrete // value: template - const double SolutionBase::width = 1./3.; + const double SolutionBase::width = 1./8.; @@ -976,12 +976,13 @@ namespace Step7 template void HelmholtzProblem::run () { - for (unsigned int cycle=0; cycle<7; ++cycle) + const unsigned int n_cycles = (refinement_mode==global_refinement)?5:9; + for (unsigned int cycle=0; cycle::cell_iterator cell = triangulation.begin (), @@ -1019,21 +1020,21 @@ namespace Step7 // After the last iteration we output the solution on the finest // grid. This is done using the following sequence of statements which we // have already discussed in previous examples. The first step is to - // generate a suitable filename (called gmv_filename here, - // since we want to output data in GMV format; we add the prefix to + // generate a suitable filename (called vtk_filename here, + // since we want to output data in VTK format; we add the prefix to // distinguish the filename from that used for other output files further // down below). Here, we augment the name by the mesh refinement // algorithm, and as above we make sure that we abort the program if // another refinement method is added and not handled by the following // switch statement: - std::string gmv_filename; + std::string vtk_filename; switch (refinement_mode) { case global_refinement: - gmv_filename = "solution-global"; + vtk_filename = "solution-global"; break; case adaptive_refinement: - gmv_filename = "solution-adaptive"; + vtk_filename = "solution-adaptive"; break; default: Assert (false, ExcNotImplemented()); @@ -1053,10 +1054,10 @@ namespace Step7 switch (fe->degree) { case 1: - gmv_filename += "-q1"; + vtk_filename += "-q1"; break; case 2: - gmv_filename += "-q2"; + vtk_filename += "-q2"; break; default: @@ -1064,10 +1065,10 @@ namespace Step7 } // Once we have the base name for the output file, we add an extension - // appropriate for GMV output, open a file, and add the solution vector to + // appropriate for VTK output, open a file, and add the solution vector to // the object that will do the actual output: - gmv_filename += ".gmv"; - std::ofstream output (gmv_filename.c_str()); + vtk_filename += ".vtk"; + std::ofstream output (vtk_filename.c_str()); DataOut data_out; data_out.attach_dof_handler (dof_handler); @@ -1100,9 +1101,9 @@ namespace Step7 // way as above. // // With the intermediate format so generated, we can then actually write - // the graphical output in GMV format: + // the graphical output: data_out.build_patches (fe->degree); - data_out.write_gmv (output); + data_out.write_vtk (output); // @sect5{Output of convergence tables} @@ -1231,6 +1232,8 @@ namespace Step7 convergence_table .evaluate_convergence_rates("L2", ConvergenceTable::reduction_rate_log2); convergence_table + .evaluate_convergence_rates("H1", ConvergenceTable::reduction_rate); + convergence_table .evaluate_convergence_rates("H1", ConvergenceTable::reduction_rate_log2); // Each of these function calls produces an additional column that is // merged with the original column (in our example the `L2' and the @@ -1347,7 +1350,19 @@ int main () std::cout << std::endl; } + { + std::cout << "Solving with Q2 elements, adaptive refinement" << std::endl + << "===========================================" << std::endl + << std::endl; + + FE_Q fe(2); + HelmholtzProblem + helmholtz_problem_2d (fe, HelmholtzProblem::adaptive_refinement); + + helmholtz_problem_2d.run (); + std::cout << std::endl; + } } catch (std::exception &exc) { -- 2.39.5