Step-90 TraceFEM tutorial #16878
Step-90 TraceFEM tutorial #16878
Conversation
|
Thanks, Vladimir. Can you upload a PDF of the documentation to help with review? You might be able to add it as a file here by dropping it inside the comment field. |
|
It seems that the build problems are related to how TrilinosWrappers and linear algebra vectors are used in step-90.cc. Should I use the LA namespace from step-40 and actually follow its code more closely to try to fix the issues? |
|
Can you fix the CMakeLists.txt to require Trilinos? See step-32. I think the problems will go away then. |
|
/rebuild |
examples/step-90/doc/intro.dox
Outdated
| This program was contributed by Vladimir Yushutin. | ||
|
|
||
| The material is based upon work partially supported by NSF. |
There was a problem hiding this comment.
Want to state your affiliation and the grant number(s) from NSF?
There was a problem hiding this comment.
@tjhei , could you please provide the appropriate grant number?
examples/step-90/doc/intro.dox
Outdated
| <a name="Intro to step-90"></a> | ||
| <h1>Introduction</h1> | ||
|
|
||
| In this example, we implement the trace finite element method (TraceFEM) in deal.II. |
There was a problem hiding this comment.
Before you jump right into it in the next sentence, I think this would be a good place to describe what TraceFEM actually is, and what kinds of problems it solves. Perhaps also to reference a paper or two.
In fact, before you describe the method itself, give an outline of the problem you want to solve, and what is difficult about it.
There was a problem hiding this comment.
Here is what I've added:
TraceFEM solves PDEs posed on a possibly evolving
Being an unfitted method, TraceFEM allows to circumvent the need of remeshing of an evolving surface if it is implicitly given by the zero contour of a level-set function. At the same time, it easily provides with a natural extension of the discrete solution which must be constructed by any time-stepping technique. Certainly, this flexibility comes with the price: one needs to design the nodes and weights for a quadrature customized for each implicit intersection of the zero level-set and the background mesh. Moreover, these intersections may be of arbitrary shape and size manifesting in the so-called ``small cut" problem and requiring addition of a stabilization form that restores well-conditioning of the problem.
There was a problem hiding this comment.
Yes, I think that's great!
In "provides with a natural extension of the discrete solution" I suppose you mean "extension in the direction perpendicular to the surface"? (If so, say "provides a natural", dropping the "with".)
"comes with the price" -> "comes with a price"
There was a problem hiding this comment.
Ahh, I see... I did mean "an extension". I rewrote it to shift the focus as follows:
At the same time, it easily provides with an extension of the discrete solution to a neighborhood of the surface which turns out to be very handy in case of non-stationary interfaces and films.
examples/step-90/doc/intro.dox
Outdated
| Two aspects are of our focus. First, the surface approximation is separated from the discretization of the surface PDE, | ||
| e.g. a $Q_2$ discrete level-set and a $Q_1$ solution are possible, but both should be based on the same bulk triangulation. | ||
| Second, we make sure that performance of TraceFEM in the parallel implementation corresponds to that of a classical | ||
| fitted FEM for a two-dimensional problem. We demonstrate how to achieve both goals by using a combination of MeshWorker |
There was a problem hiding this comment.
Since you compare with "fitted FEM", perhaps that's something that could go into the description of the method I'm asking for above as well: How does TraceFEM compare to other approaches, such as fitted FEM? What do each of these methods do, and how are they different?
There was a problem hiding this comment.
A natural alternative to TraceFEM in solving surface PDEs is the parametric surface finite element method. The latter
method relies on an explicit parametrization of the surface which may be not feasible especially for evolving interfaces with an unknown in advance shape - in this sense, TraceFEM is a technique inspired by the level-set description of interfaces. However, the parametric surface finite element method, when applicable, enjoys many well-known properties of fitted methods on flat domains provided the geometric errors - which are present for both methods - are taken under control.
examples/step-90/doc/results.dox
Outdated
| | Cycle | DOFS | EOC | Iterations | $L^2$-Error | EOC | $H^1$-Error | EOC |Stabilization| EOC | | ||
| |:-----:|:--------:|:-----:|:----------:|:-----------:|:-----:|:-----------:|:-----:|:-----------:|:-----:| | ||
| | 0 | 12370 | - | 15 | 7.6322e-02 | - | 3.6212e-01 | - | 2.2423e-01 | - | | ||
| | 1 | 49406 | -2.00 | 18 | 1.1950e-02 | 2.68 | 1.4752e-01 | 1.30 | 1.1238e-01 | 1.00 | | ||
| | 2 | 196848 | -1.99 | 19 | 1.7306e-03 | 2.79 | 7.4723e-02 | 0.98 | 6.1131e-02 | 0.88 | | ||
| | 3 | 785351 | -2.00 | 22 | 3.6276e-04 | 2.25 | 3.9329e-02 | 0.93 | 3.0185e-02 | 1.02 | | ||
| | 4 | 3136501 | -2.00 | 25 | 7.5910e-05 | 2.26 | 1.9694e-02 | 1.00 | 1.4875e-02 | 1.02 | | ||
| | 5 | 12536006 | -2.00 | 26 | 1.7279e-05 | 2.14 | 9.8443e-03 | 1.00 | 7.4067e-03 | 1.01 | | ||
| | 6 | 50122218 | -2.00 | 30 | 4.3891e-06 | 1.98 | 4.9219e-03 | 1.00 | 3.7042e-03 | 1.00 | |
There was a problem hiding this comment.
What does EOC stand for? And why would you show EOC on the number of DoFs?
There was a problem hiding this comment.
You are right, this is not a "experimental order of convergence". I'll use "Rate" instead. The message of this rate is that it corresponds to the rate of a two-dimensional problem. For a fitted 3D problem we should see something like -3.00. There was another column titled "LevelSet DOFS" with the same 2D rate -2.00 despite the fact that the mesh is three-dimensional. Would you like to see it here too?
There was a problem hiding this comment.
Ah, I see. Perhaps put that interpretation into text below the table to make it clear what you want the interpretation to be.
examples/step-90/doc/results.dox
Outdated
| In this test we refine the mesh near the surface and, as a result, the number of degrees of freedom scales in the two-dimensional fashion. | ||
| The optimal rates of error convergence in $L^2(\Gamma)$ and $H^1(\Gamma)$ norms are clearly observable. We also note | ||
| the first order convergence of the stabilization term $s_h(u_h, u_h)$. |
There was a problem hiding this comment.
Can you explain here what you show in "Stabilization"? The term as shown in the introduction is a bilinear form, not a single number.
There was a problem hiding this comment.
I've changed the title of the column from "Stabilization" to
...the first order convergence of the stabilization
examples/step-90/doc/results.dox
Outdated
| the first order convergence of the stabilization term $s_h(u_h, u_h)$. | ||
|
|
||
| <h3>Parallel scalability</h3> | ||
| In progress... |
There was a problem hiding this comment.
Let me know once you have this part too.
There was a problem hiding this comment.
First, thanks a lot for writing this program!
Here is some commentary about the program itself. I made my way down to the start of the main class, but it's late and I'll leave it here for now. I'll come back to it in a few days; go ahead and see whether you can't already address some of these comments.
As a general rule (which we have not always followed in all tutorial programs, but have found correct anyway): We've generally found that it is difficult for people to see the overall flow of a program if comments are interspersed too frequently. As a consequence, we've tried to move comments so that they cover larger blocks of code, rather than just a few lines at a time. I think step-5 is, for instance, a good example of what works well (https://dealii.org/developer/doxygen/deal.II/step_5.html#step_5-CommProg) whereas step-7 is not (https://dealii.org/developer/doxygen/deal.II/step_7.html#step_7-CommProg).
In your program, you often have 1-3 line long comments break apart larger blocks of code. That would be the right strategy if the goal were to just document a piece of code by itself (i.e., a code that is not run through doxygen to generate a website), but it's going to make it hard to read for users and keep the bigger picture in mind. Take a pass through your program and see whether perhaps some of the commentary within a function could be moved to a larger comment up front.
examples/step-90/step-90.cc
Outdated
| // We will grant to some cells empty finite element spaces FE::Nothing as done | ||
| // in step-85, | ||
| // but this time active DOFS will be only assigned to cell which are | ||
| // intersected by the surface approximation. |
There was a problem hiding this comment.
Use consistent indentation. Throughout the documentation, we try to consistently spell things "DoF" and "DoFs". Perhaps do a search-replace here too.
There was a problem hiding this comment.
Search-replace is done:)
examples/step-90/step-90.cc
Outdated
| template <class Iterator> | ||
| void cell_worker_assemble(const Iterator &cell, |
There was a problem hiding this comment.
Let's call this argument IteratorType to make clear it's a type. Do you need this to be a template parameter, though, or would typename DoFHandler<dim>::cell_iterator (or ...::active_cell_iterator) be enough?
There was a problem hiding this comment.
The stabilization should be evaluated on active cells only so typename DoFHandler<dim>::active_cell_iterator is enough.
Thank you, Wolfgang. Are you okay with the general idea/setup/test problem? |
|
Yes, sure. I haven't worked my way through the rest of the program (and was hoping that @vyushut could go through the rest and apply some of the general comments first), but I think it's a nice method to solve a nontrivial problem. With a bit more background in the introduction, this should be a nice addition to the tutorial! |
examples/step-90/doc/intro.dox
Outdated
| @f{equation*} | ||
| -\Delta_\Gamma u + u = f \qquad \text{in }\, \Gamma, | ||
| @f} | ||
| Here we added the term $u$ to the left-hand side so the problem becomes well-posed. |
There was a problem hiding this comment.
It's ok to leave it as is, but do want to point out that that's a mathematician's approach :-) From an applied perspective, the goal is to solve equations that are relevant because they describe real situations in physics, chemistry, or some other discipline. The equation is what it is, i.e., what came out of a modeling process. Whether or not we like the equation is not a relevant question, and similarly we cannot just choose to add a term to an equation to make the problem easier to solve.
That said, this is a tutorial in which we teach techniques, not solve real applications, and so I'm ok with this statement :-)
There was a problem hiding this comment.
I'm not sure if it helps, but I expanded on why the extra term is included.
We would like to solve the simplest possible problem defined on a surface, namely the Laplace--Beltrami equation, @f{equation*} -\Delta_\Gamma u + c u = f \qquad \text{in }\, \Gamma, @f} where we take $c=1$ for concreteness. We added the term $cu$ to the left-hand side so the problem becomes well-posed in the absence of any boundary; an alternative could be to take $c=0$ but impose the zero mean condition.
There was a problem hiding this comment.
Yes, I think this is great. Thanks for indulging me :-)
|
@vyushut Thanks for working through these comments. Leave a comment here (at the top level of the PR) once you've worked through everything and pushed a new version. I will then look through the remainder of the program. |
There was a problem hiding this comment.
I like the tutorial 👍 I left some comments.
Maybe @simonsticko could also review the tutorial after you have integrated @bangerth's and my comments.
examples/step-90/doc/intro.dox
Outdated
| <h3>A non-trivial surface</h3> | ||
| One of the possible bottlenecks in achieving the two-dimensional complexity while solving a surface PDE with TraceFEM | ||
| is the surface approximation. The level-set function must be defined in a three-dimensional domain and we will create | ||
| a mesh which is localized near the zero contour line, i.e near the surface, to restore the two-dimensional complexity. |
There was a problem hiding this comment.
| a mesh which is localized near the zero contour line, i.e near the surface, to restore the two-dimensional complexity. | |
| a mesh along the zero contour line to restore the two-dimensional complexity. |
There was a problem hiding this comment.
I wouldn't say "a mesh along the zero contour line" as it indicates that we know in advance which cells need to be created and refined. Instead the mesh is constructed iteratively and a good part of tutorial is devoted to it. I suggest the following:
...a mesh which is localized near the zero contour line of the level set function, i.e near the surface...
examples/step-90/doc/intro.dox
Outdated
|
|
||
| <h3>A non-trivial surface</h3> | ||
| One of the possible bottlenecks in achieving the two-dimensional complexity while solving a surface PDE with TraceFEM | ||
| is the surface approximation. The level-set function must be defined in a three-dimensional domain and we will create |
There was a problem hiding this comment.
| is the surface approximation. The level-set function must be defined in a three-dimensional domain and we will create | |
| is the surface approximation. The level-set function must be defined in a three-dimensional domain and we will create implicitly |
There was a problem hiding this comment.
I've rephrased the entire paragraph
examples/step-90/step-90.cc
Outdated
| // The normal-gradient volume stabilization form needs a bulk cell | ||
| // integration while other types of stabilization may need face | ||
| // quadratures, for example. So we check it first. | ||
| if (stabilization_scheme.needs_cell_worker()) |
There was a problem hiding this comment.
The name is chosen unfortunately!? I guess it is needs_stabilization()!? Also, I am not sure if we need a separate class for the stabilization. You only have a single implementation... Let's keep the cope simple and inline the functions here.
There was a problem hiding this comment.
You are right, there is only one implementation of the stabilization. The method needs_cell_worker() indicates that a bulk cell quadrature must be constructed in addition to the non_matching_fe_values constructed to integrate over implicitly given intersections. Hence the name.
However, there are different TraceFEM stabilizations, e.g. the so called ghost penalty term, which requires fe_values objects ( on faces AND\OR on cells depending on the fe degree) . I think that having a separate class is illuminating and it implements an important logic of separating mesh_workers that require different fe_values from the assembly().
| @@ -0,0 +1 @@ | |||
| step-85 | |||
There was a problem hiding this comment.
You also need to add the tutorial to https://github.com/dealii/dealii/blob/master/doc/doxygen/tutorial/tutorial.h.in.
Could you also add some references, e.g., to you paper?
There was a problem hiding this comment.
I've added some references.
Sure. Please tell me if/when you want more comments @vyushut . |
|
Hello @peterrum . Thank you for the comments! I will them incorporate them and let you know. |
vyushut
force-pushed
the
step86
branch
2 times, most recently
from
c4eeccd
to
68dafc8
Compare
examples/step-90/step-90.cc
Outdated
| computing_timer.print_summary(); | ||
| computing_timer.reset(); | ||
| if (cycle < convergence_levels - 1) | ||
| { | ||
| mark_intersected(); | ||
| refine(); | ||
| } | ||
| else | ||
| output_solution(); | ||
| computing_timer.print_summary(); | ||
| computing_timer.reset(); |
There was a problem hiding this comment.
Is it necessary to output timing statistics twice? I think it might be nicer to just do it once at the end of the loop.
There was a problem hiding this comment.
I started outputting it twice because, when the number of DoFs exceeds capabilities, it is the refine() method that crashes first. I would like to see the timing for solve(), assemble(), etc., for the largest #DoFs in the MPI tests, but we can safely omit the first output. Will do that.
There was a problem hiding this comment.
Yes, that happens frequently: You do certain things when you're developing the code, but then it gets removed again for the final version. (Remember that the tutorials are teaching tools, so they should really only contain what's necessary to illustrate what it is you want to teach, not be a reflection of its development history.)
|
|
||
| The results of the convergence study are shown in the following table. | ||
|
|
||
| | Cycle | DOFS | Rate | Iterations | $L^2$-Error | EOC | $H^1$-Error | EOC |$s_h^{1/2}(u_h)$| EOC | |
There was a problem hiding this comment.
Can you explain what EOC means?
There was a problem hiding this comment.
Explain in the tutorial that EOC stands for "experimental order of convergence"?
There was a problem hiding this comment.
Yes, I didn't know what the abbreviation meant :-)
|
@vyushut When you address comments, you don't have to reply to all of them. If there is no question, or if you know what to do, just do it -- no acknowledgment necessary. |
examples/step-90/step-90.cc
Outdated
| computing_timer.print_summary(); | ||
| computing_timer.reset(); | ||
| if (cycle < convergence_levels - 1) | ||
| { | ||
| mark_intersected(); | ||
| refine(); | ||
| } | ||
| else | ||
| output_solution(); | ||
| computing_timer.print_summary(); | ||
| computing_timer.reset(); |
There was a problem hiding this comment.
Yes, that happens frequently: You do certain things when you're developing the code, but then it gets removed again for the final version. (Remember that the tutorials are teaching tools, so they should really only contain what's necessary to illustrate what it is you want to teach, not be a reflection of its development history.)
|
|
||
| The results of the convergence study are shown in the following table. | ||
|
|
||
| | Cycle | DOFS | Rate | Iterations | $L^2$-Error | EOC | $H^1$-Error | EOC |$s_h^{1/2}(u_h)$| EOC | |
There was a problem hiding this comment.
Yes, I didn't know what the abbreviation meant :-)
There was a problem hiding this comment.
Separate from the other comments, could I ask you to change the names of some of the functions? The tutorial works as a teaching tool in parts because many of the programs look very much alike, allowing readers to focus on what is specific about a program. As a consequence, would you mind renaming refine->refine_grid, asemble->assemble_system, evaluate_error->compute_error.
All of my comments is small stuff. I like this program -- well done!
vyushut
force-pushed
the
step86
branch
4 times, most recently
from
e3c5a98
to
ac396f2
Compare
|
@vyushut Let us know when you're through addressing all of the change requests! |
|
The.deal.II.Library_.The.step-90.tutorial.program.pdf |
I think they look reasonable but we probably should iterate a little bit on the graphical presentation and add a little more interpretation to the results:
Maybe we can get the tutorial merged (if all code issues are addressed) and then do this as an update? |
I agree. |
ok!
I couldn't decrease the number of proc from 1000 by a factor of 4 because of the storage limitations. Instead, I ran the 200M DoFs problem with 500 proc. I can make an intermediate data point with 2000 proc to make the spacing more even. This can wait until a major rerun of parallel tests, can't it?
No, I only timed the 5 methods: assembly(), eval_errors(), refine(), setup_level_set(), solve() which are the most time consuming. In addition, eval_errors() is not shown since it, as expected, follows assembly() plots closely.
You mean " / " instead of a comma?
Nicer in which sense? Something like this?
I think it is ready, but further parallel tests will be warranted. |
I was thinking a data point for 1792 cores is missing in the right picture. Looks good otherwise. |
|
What is missing here still? I think it would be very nice if we could wrap this up this week. |
|
Hi Wolfgang! I believe it is ready. New parallel tests will come as an update. |
|
OK. If you can rearrange the commits into one or a small number of more self-contained commits, we can merge. |
…stabilized TraceFEM using adaptive refinement and meshworkers. (Squashed after the review process)
|
Thank you -- good work! |
Adds a new tutorial step-90 for the trace finite element method (TraceFEM).$Q_2$ discrete level-set and a $Q_1$ solution are possible, but both should be based on the same bulk triangulation.
Two aspects are of our focus. First, the surface approximation is separated from the discretization of the surface PDE,
e.g. a
Second, we make sure that performance of TraceFEM in the parallel implementation corresponds to that of a classical
fitted FEM for a two-dimensional problem.