Clarify construction of Jacobian in Navier-Stokes tutorial #1138

jfdev001 · 2025-02-02T14:28:29Z

In the incompressible Navier-Stokes tutorial, there are two things I noticed that might be improved:

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.
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. I also add a comment about why the nonzeros in $K$ are the nonzeros in $J$.

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.

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

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.

termi-official · 2025-02-06T17:38:49Z

docs/src/literate-tutorials/ns_vs_diffeq.jl

+# 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}
+# ```


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?

termi-official · 2025-02-06T18:00:39Z

docs/src/literate-tutorials/ns_vs_diffeq.jl

+# 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.


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

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 · 2025-02-10T18:29:02Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 89.25%. Comparing base (43db37c) to head (06fa399).
Report is 11 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1138      +/-   ##
==========================================
- Coverage   93.59%   89.25%   -4.34%     
==========================================
  Files          39       39              
  Lines        6087     6552     +465     
==========================================
+ Hits         5697     5848     +151     
- Misses        390      704     +314

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

termi-official · 2025-02-14T08:48:31Z

docs/src/literate-tutorials/ns_vs_diffeq.jl

+# * 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)


Please take note that no Newton is involved in this example.

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.

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

termi-official · 2025-02-17T11:45:59Z

Regarding the header/subheader captilization, this is not internally
consistent in the Ferrite docs anyways (maybe I open a new issue for
this?

Please go ahead. :)

I will need a few days to come back to this unfortunately, sorry for the delay.

termi-official · 2025-02-18T15:05:32Z

docs/src/literate-tutorials/ns_vs_diffeq.jl

+# 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


What is a "directional gradient"? I could not find the term in the link.

It should be directional derivative, I have pushed a change correcting this discrepancy.

termi-official

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

$$J(u) := \partial_u (K u + N(u)) = K + \partial_u N(u)$$

termi-official · 2025-02-18T15:24:33Z

docs/src/literate-tutorials/ns_vs_diffeq.jl

+#   &= \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) \\


Suggested change

# \end{bmatrix} - f(v, p) \\

# \end{bmatrix} - f(v, p) \\

is missing the eps or braces

termi-official · 2025-02-18T15:26:14Z

docs/src/literate-tutorials/ns_vs_diffeq.jl

+# \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 


Suggested change

# 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

termi-official · 2025-02-18T15:27:40Z

docs/src/literate-tutorials/ns_vs_diffeq.jl

+# 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


We are mixing up many things here and I am also not sure what these limits should be...

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., $\delta v_h$), the left side of the arrow becomes the right side. It would likely be clearer to write, "Substituting the finite element approximation $v_h$ for the viscoscity term and substituting the finite element approximation $\delta v_h$ for the analogous term in the directional derivative, ..."

jfdev001 · 2025-02-19T09:52:39Z

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

J ( u ) := ∂ u ( K u + N ( u ) ) = K + ∂ u N ( u )

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 navierstokes_jac! and navierstokes_jac_element!. The problem I am therefore attempting to solve is to more clearly motivate how one would write a function(s) with the Ferrite API defining a Jacobian that could be used with DifferentialEquations.jl solvers. Perhaps the most clear explanation would be to do the following:

State J ( u ) := ∂ u ( K u + N ( u ) ) = K + ∂ u N ( u ) at the beginning
From this, one can very clearly see that the sparsity pattern and values of the Jacobian correspond to element contributions of K. Hence, the nonzeros(J) .= nonzeros(K) in navierstokes_jac!. I could then remove the large part of the discussion of $K$ in my derivation.
Show the directional derivative derivation (see the directional derivative is equivalent to the Jacobian) since I think showing this is particularly didactic and is not shown in any other Ferrite tutorials involving a Jacobian.
From the derivation, highlight how ∂ u N ( u ) corresponds to $- \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi$ that this is reflected in the implementation. It should also be mentioned that the same trial function for $v$ is used for $\delta v$.

termi-official · 2025-02-19T11:59:23Z

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?

jfdev001 · 2025-02-19T13:25:22Z

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. ~~I am traveling for business currently, so may not have time until next week (20250-02-24) to make changes.~~

Update: Sorry for delay, Have been busy, I will make the changes the week of 2025-03-03

termi-official · 2025-02-19T13:55:46Z

No problem, we have to thank for the help to clarify the tutorial for new users!

termi-official self-requested a review February 3, 2025 09:41

termi-official reviewed Feb 6, 2025

View reviewed changes

brief outline for jacobian derivation reorganization

0f07aa3

added references section for clarity

ef28bea

jfdev001 marked this pull request as draft February 11, 2025 05:11

termi-official reviewed Feb 14, 2025

View reviewed changes

jfdev001 mentioned this pull request Feb 17, 2025

Docs: Consistent formatting for capitalization of headers #1156

Open

termi-official reviewed Feb 18, 2025

View reviewed changes

say directional derivative instead of directioal gradient!

06fa399

termi-official reviewed Feb 18, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Clarify construction of Jacobian in Navier-Stokes tutorial #1138

Clarify construction of Jacobian in Navier-Stokes tutorial #1138

jfdev001 commented Feb 2, 2025 •

edited

Loading

termi-official left a comment

termi-official Feb 6, 2025

termi-official Feb 6, 2025

jfdev001 Feb 10, 2025

codecov bot commented Feb 10, 2025 •

edited

Loading

termi-official Feb 14, 2025

jfdev001 Feb 14, 2025 •

edited

Loading

termi-official Feb 14, 2025

termi-official commented Feb 17, 2025

termi-official Feb 18, 2025

jfdev001 Feb 18, 2025 •

edited

Loading

termi-official left a comment

termi-official Feb 18, 2025

termi-official Feb 18, 2025

termi-official Feb 18, 2025

jfdev001 Feb 19, 2025

jfdev001 commented Feb 19, 2025 •

edited

Loading

termi-official commented Feb 19, 2025

jfdev001 commented Feb 19, 2025 •

edited

Loading

termi-official commented Feb 19, 2025

		# * 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)

	# 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

Clarify construction of Jacobian in Navier-Stokes tutorial #1138

Are you sure you want to change the base?

Clarify construction of Jacobian in Navier-Stokes tutorial #1138

Conversation

jfdev001 commented Feb 2, 2025 • edited Loading

termi-official left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Feb 10, 2025 • edited Loading

Codecov Report

Choose a reason for hiding this comment

jfdev001 Feb 14, 2025 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

termi-official commented Feb 17, 2025

Choose a reason for hiding this comment

jfdev001 Feb 18, 2025 • edited Loading

Choose a reason for hiding this comment

termi-official left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jfdev001 commented Feb 19, 2025 • edited Loading

termi-official commented Feb 19, 2025

jfdev001 commented Feb 19, 2025 • edited Loading

termi-official commented Feb 19, 2025

jfdev001 commented Feb 2, 2025 •

edited

Loading

codecov bot commented Feb 10, 2025 •

edited

Loading

jfdev001 Feb 14, 2025 •

edited

Loading

jfdev001 Feb 18, 2025 •

edited

Loading

jfdev001 commented Feb 19, 2025 •

edited

Loading

jfdev001 commented Feb 19, 2025 •

edited

Loading