step-86: Use reasonable time stepping for a heat equation solver. #15540
Conversation
bangerth
changed the title
|
@bangerth this is the mfem code i was referring to ealrier today. Ifunction Ijacobian |
|
Status for the night: I got step-26 converted to PETSc vectors and run (which requires #15577), but something is not right with hanging nodes. I'll get back to it tomorrow. |
I cannot see the code for step 26. Did you forget to push? |
|
@stefanozampini I misspoke. The code is in the step-86 directory (which is currently the PETSc version of step-26). |
|
The hanging node issue is because I commented out the line |
|
@stefanozampini I'm looking through the documentation. For SUNDIALS' ARKode, see https://dealii.org/developer/doxygen/deal.II/classSUNDIALS_1_1ARKode.html, the documentation presents the ODE to be solved with a mass matrix |
Petsc ts does not know anything about the mass matrix. Everything is implemented in terms of implicit residual and implicit jacobian, as per the documentation I wrote for deal.ii. see e.g here https://github.com/stefanozampini/petscopt/blob/73e4f9d50b1d216ff21d9322b5f83e8ca247ec24/src/mfemopt/modelproblems.cpp#L367 for an example using mfem |
|
The implicit form of the ode is of the form |
|
Ah, ok. Thanks for pointing me in the right direction! |
|
@stefanozampini I didn't get to this today, and am traveling for the next 3 days. I did write most of the introduction. I'll eventually get around to finishing it sometime soon I hope. |
@bangerth We (with @luca-heltai ) can continue on this in the next days if you are ok |
|
@stefanozampini Absolutely. I would love that. I'm out for the moment, please go ahead! |
|
@stefanozampini , @bangerth my stab at using the TimeStepper class |
|
@stefanozampini, this is how it should look like for both in-homogeneous BC, and time dependent rhs. Unfortunately, I can't seem to make the boundary conditions work, The main difference is that the system is solved by TS (actually by SNES), not by deal.II, and we have no control over what is passed to the solve routine. In particular, essential boundary conditions for linear problems result in a rhs term (that we discussed extensively), but I have nowhere I can tell PETSc that when this should happen (except perhaps when assembling the residual, but that's not right, since the system that is solved by the nonlinear solver is not solving for a solution, but for an update w.r.t. to the solution (i.e., it should have zero BC).
I think this is now solved, but not in the way I wanted. Using |
|
With last commit, one obtaines the movie above. This is now ok, but I discovered something (a bug? a feature?) that took some time for me to understand. |
|
Ok, we now have |
|
I'll make it work in parallel later tonight. |
examples/step-86/step-86.cc
Outdated
| @@ -90,7 +90,8 @@ namespace Step86 | |||
| DoFHandler<dim> dof_handler; | |||
| IndexSet hanging_nodes; | |||
There was a problem hiding this comment.
You can remove this now
|
I believe we are in quite good shape. @bangerth @luca-heltai Do you have any objection if I rebase against the later master and squash it into two commits? (One for TS changes, one for step-86). I think we have converged with the needed TS changes and I would like to open a separate PR for that. For step-86, I believe we only want to implement a MMS test and maybe make it nonlinear too. |
|
Go ahead with the rebase! As for nonlinear -- I'm ok with solving the same problem as in step-26 (+ having non-constant boundary values), but I then it might also be nice to make the problem more complicated by adding a nonlinearity. What do you have in mind? |
|
I guess @luca-heltai was proposing Allen-Cahn |
|
Video of the solution (not too exciting, but that's not the point) is here: https://www.youtube.com/watch?v=fhJVkcdRksM |
|
/rebuild |
There was a problem hiding this comment.
I am proposing some editing but it looks quite good to me. Thanks @bangerth
examples/step-86/doc/intro.dox
Outdated
| do by oneself. In the current case, deal.II has interfaces to two | ||
| such libraries: SUNDIALS, the *SUite of Nonlinear and DIfferential/ALgebraic | ||
| equation Solvers* (and here specifically the Runge-Kutta-type solvers | ||
| wrapped in the SUNDIALS::ARKode class), and PETSc's TS sub-packaged |
There was a problem hiding this comment.
| wrapped in the SUNDIALS::ARKode class), and PETSc's TS sub-packaged | |
| wrapped in the SUNDIALS::ARKode class), and PETSc's TS sub-package |
examples/step-86/doc/intro.dox
Outdated
| On the other hand, when using the PETScWrappers::TimeStepper class that | ||
| wraps the PETSc TS sub-package, you will find that when stating things | ||
| as a typical ODE, it does not like the presence of a mass matrix. But | ||
| it can solve ODEs that are stated in "implicit" form, and in that | ||
| case we simply bring everything to the left hand side and obtain |
There was a problem hiding this comment.
| On the other hand, when using the PETScWrappers::TimeStepper class that | |
| wraps the PETSc TS sub-package, you will find that when stating things | |
| as a typical ODE, it does not like the presence of a mass matrix. But | |
| it can solve ODEs that are stated in "implicit" form, and in that | |
| case we simply bring everything to the left hand side and obtain | |
| On the other hand, when using the PETScWrappers::TimeStepper class, | |
| we can solve ODEs that are stated in a general "implicit" form, and in that | |
| case we simply bring everything to the left hand side and obtain |
examples/step-86/doc/intro.dox
Outdated
| a small abuse of terminology. In the current case, this matrix | ||
| is $J=A + \alpha M$. If you have read through step-26, it is probably | ||
| not lost on you that this matrix appears prominently there as well -- | ||
| with $\alpha$ being a multiple of the inverse of the time step $k_n$. |
There was a problem hiding this comment.
| with $\alpha$ being a multiple of the inverse of the time step $k_n$. | |
| with $\alpha$ being a multiple of the inverse of the time step. |
examples/step-86/doc/intro.dox
Outdated
| - A callback that is called at the end of each time step (or perhaps | ||
| at other intervals) and that is provided with the current solution | ||
| and other information and that can be used to "monitor" the progress | ||
| of computations. One of the ways in which this can be used is to | ||
| output visualization data every few time steps. |
There was a problem hiding this comment.
| - A callback that is called at the end of each time step (or perhaps | |
| at other intervals) and that is provided with the current solution | |
| and other information and that can be used to "monitor" the progress | |
| of computations. One of the ways in which this can be used is to | |
| output visualization data every few time steps. | |
| - A callback that is called at the end of each time step | |
| and that is provided with the current solution | |
| and other information that can be used to "monitor" the progress | |
| of computations. One of the ways in which this can be used is to | |
| output visualization data every few time steps. |
examples/step-86/doc/intro.dox
Outdated
| the algebraic degrees of freedom are listed. At the end of each time | ||
| step (or perhaps also at other times), a second callback then needs to |
There was a problem hiding this comment.
| the algebraic degrees of freedom are listed. At the end of each time | |
| step (or perhaps also at other times), a second callback then needs to | |
| the algebraic degrees of freedom are listed. At the end of each solution stage of the solver, a second callback then needs to |
| PETScWrappers::PreconditionBoomerAMG preconditioner; | ||
| preconditioner.initialize(jacobian_matrix); | ||
| #else | ||
| PETScWrappers::PreconditionSSOR preconditioner; |
There was a problem hiding this comment.
PETSc has his native AMG ('-pc_type gamg'). It would make sense to use it here (but it is not wrapped), since SSOR is a bad choice. Missing deal.II interfaces to PETSc solvers is why I prefer command line options solvers. Once the optimal is found, then it makes sense to add a solve_with_jacobian callback
There was a problem hiding this comment.
Yes, fair point. I augmented the comment in front of the function.
| // not to refine is done in the lambda function that calls the current | ||
| // function.) | ||
| // | ||
| // Breaking things into a "mark cells" function and into a "execute mesh |
There was a problem hiding this comment.
This is because you are not considering the point of view of the solver, and you are not discussing the actual callback, the decide_for_coarsening_and_refinement. For me, it is natural to first take a decision, inform the solver about that, and then the solver can do what it needs to be done before performing the actual remesh + transfer
There was a problem hiding this comment.
I was thinking about why you structured it this way. I think the simplest way would have been if the decide_... function really only does that: Decide whether we want to do mesh refinement. Then the marking along with the actual refinement would happen in the second callback, which would also be natural (keep marking and refinement together). But I suspect that you didn't do it that way because you don't have the solution vector available in the second callback (you only have a whole bunch of vectors), and so the marking needs to happen in the first callback.
I don't know whether it would be possible to query the TS object for the current solution vector? If that were the case, we could provide that in the second callback.
There was a problem hiding this comment.
This is the relevant PETSc code https://gitlab.com/petsc/petsc/-/blob/main/src/ts/interface/ts.c?ref_type=heads#L3841. resizesetup runs the deal.ii decide callback (actually, I wrote this in PETSc to support this from deal.II). The solution vector is the first vector of the list when doing the transfer. It looked clean to me, but if you think a different approach would be cleaner, please change it.
bangerth
force-pushed
the
step-86
branch
2 times, most recently
from
c8f5385
to
5f74525
Compare
|
Thanks for your feedback, @stefanozampini ! I've made the changes as requested, plus the individual comments. |
|
-snes_monitor, with underscore
…On Mon, Jul 1, 2024, 04:07 Wolfgang Bangerth ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In examples/step-86/doc/results.dox
<#15540 (comment)>:
> + 5 linear iterations.
+ 8 linear iterations.
@stefanozampini <https://github.com/stefanozampini> I think I need help
from a PETSc expert :-) -snes-monitor has no effect here, I get no
additional output :-(
—
Reply to this email directly, view it on GitHub
<#15540 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/ABHCKANKN7O5K3WQT7GK54TZKC2WJAVCNFSM6AAAAABJRY7PN2VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCNJQGI4DCNRVGM>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
Ah, yes. This produces the following: This was useful information -- thank you for pointing to that command line flag! It is an interaction between the linear solver (which solves to an accuracy of which apparently pertain to the ODE error, not the nonlinear error. |
|
@bangerth Would you mind to fix the doxygen issue? |
|
OK, I've rebased this now and it's ready from my side. |
|
Oh yay! Thank you, @stefanozampini and @luca-heltai for your work on this program! |
step-26 uses a home-grown implicit Euler time stepper. We should really show how one can do that better, using the PETSc TS or SUNDIALS ARKode wrappers.
The current state of the PR is just a shell, but I wanted to open the PR to facilitate collaboration. Among the options I'd be open is to also replace step-52, which uses "better" time steppers, but nothing as sophisticated as the ones mentioned above. step-52 also has a number of known problems, see #7191 and #9684. Replacing step-52 is something I'll have to discuss with @Rombur tonight or tomorrow.
@stefanozampini FYI.