From: Wolfgang Bangerth Date: Thu, 8 Jun 2023 19:51:16 +0000 (-0600) Subject: Add something about postprocessing to step-4. X-Git-Tag: v9.5.0-rc1~123^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1996bd633797fd6ef21aa1769f6abe40ac3f64e7;p=dealii.git Add something about postprocessing to step-4. --- diff --git a/examples/step-4/doc/results.dox b/examples/step-4/doc/results.dox index 70e2ac1520..2d66673833 100644 --- a/examples/step-4/doc/results.dox +++ b/examples/step-4/doc/results.dox @@ -100,11 +100,323 @@ wants to see, but no more. + +

Postprocessing: What to do with the solution?

+ +This tutorial -- like most of the other programs -- principally only shows how +to numerically approximate the solution of a partial differential equation, and +then how to visualize this solution graphically. But +solving a PDE is of course not the goal in most practical applications (unless +you are a numerical methods developer and the *method* is the goal): We generally +want to solve a PDE because we want to *extract information* from it. Examples +for what people are interested in from solutions include the following: +- Let's say you solve the equations of elasticity (which we will do in step-8), + then that's presumably because you want to know about the deformation of an + elastic object under a given load. From an engineering perspective, what you + then presumably want to learn is the degree of deformation of the object, + say at a specific point; or you may want to know the maximum + [stress](https://en.wikipedia.org/wiki/Stress_(mechanics)) in order to + determine whether the applied load exceeds the safe maximal stress the + material can withstand. +- If you are solving fluid flow problems (such as in step-22, step-57, step-67, + and step-69), then you might be interested in the fluid velocity at specific + points, and oftentimes the forces the fluid exerts on the boundary of the + fluid domain. The latter is important in many applications: If the fluid + in question is the air flowing around an airplane, then we are typically + interested in the drag and lift forces on the fuselage and wings. If the + fluid is water flowing around a ship, then we typically care about + the drag force on the ship. +- If you are solving the Maxwell equations of electromagnetics, you are + typically interested in how much energy is radiated in certain directions + (say, in order to know the range of information transmission via an + antenna, or to determine the + [radar cross section](https://en.wikipedia.org/wiki/Radar_cross_section) + of planes or ships). + +The point here is that from an engineering perspective, *solving* the +PDE is only the first step. The second step is to *evaluate* the computed +solution in order to extract relevant numbers that allow us to +either *optimize a design*, or to *make decisions*. This second step is often +called "postprocessing the solution". + +This program does not solve a solid or fluid mechanics problem, so we +should try to illustrate postprocessing with something that makes sense +in the context of the equation we solve here. The Poisson equation in two +space dimensions is a model for the vertical deformation of a membrane +that is clamped at the boundary and is subject to a vertical force. For this +kind of situation, it makes sense to evaluate the *average vertical +displacement*, +@f[ + \bar u_h = \frac{\int_\Omega u_h(\mathbf x) \, dx}{|\Omega|}, +@f] +where $|\Omega| = \int_\Omega 1 \, dx$ is the area of the domain. To compute +$\bar u_h$, as usual we replace integrals over the domain by a sum of integrals +over cells, +@f[ + \int_\Omega u_h(\mathbf x) \, dx + = + \sum_K \int_K u_h(\mathbf x) \, dx, +@f] +and then integrals over cells are approximated by quadrature: +@f{align*}{ + \int_\Omega u_h(\mathbf x) \, dx + &= + \sum_K \int_K u_h(\mathbf x) \, dx, + \\ + &= + \sum_K \sum_q u_h(\mathbf x_q^K) w_q^K, +@f} +where $w_q^K$ is the weight of the $q$th quadrature point evaluated on +cell $K$. All of this is as always provided by the FEValues class -- the +entry point for all integrals in deal.II. + +The actual implementation of this is straightforward once you know how +to get the values of the solution $u$ at the quadrature points of a cell. +This functionality is provided by FEValues::get_function_values(), a +function that takes a global vector of nodal values as input and returns +a vector of function values at the quadrature points of the current cell. +Using this function, to see how it all works together you can +place the following code snippet anywhere in the program after the +solution has been computed (the `output_results()` function seems like a good +place to also do postprocessing, for example): +@code + QGauss quadrature_formula(fe.degree + 1); + FEValues fe_values(fe, + quadrature_formula, + update_values | update_JxW_values); + + std::vector solution_values(quadrature_formula.size()); + double integral_of_u = 0; + double volume_of_omega = 0; + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + fe_values.get_function_values(solution, solution_values); + + for (const unsigned int q_point : fe_values.quadrature_point_indices()) + { + integral_of_u += solution_values[q_point] * fe_values.JxW(q_point); + volume_of_omega += 1 * fe_values.JxW(q_point); + } + } + std::cout << " Mean value of u=" << integral_of_u / volume_of_omega + << std::endl; +@endcode +In this code snippet, we also compute the volume (or, since we are currently +thinking about a two-dimensional situation: the area) $|\Omega|$ by computing +the integral $|\Omega| = \int_\Omega 1 \, dx$ in exactly the same way, via +quadrature; this could be avoided by using the fact that deal.II has a function +for this, GridTools::volume(), but it is efficient to compute the two integrals +at the same time, and so that's what we do. + +This program of course also solves the same Poisson equation in three space +dimensions. In this situation, the Poisson equation is often used as a model +for diffusion of either a physical species (say, of ink in a tank of water, +or a pollutant in the air) or of energy (specifically, of thermal energy in +a solid body). In that context, the quantity +@f[ + \Phi_h = \int_{\partial\Omega} \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx +@f] +is the *flux* of this species or energy across the boundary. (In actual +physical models, one would also have to multiply the right hand side by +a diffusivity or conductivity constant, but let us ignore this here.) In +much the same way as before, we compute such integrals by splitting +it over integrals of *faces* of cells, and then applying quadrature: +@f{align*}{ + \Phi_h + &= + \int_{\partial\Omega} \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx + \\ + &= + \sum_K + \sum_{f \in \text{faces of $K$}, f\subset\partial\Omega} + \int_f \nabla u_h(\mathbf x) \cdot \mathbf n(\mathbf x) \; dx + \\ + &= + \sum_K + \sum_{f \in \text{faces of $K$}, f\subset\partial\Omega} + \sum_q \nabla u_h(\mathbf x_q^f) \cdot \mathbf n(\mathbf x_q^f) w_q^f, +@f} +where now $\mathbf x_q^f$ are the quadrature points located on face $f$, +and $w_q^f$ are the weights associated with these faces. The second +of the sum symbols loops over all faces of cell $K$, but restricted to +those that are actually at the boundary. + +This all is easily implemented by the following code that replaces the use of the +FEValues class (which is used for integrating over cells -- i.e., domain integrals) +by the FEFaceValues class (which is used for integrating over faces -- i.e., +boundary integrals): +@code + QGauss face_quadrature_formula(fe.degree + 1); + FEFaceValues fe_face_values(fe, + face_quadrature_formula, + update_gradients | update_normal_vectors | + update_JxW_values); + + std::vector> solution_gradients(face_quadrature_formula.size()); + double flux = 0; + + for (const auto &cell : dof_handler.active_cell_iterators()) + for (const auto &face : cell->face_iterators()) + if (face->at_boundary()) + { + fe_face_values.reinit(cell, face); + fe_face_values.get_function_gradients(solution, solution_gradients); + + for (const unsigned int q_point : + fe_face_values.quadrature_point_indices()) + { + flux += solution_gradients[q_point] * + fe_face_values.normal_vector(q_point) * + fe_face_values.JxW(q_point); + } + } + std::cout << " Flux=" << flux << std::endl; +@endcode + +If you add these two code snippets to the code, you will get output like the +following when you run the program: +@code +Solving problem in 2 space dimensions. + Number of active cells: 256 + Total number of cells: 341 + Number of degrees of freedom: 289 + 26 CG iterations needed to obtain convergence. + Mean value of u=1.33303 + Flux=-3.68956 +Solving problem in 3 space dimensions. + Number of active cells: 4096 + Total number of cells: 4681 + Number of degrees of freedom: 4913 + 30 CG iterations needed to obtain convergence. + Mean value of u=1.58058 + Flux=-8.29435 +@endcode + +This makes some sense: If you look, for example, at the 2d output above, +the solution varies between values of 1 and 2, but with a larger part of the +solution closer to one than two; so an average value of 1.33 for the mean value +is reasonable. For the flux, recall that $\nabla u \cdot \mathbf n$ is the +directional derivative in the normal direction -- in other words, how the +solution changes as we move from the interior of the domain towards the +boundary. If you look at the 2d solution, you will realize that for most parts +of the boundary, the solution *decreases* as we approach the boundary, so the +normal derivative is negative -- so if we integrate along the boundary, we +should expect (and obtain!) a negative value. + + +

Possibilities for extensions

+There are many ways with which one can play with this program. The simpler +ones include essentially all the possibilities already discussed in the +Possibilities for extensions in the documentation of step 3, +except that you will have to think about whether something now also applies +to the 3d case discussed in the current program. + +It is also worthwhile considering the postprocessing options discussed +above. The documentation states two numbers (the mean value and the +normal flux) for both the 2d and 3d cases. Can we trust these +numbers? We have convinced ourselves that at least the mean value +is reasonable, and that the sign of the flux is probably correct. +But are these numbers accurate? + +A general rule is that we should never trust a number unless we have +verified it in some way. From the theory of finite element methods, +we know that as we make the mesh finer and finer, the numerical +solution $u_h$ we compute here must converge to the exact solution +$u$. As a consequence, we also expect that $\bar u_h \rightarrow \bar u$ +and $\Phi_h \rightarrow \Phi$, but that does not mean that for any +given mesh $\bar u_h$ or $\Phi_h$ are particularly accurate approximations. + +To test this kind of thing, we have already considered the convergence of +a point value in step-3. We can do the same here by selecting how many +times the mesh is globally refined in the `make_grid()` function of this +program. For the mean value of the solution, we then get the following +numbers: + + + + + + + + + + + +
# of refinements$\bar u_h$ in 2d$\bar u_h$ in 3d
4 1.33303 1.58058
5 1.33276 1.57947
6 1.3327 1.5792
7 1.33269 1.57914
8 1.33268
9 1.33268
+I did not have the patience to run the last two values for the 3d case -- +one needs quite a fine mesh for this, with correspondingly long run times. +But we can be reasonably assured that values around 1.33 (for the 2d case) +and 1.58 (for the 3d case) are about right -- and at least for engineering +applications, three digits of accuracy are good enough. + +The situation looks very different for the flux. Here, we get results +such as the following: + + + + + + + + + + + + + +
# of refinements$\Phi_h$ in 2d$\Phi_h$ in 3d
4 -3.68956 -8.29435
5 -4.90147 -13.0691
6 -5.60745 -15.9171
7 -5.99111 -17.4918
8 -6.19196
9 -6.29497
10 -6.34721
11 -6.37353
+So this is not great. For the 2d case, we might infer that perhaps +a value around -6.4 might be right if we just refine the mesh enough -- +though 11 refinements already leads to some 4,194,304 cells. In any +case, the first number (the one shown in the beginning where we +discussed postprocessing) was off by almost a factor of 2! + +For the 3d case, the last number shown was on a mesh with 2,097,152 +cells; the next one would have had 8 times as many cells. In any case, the +numbers mean that we can't even be sure +that the first digit of that last number is correct! In other words, +it was worth checking, or we would have just believed all of these +numbers. In fact, that last column isn't even doing a particularly +good job convincing that the code might be correctly implemented. + +If you keep reading through the tutorial, you will find many ways +to make these sorts of computations more accurate and to come to +believe that the flux actually does converge to a correct value. +For example, we can dramatically increase the accuracy of the computation +by using adaptive mesh refinement (step-6) near the boundary, and +in particular by using higher polynomial degree finite elements (also +step-6, but also step-7). Using the latter, using cubic elements +(polynomial degree 3), we can actually compute the flux pretty +accurately even in 3d: $\Phi_h=-19.0148$ with 4 global refinement steps, +and $\Phi_h=-19.1533$ with 5 refinement steps. These numbers are already +pretty close together and give us a reasonable idea of the first +two correct digits of the "true" answer. -Essentially the possibilities for playing around with the program are the same -as for the previous one, except that they will now also apply to the 3d -case. For inspiration read up on possible extensions in the documentation of step 3. +@note We would be remiss to not also comment on the fact that there + are good theoretical reasons why computing the flux accurately + appears to be so much more difficult than the average value. + This has to do with the fact that finite element theory + provides us with the estimate + $\|u-u_h\|_{L_2(\Omega)} \le C h^2 \|\nabla^2u\|_{L_2(\Omega)}$ + when using the linear elements this program uses -- that is, for + every global mesh refinement, $h$ is reduced by a factor of two + and the error goes down by a factor of 4. Now, the $L_2$ error is + not equivalent to the error in the mean value, but the two are + related: They are both integrals over the domain, using the *value* + of the solution. We expect the mean value to converge no worse than + the $L_2$ norm of the error. At the same time, theory also provides + us with this estimate: + $\|\nabla (u-u_h)\|_{L_2(\partial\Omega)} \le + C h^{1/2} \|\nabla^2u\|_{L_2(\Omega)}$. The move from values to + gradients reduces the convergence rates by one order, and the move + from domain to boundary by another half order. Here, then, each + refinement step reduces the error not by a factor of 4 any more, + by only by a factor of $\sqrt{2} \approx 1.4$. It takes a lot + of global refinement steps to reduce the error by, say, a factor + ten or hundred, and this is reflected in the very slow convergence + evidenced by the table.