Clarify construction of Jacobian in Navier-Stokes tutorial

docs/src/literate-tutorials/ns_vs_diffeq.jl

@@ -27,7 +27,7 @@ nothing                    #hide
 # discretization technique - in our case FEM. The mass matrix form of ODEs and DAEs
 # is given as:
 # ```math
-#   M(t) \mathrm{d}_t u = f(u,t)
+#   M(t) \mathrm{d}_t u = f(u,t),
 # ```
 # where $M$ is a possibly time-dependent and not necessarily invertible mass matrix,
 # $u$ the vector of unknowns and $f$ the right-hand-side (RHS). For us $f$ can be interpreted as
@@ -41,8 +41,8 @@ nothing                    #hide
 # The incompressible Navier-Stokes equations can be stated as the system
 # ```math
 #  \begin{aligned}
-#    \partial_t v &= \underbrace{\nu \Delta v}_{\text{viscosity}} - \underbrace{(v \cdot \nabla) v}_{\text{advection}} - \underbrace{\nabla p}_{\text{pressure}} \\
-#               0 &= \underbrace{\nabla \cdot v}_{\text{incompressibility}}
+#    \partial_t v &= \underbrace{\nu \Delta v}_{\text{viscosity}} - \underbrace{(v \cdot \nabla) v}_{\text{advection}} - \underbrace{\nabla p,}_{\text{pressure}} \\
+#               0 &= \underbrace{\nabla \cdot v,}_{\text{incompressibility}}
 #  \end{aligned}
 # ```
 # where $v$ is the unknown velocity field, $p$ the unknown pressure field,
@@ -59,23 +59,23 @@ nothing                    #hide
 #  \begin{bmatrix}
 #      4 v_{in}(t) y (0.41-y)/0.41^2 \\
 #      0
-#  \end{bmatrix}
+#  \end{bmatrix},
 # ```
 # where $v_{in}(t) = \text{clamp}(t, 0.0, 1.5)$. With a dynamic viscosity of $\nu = 0.001$
 # this is enough to induce turbulence behind the cylinder which leads to vortex shedding. The top and bottom of our
 # channel have no-slip conditions, i.e. $v = [0,0]^{\textrm{T}}$, while the right boundary has the do-nothing boundary condition
 # $\nu \partial_{\textrm{n}} v - p n = 0$ to model outflow. With these boundary conditions we can choose the zero solution as a
 # feasible initial condition.
 #
-# ### Derivation of semi-discrete weak form
+# ### [Derivation of semi-discrete weak form](@id weak-form-derivation)
 #
 # By multiplying test functions $\varphi$ and $\psi$ from a suitable test function space on the strong form,
 # followed by integrating over the domain and applying partial integration to the pressure and viscosity terms
 # we can obtain the following weak form
 # ```math
 #  \begin{aligned}
-#    \int_\Omega \partial_t v \cdot \varphi &= - \int_\Omega \nu \nabla v : \nabla \varphi - \int_\Omega (v \cdot \nabla) v \cdot \varphi + \int_\Omega p (\nabla \cdot \varphi) + \int_{\partial \Omega_{N}} \underbrace{(\nu \partial_n v - p n )}_{=0} \cdot \varphi \\
-#                                  0 &= \int_\Omega (\nabla \cdot v) \psi
+#    \int_\Omega \partial_t v \cdot \varphi &= - \int_\Omega \nu \nabla v : \nabla \varphi - \int_\Omega (v \cdot \nabla) v \cdot \varphi + \int_\Omega p (\nabla \cdot \varphi) + \int_{\partial \Omega_{N}} \underbrace{(\nu \partial_n v - p n )}_{=0} \cdot \varphi, \\
+#                                  0 &= \int_\Omega (\nabla \cdot v) \psi,
 #  \end{aligned}
 # ```
 # for all possible test functions from the suitable space.
@@ -105,13 +105,128 @@ nothing                    #hide
 #  \begin{bmatrix}
 #      N(\hat{v}, \hat{v}, \hat{\varphi}) \\
 #      0
-#  \end{bmatrix}
+#  \end{bmatrix}.
 # ```
 # Here $M$ is the singular block mass matrix, $K$ is the discretized Stokes operator and $N$ the nonlinear advection term, which
 # is also called trilinear form. $\hat{v}$ and $\hat{p}$ represent the time-dependent vectors of nodal values of the discretizations
 # of $v$ and $p$ respectively, while $\hat{\varphi}$ is the choice for the test function in the discretization. The hats are dropped
 # in the implementation and only stated for clarity in this section.
 #
+# ### [Derivation of the Jacobian](@id jacobian-derivation)
+# 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](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#stiff)
+# $J$ of $f(u, t) \equiv f(v, p, t)$ where $f(u, t)$ is just the RHS of the
+# mass matrix form given by $f(u, t) = Ku + N(u)$. By [definition](https://nl.mathworks.com/help/matlab/ref/odeset.html#d126e1203285),
+# ```math
+# J = \frac{\partial f}{\partial u} = \frac{\partial}{\partial u}(K u + N(u)) = K + \frac{\partial N}{\partial u}.
+# ```
+# It is simple to see that the Jacobian and the discretized Stokes operator $K$
+# share the same sparsity pattern since they share the same relation between
+# trial and test functions. This implies that once $K$ is assembled, its
+# 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
+# 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
+#
+# ```math
+#     - \sum_{i}^{N} (\int_{\Omega} ((\varphi_i \cdot \nabla) v + (v \cdot \nabla) \varphi_i) \cdot \varphi_j) \delta \hat{v}_i,
+# ```
+#
+# on a mesh with characteristic element size $h$ and $N$ trial functions.
+#
+# With the Jacobian accounting for both the linear and nonlinear contributions
+# of the semi-discrete weak form, we can easily use
+# [ODE solvers from DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).
+#
+# !!! details "Extra details on the Jacobian"
+#     Since the [directional derivative](https://en.wikipedia.org/wiki/Directional_derivative)
+#     is [equivalent to the Jacobian](https://math.stackexchange.com/questions/3191003/directional-derivative-and-jacobian-matrix),
+#     we take the directional derivative of $f(u, t) \equiv f(v, p, t)$ along
+#     $\delta u \equiv (\delta v, \delta p)$ (similar to [deal.ii: step 57](https://www.dealii.org/current/doxygen/deal.II/step_57.html))
+#     to explicitly determine the sparsity pattern
+#     and values for the Jacobian. Note that $f$ is a vector function with
+#     components $f_1$ and $f_2$ corresponding to
+#     the RHS of the first and second equations, respectively, of the
+#     weak form defined in the previous [section](@ref weak-form-derivation).
+#     By definition, and omitting $t$, the directional derivative is given by
+#
+#     ```math
+#     \begin{aligned}
+#       J &= \nabla f(v, p) \cdot (\delta v, \delta p) \\
+#       &= \lim_{\epsilon \to 0 } \frac{1}{\epsilon}(f(v + \epsilon \delta v, p + \epsilon \delta p ) - f(v, p)) \\
+#       &= \lim_{\epsilon \to 0}  \frac{1}{\epsilon} \left( \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} - \begin{bmatrix}
+#       f_1(v,p) \\
+#       f_2(v,p)
+#       \end{bmatrix} \right) \\
+#       &= \lim_{\epsilon \to 0} \frac{1}{\epsilon} \left( \begin{bmatrix}
+#       - \int_{\Omega} \cancel{\nu \nabla v : \nabla \varphi} + \nu \epsilon \nabla \delta v : \nabla \varphi - \int_{\Omega} (\cancel{(v \cdot \nabla) v} + (\epsilon \delta v \cdot \nabla) v + (v \cdot \nabla) \epsilon \delta v + (\epsilon \delta v \cdot \nabla) \epsilon \delta v) \cdot \varphi + \int_{\Omega} \cancel{p (\nabla \cdot \varphi)} + \epsilon \delta p (\nabla \cdot \varphi) \\
+#       \int_{\Omega} \cancel{(\nabla \cdot v)\psi} + (\nabla \cdot \epsilon \delta v)\psi
+#       \end{bmatrix} - \begin{bmatrix}
+#       \cancel{f_1(v,p)} \\
+#       \cancel{f_2(v,p)}
+#       \end{bmatrix} \right)\\
+#       &= \lim_{\epsilon \to 0} \frac{1}{\cancel{\epsilon}} \begin{bmatrix}
+#       - \int_{\Omega} \nu \cancel{\epsilon} \nabla \delta v : \nabla \varphi - \int_{\Omega} ((\cancel{\epsilon} \delta v \cdot \nabla) v + (v \cdot \nabla) \cancel{\epsilon} \delta v + \cancel{\epsilon^2 (\delta v \cdot \nabla) \delta v}) \cdot \varphi + \int_{\Omega} \cancel{\epsilon} \delta p (\nabla \cdot \varphi) \\
+#       \int_{\Omega} (\nabla \cdot \cancel{\epsilon} \delta v) \psi
+#       \end{bmatrix} \\
+#       &=  \begin{bmatrix}
+#       - \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{bmatrix}.
+#     \end{aligned}
+#     ```
+#
+#     The term in the Jacobian
+#
+#     ```math
+#     - \int_{\Omega} ((\delta v \cdot \nabla) v + (v \cdot \nabla) \delta v) \cdot \varphi,
+#     ```
+#
+#     corresponding to the nonlinear advection term is of particular interest since
+#     this integrand contains terms that evolve over time and thus these terms
+#     cannot be formulated as constant coefficients.
+#     To make this explicit, we construct a finite element approximation of
+#     $\delta v$ and $\delta p$ on a mesh with characteristic element size $h$
+#     and $N$ trial functions $\varphi$. An important note here is
+#     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,
+#     ```
+#     for a specific point $\mathbf{x}$ and nodal trial functions $\varphi_i$
+#     with unknown nodal values $\delta \hat{v}_i$. Substituting $\delta v_h$
+#     into the corresponding nonlinear advection term above, we have
+#
+#     ```math
+#     - \sum_{i}^{N} (\int_{\Omega} ((\varphi_i \cdot \nabla) v + (v \cdot \nabla) \varphi_i) \cdot \varphi_j) \delta \hat{v}_i,
+#     ```
+#     and we implement a function for the terms in this integrand.
+#     With this function for the nonlinear term and $K$, we can fully describe the
+#     Jacobian in the manner required
+#     by the [ODEFunction](https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction) used by the solver.
+#
+#     Note that the operation
+#
+#     ```math
+#     g(\xi, \theta) := (\xi \cdot \nabla) \theta
+#     ```
+#
+#     for vectors $\xi$ and $\theta$ can be implemented as either
+#
+#     ```math
+#     g(\xi, \theta) := \begin{cases}
+#       \xi \cdot (\nabla \theta)^{\text{T}},\ \text{or} \\
+#       \nabla \theta \cdot \xi.
+#     \end{cases}
+#     ```
+#
+#     We use the first of these two definitions in the implementation.
 #
 # ## Commented implementation
 #
@@ -363,23 +478,24 @@ K = assemble_stokes_matrix(cellvalues_v, cellvalues_p, ν, K, dh);
 u₀ = zeros(ndofs(dh))
 apply!(u₀, ch);
 
-# DifferentialEquations assumes dense matrices by default, which is not
-# feasible for semi-discretization of finite element models. We communicate
-# that a sparse matrix with specified pattern should be utilized through the
-# `jac_prototyp` argument. It is 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.
+# When using ODE solvers from DifferentialEquations for stiff problems, the
+# solvers assume the Jacobian is a dense matrix by default. This is not
+# feasible for the semi-discretization of finite element models, and it is much
+# more efficient to provide the sparsity pattern and construction of the Jacobian.
+# From the derivation of the jacobian [section](@ref jacobian-derivation),
+# we communicate that a sparse matrix with a specified pattern should be utilized
+# through the `jac_prototype` argument.
 jac_sparsity = sparse(K);
 
 # To apply the nonlinear portion of the Navier-Stokes problem we simply hand
-# over the dof handler and cell values to the right-hand-side (RHS) as a parameter.
+# over the dof handler and cell values to the RHS as a parameter.
 # Furthermore the pre-assembled linear part, our Stokes operator (which is time independent)
 # is passed to save some additional runtime. To apply the time-dependent Dirichlet BCs, we
 # also need to hand over the constraint handler.
 # The basic idea to apply the Dirichlet BCs consistently is that we copy the
-# current solution `u`, apply the Dirichlet BCs on the copy, evaluate the
+# current solution `u`, apply the Dirichlet BCs on the copy, and evaluate the
 # discretized RHS of the Navier-Stokes equations with this vector.
-# Furthermore we pass down the Jacobian assembly manually. For the Jacobian we eliminate all
+# Furthermore, we pass down the Jacobian assembly manually. For the Jacobian we eliminate all
 # rows and columns associated with constrained dofs. Also note that we eliminate the mass
 # matrix beforehand in a similar fashion. This decouples the time evolution of the constrained
 # dofs from the true unknowns. The correct solution is enforced by utilizing step and

docs/src/tutorials/index.md

@@ -145,7 +145,7 @@ for the time-integration.
 
 ---
 
-#### [Tutorial 10: Reactive surface](@ref tutorial-reactive-surface)
+#### [Tutorial 11: Reactive surface](@ref tutorial-reactive-surface)
 
 In this tutorial a reaction diffusion system on a sphere surface embedded in 3D is solved.
 Ferrite is used to assemble the diffusion operators and the mass matrices. The problem is
@@ -155,7 +155,7 @@ solved by using the usual first order reaction diffusion operator splitting.
 
 ---
 
-#### [Tutorial 11: Linear shell](@ref tutorial-linear-shell)
+#### [Tutorial 12: Linear shell](@ref tutorial-linear-shell)
 
 In this tutorial a linear shell element formulation is set up as a two-dimensional domain
 embedded in three-dimensional space. This will teach, and perhaps inspire, you on how
@@ -166,7 +166,7 @@ Ferrite.
 
 ---
 
-#### [Tutorial 12: Discontinuous Galerkin heat equation](@ref tutorial-dg-heat-equation)
+#### [Tutorial 13: Discontinuous Galerkin heat equation](@ref tutorial-dg-heat-equation)
 
 This tutorial guides you through the process of solving the linear stationary heat equation
 (i.e. Poisson's equation) on a unit square with inhomogeneous Dirichlet and Neumann boundary