-
Notifications
You must be signed in to change notification settings - Fork 708
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add something about postprocessing to step-4. #15329
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change | ||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
@@ -100,11 +100,363 @@ wants to see, but no more. | |||||||||||||||
|
||||||||||||||||
|
||||||||||||||||
|
||||||||||||||||
<a name="postprocessing"></a> | ||||||||||||||||
<h3> Postprocessing: What to do with the solution? </h3> | ||||||||||||||||
|
||||||||||||||||
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<dim> quadrature_formula(fe.degree + 1); | ||||||||||||||||
FEValues<dim> fe_values(fe, | ||||||||||||||||
quadrature_formula, | ||||||||||||||||
update_values | update_JxW_values); | ||||||||||||||||
|
||||||||||||||||
std::vector<double> 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. (We could avoid having to compute $|\Omega|$ by hand here, using the | ||||||||||||||||
fact that deal.II has a function for this: GridTools::volume(). That said, | ||||||||||||||||
it is efficient to compute the two integrals | ||||||||||||||||
concurrently in the same loop, 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<dim - 1> face_quadrature_formula(fe.degree + 1); | ||||||||||||||||
FEFaceValues<dim> fe_face_values(fe, | ||||||||||||||||
face_quadrature_formula, | ||||||||||||||||
update_gradients | update_normal_vectors | | ||||||||||||||||
update_JxW_values); | ||||||||||||||||
|
||||||||||||||||
std::vector<Tensor<1, dim>> 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. | ||||||||||||||||
|
||||||||||||||||
|
||||||||||||||||
|
||||||||||||||||
<a name="extensions"></a> | ||||||||||||||||
<h3>Possibilities for extensions</h3> | ||||||||||||||||
|
||||||||||||||||
There are many ways with which one can play with this program. The simpler | ||||||||||||||||
ones include essentially all the possibilities already discussed in the | ||||||||||||||||
<a href="step_3.html#extensions" target="body">Possibilities for extensions in the documentation of step 3</a>, | ||||||||||||||||
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: | ||||||||||||||||
<table align="center" class="doxtable"> | ||||||||||||||||
<tr> <th># of refinements</th> | ||||||||||||||||
<th>$\bar u_h$ in 2d</th> | ||||||||||||||||
<th>$\bar u_h$ in 3d</th> | ||||||||||||||||
</tr> | ||||||||||||||||
<tr> <td>4</td> <td>1.33303</td> <td>1.58058</td> </tr> | ||||||||||||||||
<tr> <td>5</td> <td>1.33276</td> <td>1.57947</td> </tr> | ||||||||||||||||
<tr> <td>6</td> <td>1.3327</td> <td>1.5792</td> </tr> | ||||||||||||||||
<tr> <td>7</td> <td>1.33269</td> <td>1.57914</td> </tr> | ||||||||||||||||
<tr> <td>8</td> <td>1.33268</td> <td></td> </tr> | ||||||||||||||||
<tr> <td>9</td> <td>1.33268</td> <td></td> </tr> | ||||||||||||||||
</table> | ||||||||||||||||
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. | ||||||||||||||||
Comment on lines
+352
to
+353
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think this statement is repeated further down (where I suggest to delete it altogether), and I do not like this. I also think that the problem with 16 million unknowns will hardly fit in memory, whereas the next one with 135m DoFs will certainly not in a typical laptop as it would take around 80GB for just the matrix. I suggest the following:
Suggested change
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I have wanted to add something about solvers elsewhere (https://dealii.org/developer/doxygen/deal.II/step_6.html#Solversandpreconditioners is the right place). So let me skip this suggestion here and instead pick it up in step-6. |
||||||||||||||||
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: | ||||||||||||||||
<table align="center" class="doxtable"> | ||||||||||||||||
<tr> <th># of refinements</th> | ||||||||||||||||
<th>$\Phi_h$ in 2d</th> | ||||||||||||||||
<th>$\Phi_h$ in 3d</th> | ||||||||||||||||
</tr> | ||||||||||||||||
<tr> <td>4</td> <td>-3.68956</td> <td>-8.29435</td> </tr> | ||||||||||||||||
<tr> <td>5</td> <td>-4.90147</td> <td>-13.0691</td> </tr> | ||||||||||||||||
<tr> <td>6</td> <td>-5.60745</td> <td>-15.9171</td> </tr> | ||||||||||||||||
<tr> <td>7</td> <td>-5.99111</td> <td>-17.4918</td> </tr> | ||||||||||||||||
<tr> <td>8</td> <td>-6.19196</td> <td></td> </tr> | ||||||||||||||||
<tr> <td>9</td> <td>-6.29497</td> <td></td> </tr> | ||||||||||||||||
<tr> <td>10</td> <td>-6.34721</td> <td></td> </tr> | ||||||||||||||||
<tr> <td>11</td> <td>-6.37353</td> <td></td> </tr> | ||||||||||||||||
</table> | ||||||||||||||||
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 us that the code might be correctly implemented. | ||||||||||||||||
|
||||||||||||||||
If you keep reading through the other tutorial programs, 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 its 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. | ||||||||||||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This discussion shows convincingly is that one should not use polynomial degree 1 when trying to evaluate the gradient, which is then a piecewise constant (well, not entirely do to bilinear shape functions), because piecewise constants are way to inefficient. I like the fact that cubic polynomials give much better results for derivatives. And in fact, the correct result for this quantity is -6.4 in 2d and -19.2 in 3d (analytically), because the right hand side is |
||||||||||||||||
|
||||||||||||||||
@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. On the other hand, for cubic elements (i.e., | ||||||||||||||||
polynomial degree 3), we would get | ||||||||||||||||
$\|u-u_h\|_{L_2(\Omega)} \le C h^4 \|\nabla^4u\|_{L_2(\Omega)}$ | ||||||||||||||||
and after reduction by 1.5 orders, we would still have | ||||||||||||||||
$\|\nabla (u-u_h)\|_{L_2(\partial\Omega)} \le | ||||||||||||||||
C h^{2+1/2} \|\nabla^4u\|_{L_2(\Omega)}$. This rate, | ||||||||||||||||
${\cal O}(h^{2.5})$ is still quite rapid, and it is perhaps not | ||||||||||||||||
surprising that we get much better answers with these higher | ||||||||||||||||
order elements. This also illustrates that when trying to | ||||||||||||||||
approximate anything that relates to a gradient of the solution, | ||||||||||||||||
using linear elements (polynomial degree one) is really not a | ||||||||||||||||
good choice at all. | ||||||||||||||||
|
||||||||||||||||
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 <a href="step_3.html#extensions" | ||||||||||||||||
target="body">possible extensions in the documentation of step 3</a>. | ||||||||||||||||
@note In this very specific case, it turns out that we can actually | ||||||||||||||||
compute the exact value of $\Phi$. This is because for the Poisson | ||||||||||||||||
equation we compute the solution of here, $-\Delta u = f$, we can | ||||||||||||||||
integrate over the domain, $-\int_\Omega \Delta u = \int_\Omega f$, | ||||||||||||||||
and then use that $\Delta = \text{div}\;\text{grad}$; this allows | ||||||||||||||||
us to use the | ||||||||||||||||
divergence theorem followed by multiplying by minus one to find | ||||||||||||||||
$\int_{\partial\Omega} \nabla u \cdot n = -\int_\Omega f$. The | ||||||||||||||||
left hand side happens to be $\Phi$. For the specific right | ||||||||||||||||
hand side $f(x_1,x_2)=4(x_1^4+x_2^4)$ we use in 2d, we then | ||||||||||||||||
get $-\int_\Omega f = -\int_{-1}^{1} \int_{-1}^{1} 4(x_1^4+x_2^4) \; dx_2\; dx_1 | ||||||||||||||||
= -16 \left[\int_{-1}^{1} x^4 \; dx\right] = -16\times\frac 25$, | ||||||||||||||||
which has a numerical value of exactly -6.4 -- right on with our | ||||||||||||||||
guess above. In 3d, we can do the same and get that the exact | ||||||||||||||||
value is | ||||||||||||||||
$-\int_\Omega f = | ||||||||||||||||
-\int_{-1}^{1} \int_{-1}^{1} \int_{-1}^{1} 4(x_1^4+x_2^4+x_3^4) \; dx_3 \; dx_2\; dx_1 | ||||||||||||||||
= -48\times\frac 25=-19.2$. What we found with cubic elements | ||||||||||||||||
is then quite close to this exact value. Of course, in practice | ||||||||||||||||
we almost never have exact values to compare with: If we could | ||||||||||||||||
compute something on a piece of paper, we wouldn't have to solve | ||||||||||||||||
the PDE numerically. But these sorts of situations make for excellent | ||||||||||||||||
test cases that help us verify that our numerical solver works | ||||||||||||||||
correctly. In many other cases, the literature contains | ||||||||||||||||
numbers where others have already computed an answer accurately | ||||||||||||||||
using their own software, and these are also often useful to | ||||||||||||||||
compare against in verifying the correctness of our codes. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I did not check the equation in detail, but shouldn't we simply get
int_\Omega -right_hand_side(x) dx
for the value of the flux? This is a consequence of the Gauss divergence theorem / integration by parts, as the integral of the flux over the boundary is equal to the integral of the Laplacian in the interior, which in turn is equal to the negative right hand side.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You're right. I did not think of that!