Skip to content
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

Merged
merged 3 commits into from Jun 13, 2023
Merged

Conversation

bangerth
Copy link
Member

@bangerth bangerth commented Jun 8, 2023

Fixes #15309.

Copy link
Member

@kronbichler kronbichler left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a good description. I have a few suggestions below. The most important is the one about the available analytical result for the flux in this particular case: I think it makes sense to mention knowledge a mathematically educated reader might have, but at the same time the example is representative of many much more complicated evaluation tasks where no analytical expression is given.

examples/step-4/doc/results.dox Outdated Show resolved Hide resolved
examples/step-4/doc/results.dox Outdated Show resolved Hide resolved
Comment on lines 210 to 211
for this, GridTools::volume(), but it is efficient to compute the two integrals
at the same time, and so that's what we do.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not sure if efficiency is the most important argument here. To me, the educative aspect, showing how the ingredients are computed as a means to actually understand things, is most important, rather than relying on a black-box function. So I would simply remove the last clause starting from "but".

Comment on lines +299 to +307
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.
Copy link
Member

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.

Copy link
Member Author

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!

examples/step-4/doc/results.dox Outdated Show resolved Hide resolved
Comment on lines +351 to +353
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.
Copy link
Member

Choose a reason for hiding this comment

The 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
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.
We did not run the solver for the last two values for the 3d case,
because as the mesh is refined the number of unknowns grows quickly.
Then, the run times would be unacceptably high with the basic iterative
solver chosen here. In fact, it is not optimal for Poisson-type problems;
for better solvers, see, e.g., step-16 or step-37.

Copy link
Member Author

Choose a reason for hiding this comment

The 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.

examples/step-4/doc/results.dox Outdated Show resolved Hide resolved
Comment on lines 385 to 386
numbers. In fact, that last column isn't even doing a particularly
good job convincing that the code might be correctly implemented.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While I entirely agree with the message given here, I think this particular case is special because there is an easy-to-compute reference result as I wrote in my comment above. I guess one could add that, if available, mathematical identities of results of benchmarks from the literature should be consulted. Solving PDEs is not really a new field, and readers should learn that many people have thought about similar cases before.

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.
Copy link
Member

Choose a reason for hiding this comment

The 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 4(x^4+y^4+z^4), which can be integrated analytically from the 1d result int_{-1}^1 4 x^4 dx = 8/5.

examples/step-4/doc/results.dox Outdated Show resolved Hide resolved
@bangerth
Copy link
Member Author

Thanks for the review, @kronbichler! I've merged the minor edits into the previous commits and left the new material to the new third commit. I did not even think about computing the exact value of the flux -- thanks for pointing this out. Can you take another look?

Copy link
Member

@kronbichler kronbichler left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice, I like this discussion, as it covers many facets in finite element programs. The solver topic can be postponed indeed.

@kronbichler kronbichler merged commit be577f0 into dealii:master Jun 13, 2023
11 of 14 checks passed
@bangerth bangerth deleted the 4 branch June 13, 2023 14:10
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Talk about actually doing something with the solution computed in tutorials
2 participants