Clarify construction of Jacobian in Navier-Stokes tutorial#1138
Clarify construction of Jacobian in Navier-Stokes tutorial#1138jfdev001 wants to merge 15 commits intoFerrite-FEM:masterfrom
Conversation
In the incompressible Navier-Stokes tutorial, there are two things wish I noticed that might be improved: P1. The construction of the Jacobian is a bit hand-wavy in the sense that there is a claim that the Jacobian and stiffness matrix share the same sparsity pattern. P2. It is not necessarily clear what motivates the `navierstokes_jac_element!` function. To address these points, I have added a small derivation for the analytical form of the Jacobian. In particular, I note that the nonlinear advection term must be handled explicitly when constructing the Jacobian. The functions `navierstokes_jac!` and `navierstokes_jac_element!` already implement what I describe in the derivation, but to make things explicit, I believe the derivation is worth including. The derivation is currently in the section `Commented Implementation > Solution of the semi-discretized system via DifferentialEquations.j` and may be more appropriate in the theory section.
termi-official
left a comment
There was a problem hiding this comment.
Thanks for this proposal.
The construction of the Jacobian is a bit hand-wavy in the sense that there is a claim that the Jacobian and stiffness matrix share the same sparsity pattern.
From reading over your suggestion, your actual question here is: How do we compute the Jacobian for the right hand side? (See comments in the diff)
It is not necessarily clear what motivates the navierstokes_jac_element! function.
Good point. We can probably clarify be renaming this function.
| # through the `jac_prototype` argument. The sparsity pattern and values for the | ||
| # Jacobian can be shown by taking the directional gradient (see also [deal.ii: step 57](https://www.dealii.org/current/doxygen/deal.II/step_57.html)) of | ||
| # $F(\bold{x})$ along $\delta \bold{x}$ where $F(\bold{x}) = F(v, p)$ is the | ||
| # semi-discrete weak form of the incompressible Navier-Stokes equations such that | ||
| # ```math | ||
| # F(v, p) = \begin{pmatrix} | ||
| # F_1 \\ | ||
| # F_2 | ||
| # \end{pmatrix} = \begin{pmatrix} | ||
| # - \int_\Omega \nu \nabla v : \nabla \varphi - \int_\Omega (v \cdot \nabla) v \cdot \varphi + \int_\Omega p (\nabla \cdot \varphi) \\ | ||
| # \int_{\Omega} (\nabla \cdot v) \psi | ||
| # \end{pmatrix}. | ||
| # ``` | ||
| # The variable $t$ is omitted as it is not relevant to the directional gradient. | ||
| # The directional gradient is thus given by | ||
| # ```math | ||
| # \begin{aligned} | ||
| # \nabla F(v, p) \cdot (\delta v, \delta p) &= \lim_{\epsilon \to 0 } \frac{1}{\epsilon} (F(v + \delta v, p + \delta p ) - F(v, p)) \\ | ||
| # &= \lim_{\epsilon \to 0} \frac{1}{\epsilon} \begin{pmatrix} | ||
| # F_1(v + \epsilon \delta v, p + \epsilon \delta p) - F_1(v,p) \\ | ||
| # F_2(v + \epsilon \delta v) - F_2(v) | ||
| # \end{pmatrix} \\ | ||
| # &= \lim_{\epsilon \to 0} \frac{1}{\epsilon} \begin{pmatrix} | ||
| # - \int_\Omega \nu \nabla (v + \epsilon \delta v) : \nabla \varphi - \int_\Omega ((v + \epsilon \delta v) \cdot \nabla) (v + \epsilon \delta v) \cdot \varphi + \int_\Omega (p + \epsilon \delta p)(\nabla \cdot \varphi) - F_1 \\ | ||
| # \int_{\Omega} (\nabla \cdot (v + \epsilon \delta v)) \psi - F_2 | ||
| # \end{pmatrix} \\ | ||
| # &= \lim_{\epsilon \to 0} \frac{1}{\epsilon} \begin{pmatrix} | ||
| # - \int_{\Omega} \nu \epsilon \nabla \delta v : \nabla \varphi - \int_{\Omega} (\epsilon \delta v \cdot \nabla v + v \cdot \epsilon \nabla \delta v + \epsilon^2 \delta v \cdot \nabla \delta v) \cdot \varphi + \int_{\Omega} \epsilon \delta p (\nabla \cdot \varphi) \\ | ||
| # - \int_{\Omega} (\nabla \cdot \epsilon \delta v) \psi | ||
| # \end{pmatrix} \\ | ||
| # &= \begin{pmatrix} | ||
| # - \int_{\Omega} \nu \nabla \delta v : \nabla \varphi - \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi + \int_{\Omega} \delta p (\nabla \cdot \varphi) \\ | ||
| # - \int_{\Omega} (\nabla \cdot \delta v) \psi | ||
| # \end{pmatrix}. | ||
| # \end{aligned} | ||
| # ``` |
There was a problem hiding this comment.
Some things here. First, the notation is totally disconnected from the theory section. Second, it makes more sense to add a new section in the theory section named "derivation of the Jacobian".
Can you move the section up (above "Commented implementation") and connect the notations properly?
| # The directional gradient of the nonlinear advection term is the only term that | ||
| # does not match the pattern described by the stiffness matrix. This implies | ||
| # that local contributions from the nonlinear elements | ||
| # ```math | ||
| # - \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi | ||
| # ``` | ||
| # must be added to the Jacobian explicitly, and $v$ must be a function value | ||
| # while $\nabla v$ must be a function gradient. | ||
| # Conversely, | ||
| # $\delta v$ and $\delta p$ are approximated by $\varphi$ and $\psi$, | ||
| # respectively, so the nonzero values of $K$ *are* the local contributions of the | ||
| # linear terms and can be used to populate the Jacobian directly. | ||
| # It is therefore simple to see that the Jacobian and the stiffness matrix share | ||
| # the same sparsity pattern, since they share the same relation between trial | ||
| # and test functions. |
There was a problem hiding this comment.
If it is still necessary after fixing the comment above, can you rephrase this? The continuous and discrete notation/theory is mixed up (plus notation).
There was a problem hiding this comment.
I will address both of your comments (#1138 (comment) and #1138 (comment)) this week/weekend (2025-02-10). Just confirming that I am actively working on it.
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## master #1138 +/- ##
==========================================
- Coverage 93.75% 93.61% -0.14%
==========================================
Files 39 39
Lines 6241 6267 +26
==========================================
+ Hits 5851 5867 +16
- Misses 390 400 +10 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
| # * The solution of nonlinear problems with newton-like methods (deal.ii, bueller) | ||
| # * The Jacobian's relevance, provide the newton update equation (deal.ii, bueller, heath) |
There was a problem hiding this comment.
Please take note that no Newton is involved in this example.
There was a problem hiding this comment.
Good point. I see now that the solver is Rodas5P, which relies on the jacobian (see eq 1.3 in Steinebach G. Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks within the Julia Differential Equations package. In: BIT Numerical Mathematics, 63(2), 2023 from OrdinaryDiffEqRosenbrock.Rodas5P) but is not a newton method. I will correct this.
There was a problem hiding this comment.
We do not need to go over board with details here, as we really just want to showcase here that there is not too much to do to use a wide range of solvers from DiffEq.
* Previously, the notation I used to describe the jacobian derivation did not match with what was introduced in the tutorial. I have updated my notation to more accuately reflect the notation that was already introduced. * Previously, the continuous/discrete notation/theory was unclear. I have explicitly added how the finite element approximation of the jacobian can be used to fill its values and sparsity pattern. * Stylistically, the headers/subheaders were not following a uniform case system (i.e., sometimes all first letter capital, sometimes not). To conform with the case system used by (some) other tutorials, only the first letter is capitalized. * Removed the reference to Newton's method as it is not relevant Regarding the header/subheader captilization, this is not internally consistent in the Ferrite docs anyways (maybe I open a new issue for this?). For example, [Stokes Flow Tutorial](https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/stokes-flow/) does not capitalize the first letters of the `## My header` while [Transient Heat Equation Tutorial](https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/transient_heat_equation/) does `## My Header`.
Please go ahead. :) I will need a few days to come back to this unfortunately, sorry for the delay. |
| # To enable the use of a wide range of solvers, it is more efficient to provide | ||
| # information to the solver about the sparsity pattern and values of the Jacobian | ||
| # of $f(u, t)$. The sparsity pattern and values for the Jacobian can be shown by | ||
| # taking the [directional gradient](https://en.wikipedia.org/wiki/Directional_derivative) of |
There was a problem hiding this comment.
What is a "directional gradient"? I could not find the term in the link.
There was a problem hiding this comment.
It should be directional derivative, I have pushed a change correcting this discrepancy.
termi-official
left a comment
There was a problem hiding this comment.
Okay, maybe we need to go one step back here. What exactly is the problem you are trying to solve here? Maybe I am missing something, but if you want to show how the Jacobian and Stokes operator relate to each other, then we can compute the Jacobian just directly without iterating over theory on the function space again.
The Jacobian here is per definition really just
| # &= \lim_{\epsilon \to 0} \frac{1}{\epsilon} \begin{bmatrix} | ||
| # - \int_\Omega \nu \nabla (v + \epsilon \delta v) : \nabla \varphi - \int_\Omega ((v + \epsilon \delta v) \cdot \nabla) (v + \epsilon \delta v) \cdot \varphi + \int_\Omega (p + \epsilon \delta p)(\nabla \cdot \varphi) \\ | ||
| # \int_{\Omega} (\nabla \cdot (v + \epsilon \delta v)) \psi | ||
| # \end{bmatrix} - f(v, p) \\ |
There was a problem hiding this comment.
| # \end{bmatrix} - f(v, p) \\ | |
| # \end{bmatrix} - f(v, p) \\ |
is missing the eps or braces
| # \delta p_h(\mathbf{x}) &= \sum_{i}^{N} \psi_i(\mathbf{x}) \delta \hat{p}_i. | ||
| # \end{aligned} | ||
| # ``` | ||
| # on a mesh with element size $h$ and $N$ trial (aka nodal shape) functions |
There was a problem hiding this comment.
| # on a mesh with element size $h$ and $N$ trial (aka nodal shape) functions | |
| # on a mesh with characteristic element size $h$ and $N$ trial (aka nodal shape) functions |
| # We now show that the finite element contributions for *part* of the Jacobian | ||
| # exactly match the contributions for the Stokes operator $K$. When substituting | ||
| # the finite element approximation for the viscosity term and comparing this | ||
| # with the analagous term in the directional gradient, | ||
| # ```math | ||
| # \begin{aligned} | ||
| # - \int_\Omega \nu \nabla \delta v : \nabla \varphi & \xrightarrow{\delta v_h} - \sum_{i}^{N} (\int_{\Omega} \nu \nabla \varphi_i : \nabla \varphi_j) \delta \hat{v}_i \\ | ||
| # - \int_\Omega \nu \nabla v : \nabla \varphi & \xrightarrow{v_h} - \sum_{i}^{N} (\int_{\Omega} \nu \nabla \varphi_i : \nabla \varphi_j) \hat{v}_i \\ | ||
| # \end{aligned} | ||
| # ``` | ||
| # we see that the coefficients (i.e., the terms in the integrand on the right | ||
| # hand side of the arrow) for the nodal unknowns $\delta \hat{v}_i$ and $\hat{v}_i$ | ||
| # are the same. The same can be shown for the pressure term and the |
There was a problem hiding this comment.
We are mixing up many things here and I am also not sure what these limits should be...
There was a problem hiding this comment.
This is just a notational error. I am not taking a limit. I wanted to indicate that when one substitutes the expression over the arrow (e.g.,
This is definitely the most succinct way of communicating what the Jacobian is, but given that it's a tutorial, I wanted to provide lower level details that correspond more directly to the function
|
|
Thanks for elaborating. So, we can probably do it similar as in the hyperelasticity tutorial here: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/hyperelasticity/ , i.e. just having the main points in the text and derivation details on the underlying function space in spoilers. What do you think? |
|
Thank you for working with me on this. It is my first proper open-source contribution beyond trivial documentation PRs, so I appreciate your patience! I think the style in https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/hyperelasticity/ is perfect. I will rewrite my PR in that spoilers style. Update: Sorry for delay, Have been busy, I will make the changes by the end of 2025-03. |
|
No problem, we have to thank for the help to clarify the tutorial for new users! |
Also added punctuation to maths equations in `ns_vs_diffeq.jl` Updated `tutorials/index.md` with the correct numbering of tutorials.
|
Two comments that are unrelated to this direct PR but which might warrant opening separate issues related to the Navier-Stokes tutorial:
|
|
… on implementing this operation
|
@termi-official Any thoughts on my comments here #1138 (comment) ? Regarding opening another issue/PR for
if we include this in the current Navier-Stokes tutorial, it would eliminate the "need" to include some steps in the "how-to" mentioned in #1139. Incidentally, the current PR is ready for review. Is there a better way to indicate that a PR is ready for review other than just converting the draft PR to a normal PR? I see that, for example, #1183 has the awaiting review label, but it doesn't look like I can assign such a label. |
This is correct. The label is for us to notify Fredrik that we need him. |
termi-official
left a comment
There was a problem hiding this comment.
There are multiple inconsistencies. Can you check the comments I left?
| # $\frac{\delta N}{\delta u}$ is therefore given by | ||
| # | ||
| # ```math | ||
| # - \sum_{i}^{N} (\int_{\Omega} ((\varphi_i \cdot \nabla) v + (v \cdot \nabla) \varphi_i) \cdot \varphi_j) \delta \hat{v}_i, |
There was a problem hiding this comment.
I think we should specify this further, as this expression is not correct in its dimension. The left side is dim(u) and the right side is dim(v). Also the continuous and discrete variables are mixed.
| # - \sum_{i}^{N} (\int_{\Omega} ((\varphi_i \cdot \nabla) v + (v \cdot \nabla) \varphi_i) \cdot \varphi_j) \delta \hat{v}_i, | |
| # \frac{\partial N}{\partial u} = - \sum_{i}^{N} (\int_{\Omega} ((\varphi_i \cdot \nabla) v + (v \cdot \nabla) \varphi_i) \cdot \varphi_j) \cdot \delta \hat{v}_i, |
| # sparsity values and pattern can be used to initialize the Jacobian. | ||
| # | ||
| # To define $\frac{\partial N}{\partial u}$, we take the directional derivative | ||
| # of $f$ along the vector $(\delta v, \delta p)$ and |
There was a problem hiding this comment.
I think we are mixing the derivative at the continuous weak form with the derivative of the discrete problem.
| # of $f$ along the vector $(\delta v, \delta p)$ and | ||
| # define finite element approximations of $\delta v$ and $\delta p$ with trial | ||
| # functions $\varphi$ and $\psi$, respectively. The discrete form of | ||
| # $\frac{\delta N}{\delta u}$ is therefore given by |
There was a problem hiding this comment.
| # $\frac{\delta N}{\delta u}$ is therefore given by | |
| # $\frac{\partial N}{\partial u}$ is therefore given by |
I also think I see where the confusion comes from. We usually use
There was a problem hiding this comment.
Working on this now. Need to make more adjustments to 54fbfbd since it's still off, please ignore it since I need to rethink it.
There was a problem hiding this comment.
So to be clear with the notation, I am assuming the following:
- The test functions
$\varphi$ and$\psi$ for$v$ and$p$ , respectively, are$\varphi \in H^{1}_{0}$ and$\psi \in L^2(\Omega)$ where
where
2. The finite element approximation of
3. We restrict the test functions to the finite set of basis functions
This implies to me that the only place I need to use the hat notation (e.g.,
would be correct, yes? I know I still need to resolve the dimensionality comment #1138 (comment)!
The use of
There was a problem hiding this comment.
the right hand side of ... would be correct, yes?
No. That expression does not even make sense. You have matrix with (N+Mx(N+M) entries = vector with N entries. If you really want to learn how to do this stuff, then I can just highly recommend to take some time, a pen and some paper and do the derivation step by step yourself without peeking into the tutorial.
There was a problem hiding this comment.
You're completely right. I am doing some more reading and re-writing. Thank you for the feedback, I am on vacation now and will send some corrections soon as a comment (not as a commit) :))
There was a problem hiding this comment.
Sorry, took me really forever to come back to this one.
Indeed. We typically use the continuous form in the analysis of the equations and you can find often that also the weak form is linearized directly (in contrast to linearizing the discretized form). IIRC both approaches should give you the same result for consistency reasons. I usually prefer the discretized form to be used directly, as it is easier to digest for people outside formal training in math, as we only have to deal with the "simple" finite dimensional case. I am fine also with the explanation you give tho. It comes down to a matter of "taste" anyway.
Under these conditions, the Gateaux derivative is the continuous linear operator represented by the Jacobian matrix---this definition is on pages 57-58 from The Differentiation and Integration of Nonlinear Operators in Nonlinear Functional Analysis and Applications (Rall 1971).
Not quite. In your reference Tapia states this for the finite dimensional (i.e. discrete) case, not for the infinite-dimensional (i.e. continuous) case.
[...] common to simply use autodiff [...]
We could probably just use autodiff here, this is correct. I do not know if the performance will be on par tho and if the resulting system will be correctly constrained (Dirichlet boundary conditions). Have you by any chance tried if it works right now?
There was a problem hiding this comment.
No problem! Thanks for the feedback! I did not try with autodiff. I'll come back to this PR during the holidays (after 2025-12-23), had to step away for a while, thanks for your patience.
There was a problem hiding this comment.
If you have written something somewhere here, I have gotten some notificaiton on Github, but I think the comment got gobbled up by some Github outage.
There was a problem hiding this comment.
Hi @termi-official, I've been quite caught up with work/application season right now, so will have to step away from Ferrite longer than I anticipated. I am not sure when I would be able to return to this PR. I can either leave it as a draft, or close it. Lmk what you think.
Unfortunately, same thing applies to the #1139. I would recommend unassigning me from that one so that other contributors might be able to tackle it if they want.
I appreciate all your assistance with this PR up until this point! I will try to return to it at some point, but don't know a timeline yet.
There was a problem hiding this comment.
I can fully understand. And no worries, I won't run away in the meantime. :)
| # that the trial functions are the same as in the previous [section](@ref weak-form-derivation). The finite element approximation for $\delta v$ is thus given by | ||
| # | ||
| # ```math | ||
| # \delta v_h(\mathbf{x}) = \sum_{i}^{N} \varphi_i(\mathbf{x}) \delta \hat{v}_i, |
There was a problem hiding this comment.
I think we never used
There was a problem hiding this comment.
Right!
may be less clear. It is therefore not necessary to define
* Specify that K has zeros but that they are ignored in the construction of the jacobian * Use word "derive" not "define" the expression for the derivative of the nonlinear term (i.e., advection) w.r.t state variables u * notation fix: \delta change to \partial * more specific: let dN/du = rhs(dN/du) instead of just having the right hand side with no left hand side Co-authored-by: Dennis Ogiermann <termi-official@users.noreply.github.com>
In the incompressible Navier-Stokes tutorial, there are two things I noticed that might be improved:
navierstokes_jac_element!function.To address these points, I have added a small derivation for the analytical form of the Jacobian. In particular, I note that the nonlinear advection term must be handled explicitly when constructing the Jacobian. The functions$K$ are the nonzeros in $J$ .
navierstokes_jac!andnavierstokes_jac_element!already implement what I describe in the derivation, but to make things explicit, I believe the derivation is worth including. I also add a comment about why the nonzeros inThe derivation is currently in the section
Commented Implementation > Solution of the semi-discretized system via DifferentialEquations.jand may be more appropriate in the theory section.