From 4502c48cc983b02010297472bf3f9b7ee393a658 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Sun, 2 Feb 2025 15:11:10 +0100 Subject: [PATCH 01/14] Clarify construction of Jacobian in Navier-Stokes tutorial 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. --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 60 +++++++++++++++++++-- 1 file changed, 55 insertions(+), 5 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 67871b872e..03f5b96c0c 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -364,11 +364,61 @@ 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. +# feasible for the semi-discretization of finite element models. Since we use +# an implicit method to solve the ODEs, it is more efficient to provide +# information to the solver about the sparsity pattern and values of the Jacobian. +# We communicate that a sparse matrix with specified pattern should be utilized +# 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} +# ``` +# 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. jac_sparsity = sparse(K); # To apply the nonlinear portion of the Navier-Stokes problem we simply hand From 0f07aa3778e6de5aba19a93d655913e83e77795d Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Mon, 10 Feb 2025 19:24:38 +0100 Subject: [PATCH 02/14] brief outline for jacobian derivation reorganization --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 03f5b96c0c..f5f3e31621 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -113,6 +113,14 @@ nothing #hide # in the implementation and only stated for clarity in this section. # # +# ### Derivation of the Jacobian +# * 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) +# * The derivation of the jacobian (deal.ii) +# * Relevance to the discrete form (Whitely discrete, Bubnov-Galerkin approach and so $K^{ele}$ same for $J^{ele}$ for linear elements) +# * Add references section +# * Further reading +# # ## Commented implementation # # Now we solve the problem with Ferrite and [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). What follows is a program spliced with comments. From ef28beabf9b04e98d460e73736314119a539cd35 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Mon, 10 Feb 2025 20:32:48 +0100 Subject: [PATCH 03/14] added references section for clarity --- docs/src/assets/references.bib | 20 ++++++++++++++++++++ docs/src/literate-tutorials/ns_vs_diffeq.jl | 7 ++++++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index c240c43833..5d05815b74 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -190,3 +190,23 @@ @article{CroRav:1973:cnf year={1973}, publisher={EDP Sciences} } +@inbook{Bueler2021_Ch4NonlinearEquationsByNewtonsMethod, +author = {Bueler, Ed}, +title = {Chapter 4: Nonlinear equations by Newton's method}, +url = {https://epubs.siam.org/doi/abs/10.1137/1.9781611976311.ch4}, +booktitle = {PETSc for Partial Differential Equations: Numerical Solutions in C and Python}, +year = {2021}, +pages = {67-93}, +publisher = {Society for Industrial and Applied Mathematics} +} +@inbook{Whiteley2017_Ch5NonlinearBoundaryValueProblems, +title = {Nonlinear {Boundary} {Value} {Problems}}, +isbn = {978-3-319-49971-0}, +url = {https://doi.org/10.1007/978-3-319-49971-0_5}, +booktitle = {Finite {Element} {Methods}: {A} {Practical} {Guide}}, +publisher = {Springer International Publishing}, +author = {Whiteley, Jonathan}, +year = {2017}, +pages = {81-101}, +} + diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index f5f3e31621..459573db60 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -377,7 +377,7 @@ apply!(u₀, ch); # information to the solver about the sparsity pattern and values of the Jacobian. # We communicate that a sparse matrix with specified pattern should be utilized # 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 +# Jacobian can be shown by taking the directional gradient 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 @@ -701,3 +701,8 @@ end #md # ```julia #md # @__CODE__ #md # ``` + +# ## Further Reading +# * [deal.ii: step 57](https://www.dealii.org/current/doxygen/deal.II/step_57.html) +# * Chapter 4: Nonlinear Equations by Newton's Method in Bueler (2021) [Bueler2021_Ch4NonlinearEquationsByNewtonsMethod](@cite) +# * Nonlinear Boundary Value Problems in Whiteley (2017) [Whiteley2017_Ch5NonlinearBoundaryValueProblems](@cite) From c261e184f25bcdbdb7514a4a2a14f36f818c8830 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Sun, 16 Feb 2025 23:41:02 +0100 Subject: [PATCH 04/14] Fixing notation and rephrasing for jacobian derivation * 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`. --- docs/src/assets/references.bib | 9 -- docs/src/literate-tutorials/ns_vs_diffeq.jl | 163 +++++++++++--------- 2 files changed, 93 insertions(+), 79 deletions(-) diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index 5d05815b74..7a69769479 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -190,15 +190,6 @@ @article{CroRav:1973:cnf year={1973}, publisher={EDP Sciences} } -@inbook{Bueler2021_Ch4NonlinearEquationsByNewtonsMethod, -author = {Bueler, Ed}, -title = {Chapter 4: Nonlinear equations by Newton's method}, -url = {https://epubs.siam.org/doi/abs/10.1137/1.9781611976311.ch4}, -booktitle = {PETSc for Partial Differential Equations: Numerical Solutions in C and Python}, -year = {2021}, -pages = {67-93}, -publisher = {Society for Industrial and Applied Mathematics} -} @inbook{Whiteley2017_Ch5NonlinearBoundaryValueProblems, title = {Nonlinear {Boundary} {Value} {Problems}}, isbn = {978-3-319-49971-0}, diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 459573db60..f2bbf0a6fd 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -67,7 +67,7 @@ nothing #hide # $\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 @@ -112,14 +112,87 @@ nothing #hide # 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 +# 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 +# $f(u, t) \equiv f(v, p, t)$ along $\delta u \equiv (\delta v, \delta p)$. It is useful +# to frame the Jacobian of $f(u, t)$ in the context of the weak form described in +# the previous [section](@ref weak-form-derivation). The motivation for this +# framing will later become clear when determining the values and sparsity +# pattern of the Jacobian. By definition, and omitting $t$, the directional +# gradient is given by +# +# ```math +# \begin{aligned} +# \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} \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) \\ +# &= \lim_{\epsilon \to 0} \frac{1}{\epsilon} \begin{bmatrix} +# - \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{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} +# ``` # -# ### Derivation of the Jacobian -# * 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) -# * The derivation of the jacobian (deal.ii) -# * Relevance to the discrete form (Whitely discrete, Bubnov-Galerkin approach and so $K^{ele}$ same for $J^{ele}$ for linear elements) -# * Add references section -# * Further reading +# To determine the values and sparsity pattern of the Jacobian given +# $\nabla f(v, p) \cdot (\delta v, \delta p)$, we now consider the +# finite element approximation of $\delta v$ and $\delta p$ +# ```math +# \begin{aligned} +# \delta v_h(\mathbf{x}) &= \sum_{i}^{N} \varphi_i(\mathbf{x}) \delta \hat{v}_i \\ +# \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 +# $\varphi$ and $\psi$. An important note here is that the trial functions +# are the same as in the previous [section](@ref weak-form-derivation) where +# the finite element approximation of $v$ and $p$ is given by +# ```math +# \begin{aligned} +# v_h(\mathbf{x}) &= \sum_{i}^{N} \varphi_i(\mathbf{x}) \hat{v}_i \\ +# p_h(\mathbf{x}) &= \sum_{i}^{N} \psi_i(\mathbf{x}) \hat{p}_i. +# \end{aligned} +# ``` +# 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 +# incompressibility terms. This implies that once $K$ is assembled, its +# sparsity values and pattern can be used to initialize the Jacobian. +# In other words, it is simple to see that the Jacobian and the discretized +# Stokes operator share the same sparsity pattern since they share the same +# relation between trial and test functions. +# +# However, the Jacobian cannot be assembled using the values and pattern from $K$ +# alone since the nonlinear advection term is not accounted for by $K$. +# Substituting the finite element approximation into the term resulting from the +# directional gradient of the nonlinear advection term, we have +# ```math +# - \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi \xrightarrow{\delta v_h} - \sum_{i}^{N} (\int_{\Omega} (\varphi_i \cdot \nabla v + v \cdot \nabla \varphi_i) \cdot \varphi_j) \delta \hat{v}_i. +# ``` +# Since it is clear that the values of the Jacobian corresponding to the +# directional gradient of the nonlinear advection term rely on the velocity +# values $v$ and velocity gradient $\nabla v$ that change during the solving +# phase, the finite element contributions to the Jacobian for this term are +# calculated separately. + # # ## Commented implementation # @@ -191,7 +264,7 @@ gmsh.model.mesh.generate(dim) grid = togrid() Gmsh.finalize(); -# ### Function Space +# ### Function space # To ensure stability we utilize the Taylor-Hood element pair Q2-Q1. # We have to utilize the same quadrature rule for the pressure as for the velocity, because in the weak form the # linear pressure term is tested against a quadratic function. @@ -247,7 +320,7 @@ add!(ch, inflow_bc); close!(ch) update!(ch, 0.0); -# ### Linear System Assembly +# ### Linear system assembly # Next we describe how the block mass matrix and the Stokes matrix are assembled. # # For the block mass matrix $M$ we remember that only the first equation had a time derivative @@ -371,62 +444,13 @@ 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 the semi-discretization of finite element models. Since we use -# an implicit method to solve the ODEs, it is more efficient to provide -# information to the solver about the sparsity pattern and values of the Jacobian. -# We communicate that a sparse matrix with specified pattern should be utilized -# through the `jac_prototype` argument. The sparsity pattern and values for the -# Jacobian can be shown by taking the directional gradient 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} -# ``` -# 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. +# 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 @@ -435,9 +459,9 @@ jac_sparsity = sparse(K); # 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 @@ -702,7 +726,6 @@ end #md # @__CODE__ #md # ``` -# ## Further Reading +# ## Further reading # * [deal.ii: step 57](https://www.dealii.org/current/doxygen/deal.II/step_57.html) -# * Chapter 4: Nonlinear Equations by Newton's Method in Bueler (2021) [Bueler2021_Ch4NonlinearEquationsByNewtonsMethod](@cite) # * Nonlinear Boundary Value Problems in Whiteley (2017) [Whiteley2017_Ch5NonlinearBoundaryValueProblems](@cite) From 06fa3995c635ae705685cb1cdde5b78654b3fec3 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Tue, 18 Feb 2025 16:24:48 +0100 Subject: [PATCH 05/14] say directional derivative instead of directioal gradient! --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index f2bbf0a6fd..8db9cc674f 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -116,7 +116,7 @@ nothing #hide # 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 +# taking the [directional derivative](https://en.wikipedia.org/wiki/Directional_derivative) of # $f(u, t) \equiv f(v, p, t)$ along $\delta u \equiv (\delta v, \delta p)$. It is useful # to frame the Jacobian of $f(u, t)$ in the context of the weak form described in # the previous [section](@ref weak-form-derivation). The motivation for this @@ -164,7 +164,7 @@ nothing #hide # 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, +# with the analagous term in the directional derivative, # ```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 \\ @@ -183,12 +183,12 @@ nothing #hide # However, the Jacobian cannot be assembled using the values and pattern from $K$ # alone since the nonlinear advection term is not accounted for by $K$. # Substituting the finite element approximation into the term resulting from the -# directional gradient of the nonlinear advection term, we have +# directional derivative of the nonlinear advection term, we have # ```math # - \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi \xrightarrow{\delta v_h} - \sum_{i}^{N} (\int_{\Omega} (\varphi_i \cdot \nabla v + v \cdot \nabla \varphi_i) \cdot \varphi_j) \delta \hat{v}_i. # ``` # Since it is clear that the values of the Jacobian corresponding to the -# directional gradient of the nonlinear advection term rely on the velocity +# directional derivative of the nonlinear advection term rely on the velocity # values $v$ and velocity gradient $\nabla v$ that change during the solving # phase, the finite element contributions to the Jacobian for this term are # calculated separately. From 021bb32677c28fefac52ead9ed8f05946f7a9843 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Sun, 2 Mar 2025 17:28:17 +0100 Subject: [PATCH 06/14] init reformatting with extra details, TODO: come back to this --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 84 +++++++++++++-------- 1 file changed, 53 insertions(+), 31 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 8db9cc674f..534a152e0d 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -115,32 +115,56 @@ nothing #hide # ### [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 -# of $f(u, t)$. The sparsity pattern and values for the Jacobian can be shown by -# taking the [directional derivative](https://en.wikipedia.org/wiki/Directional_derivative) of -# $f(u, t) \equiv f(v, p, t)$ along $\delta u \equiv (\delta v, \delta p)$. It is useful -# to frame the Jacobian of $f(u, t)$ in the context of the weak form described in -# the previous [section](@ref weak-form-derivation). The motivation for this -# framing will later become clear when determining the values and sparsity -# pattern of the Jacobian. By definition, and omitting $t$, the directional -# gradient is given by -# +# $J$ of $f(u, t) \equiv f(v, p, t)$. By definition, # ```math -# \begin{aligned} -# \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} \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) \\ -# &= \lim_{\epsilon \to 0} \frac{1}{\epsilon} \begin{bmatrix} -# - \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{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} -# ``` +# 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. +# +# TODO: Add $u$ to semi-discrete weak form as vector of unknowns syntax +# +# !!! 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 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 right handside 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} (\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} +# ``` # # To determine the values and sparsity pattern of the Jacobian given # $\nabla f(v, p) \cdot (\delta v, \delta p)$, we now consider the @@ -151,7 +175,7 @@ nothing #hide # \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 +# on a mesh with characteristic element size $h$ and $N$ trial (aka nodal shape) functions # $\varphi$ and $\psi$. An important note here is that the trial functions # are the same as in the previous [section](@ref weak-form-derivation) where # the finite element approximation of $v$ and $p$ is given by @@ -164,7 +188,7 @@ nothing #hide # 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 derivative, +# with the analogous term in the directional derivative, # ```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 \\ @@ -176,9 +200,6 @@ nothing #hide # are the same. The same can be shown for the pressure term and the # incompressibility terms. This implies that once $K$ is assembled, its # sparsity values and pattern can be used to initialize the Jacobian. -# In other words, it is simple to see that the Jacobian and the discretized -# Stokes operator share the same sparsity pattern since they share the same -# relation between trial and test functions. # # However, the Jacobian cannot be assembled using the values and pattern from $K$ # alone since the nonlinear advection term is not accounted for by $K$. @@ -727,5 +748,6 @@ end #md # ``` # ## Further reading -# * [deal.ii: step 57](https://www.dealii.org/current/doxygen/deal.II/step_57.html) +# TODO: Move refs (also odeset ref) to details +# * [deal.ii: step 57]( # * Nonlinear Boundary Value Problems in Whiteley (2017) [Whiteley2017_Ch5NonlinearBoundaryValueProblems](@cite) From db39dbbfadcad6c2a48b75ed63f09856d166c0b2 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Tue, 18 Mar 2025 19:02:46 +0100 Subject: [PATCH 07/14] WIP: more structure to extra details section and TODOs added --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 114 +++++++++++--------- 1 file changed, 64 insertions(+), 50 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 534a152e0d..853b0a60d8 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -113,41 +113,44 @@ nothing #hide # 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 -# $J$ of $f(u, t) \equiv f(v, p, t)$. By definition, +# 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 right hand side 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. +# ``` +# 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. # -# TODO: Add $u$ to semi-discrete weak form as vector of unknowns syntax +# TODO: Some brief summary here about how dN/du is computed, then go into extra +# details... # -# !!! 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 determine the sparsity pattern -# and values for the Jacobian. Note that $f$ is a vector function with +# !!! 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 right handside of the first and second equations, respectively, of the +# the right handside 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 -# +# 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 +# \int_{\Omega} (\nabla \cdot (v + \epsilon \delta v)) \psi # \end{bmatrix} - \begin{bmatrix} # f_1(v,p) \\ # f_2(v,p) -# \end{bmatrix} \right) \\ +# \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} (\nabla \cdot v)\psi + (\nabla \cdot \epsilon \delta v)\psi @@ -156,18 +159,35 @@ nothing #hide # \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} \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} \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} # ``` # -# To determine the values and sparsity pattern of the Jacobian given -# $\nabla f(v, p) \cdot (\delta v, \delta p)$, we now consider the +# 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$ and $\psi$. An important note here is +# that the trial functions are the same as in the previous [section](@ref weak-form-derivation). +# +# TODO: Brief insert here about finite element discretization since this +# will show how dN/du is implemented in the code... +# +# To determine the values and sparsity pattern of the Jacobian given +# $\nabla f(v, p) \cdot (\delta v, \delta p)$, we now consider the # finite element approximation of $\delta v$ and $\delta p$ # ```math # \begin{aligned} @@ -175,8 +195,8 @@ nothing #hide # \delta p_h(\mathbf{x}) &= \sum_{i}^{N} \psi_i(\mathbf{x}) \delta \hat{p}_i. # \end{aligned} # ``` -# on a mesh with characteristic element size $h$ and $N$ trial (aka nodal shape) functions -# $\varphi$ and $\psi$. An important note here is that the trial functions +# on a mesh with characteristic element size $h$ and $N$ trial (aka nodal shape) functions +# $\varphi$ and $\psi$. An important note here is that the trial functions # are the same as in the previous [section](@ref weak-form-derivation) where # the finite element approximation of $v$ and $p$ is given by # ```math @@ -185,33 +205,32 @@ nothing #hide # p_h(\mathbf{x}) &= \sum_{i}^{N} \psi_i(\mathbf{x}) \hat{p}_i. # \end{aligned} # ``` -# 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 +# 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 analogous term in the directional derivative, # ```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 \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} +# \end{aligned} # ``` -# we see that the coefficients (i.e., the terms in the integrand on the right +# 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 -# incompressibility terms. This implies that once $K$ is assembled, its -# sparsity values and pattern can be used to initialize the Jacobian. +# incompressibility terms. # -# However, the Jacobian cannot be assembled using the values and pattern from $K$ -# alone since the nonlinear advection term is not accounted for by $K$. -# Substituting the finite element approximation into the term resulting from the +# However, the Jacobian cannot be assembled using the values and pattern from $K$ +# alone since the nonlinear advection term is not accounted for by $K$. +# Substituting the finite element approximation into the term resulting from the # directional derivative of the nonlinear advection term, we have # ```math # - \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi \xrightarrow{\delta v_h} - \sum_{i}^{N} (\int_{\Omega} (\varphi_i \cdot \nabla v + v \cdot \nabla \varphi_i) \cdot \varphi_j) \delta \hat{v}_i. # ``` -# Since it is clear that the values of the Jacobian corresponding to the +# Since it is clear that the values of the Jacobian corresponding to the # directional derivative of the nonlinear advection term rely on the velocity -# values $v$ and velocity gradient $\nabla v$ that change during the solving -# phase, the finite element contributions to the Jacobian for this term are +# values $v$ and velocity gradient $\nabla v$ that change during the solving +# phase, the finite element contributions to the Jacobian for this term are # calculated separately. # @@ -465,13 +484,13 @@ K = assemble_stokes_matrix(cellvalues_v, cellvalues_p, ν, K, dh); u₀ = zeros(ndofs(dh)) apply!(u₀, ch); -# When using ODE solvers from DifferentialEquations for stiff problems, the -# solvers assume the Jacobian is a dense matrix by default. This is not +# 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. +# 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. +# 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 @@ -746,8 +765,3 @@ end #md # ```julia #md # @__CODE__ #md # ``` - -# ## Further reading -# TODO: Move refs (also odeset ref) to details -# * [deal.ii: step 57]( -# * Nonlinear Boundary Value Problems in Whiteley (2017) [Whiteley2017_Ch5NonlinearBoundaryValueProblems](@cite) From ef6f28560b2d77d69abc1efec29861b86fc81940 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Sat, 29 Mar 2025 11:10:24 +0100 Subject: [PATCH 08/14] extra details and clarifications on jacobian derivation Also added punctuation to maths equations in `ns_vs_diffeq.jl` Updated `tutorials/index.md` with the correct numbering of tutorials. --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 103 ++++++++------------ docs/src/tutorials/index.md | 6 +- 2 files changed, 42 insertions(+), 67 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 853b0a60d8..e62ac9cad1 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/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,7 +59,7 @@ 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 @@ -74,8 +74,8 @@ nothing #hide # 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,7 +105,7 @@ 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 @@ -115,7 +115,7 @@ nothing #hide # ### [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 right hand side of the +# $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}. @@ -125,8 +125,20 @@ nothing #hide # trial and test functions. This implies that once $K$ is assembled, its # sparsity values and pattern can be used to initialize the Jacobian. # -# TODO: Some brief summary here about how dN/du is computed, then go into extra -# details... +# 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=1}^{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) @@ -136,7 +148,7 @@ nothing #hide # 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 right handside of the first and second equations, respectively, of the +# 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 # @@ -153,7 +165,7 @@ nothing #hide # \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} (\nabla \cdot v)\psi + (\nabla \cdot \epsilon \delta v)\psi +# \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)} @@ -180,59 +192,22 @@ nothing #hide # 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$ and $\psi$. An important note here is -# that the trial functions are the same as in the previous [section](@ref weak-form-derivation). -# -# TODO: Brief insert here about finite element discretization since this -# will show how dN/du is implemented in the code... +# 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 # -# To determine the values and sparsity pattern of the Jacobian given -# $\nabla f(v, p) \cdot (\delta v, \delta p)$, we now consider the -# finite element approximation of $\delta v$ and $\delta p$ -# ```math -# \begin{aligned} -# \delta v_h(\mathbf{x}) &= \sum_{i}^{N} \varphi_i(\mathbf{x}) \delta \hat{v}_i \\ -# \delta p_h(\mathbf{x}) &= \sum_{i}^{N} \psi_i(\mathbf{x}) \delta \hat{p}_i. -# \end{aligned} -# ``` -# on a mesh with characteristic element size $h$ and $N$ trial (aka nodal shape) functions -# $\varphi$ and $\psi$. An important note here is that the trial functions -# are the same as in the previous [section](@ref weak-form-derivation) where -# the finite element approximation of $v$ and $p$ is given by -# ```math -# \begin{aligned} -# v_h(\mathbf{x}) &= \sum_{i}^{N} \varphi_i(\mathbf{x}) \hat{v}_i \\ -# p_h(\mathbf{x}) &= \sum_{i}^{N} \psi_i(\mathbf{x}) \hat{p}_i. -# \end{aligned} -# ``` -# 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 analogous term in the directional derivative, -# ```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 -# incompressibility terms. +# ```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 # -# However, the Jacobian cannot be assembled using the values and pattern from $K$ -# alone since the nonlinear advection term is not accounted for by $K$. -# Substituting the finite element approximation into the term resulting from the -# directional derivative of the nonlinear advection term, we have -# ```math -# - \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi \xrightarrow{\delta v_h} - \sum_{i}^{N} (\int_{\Omega} (\varphi_i \cdot \nabla v + v \cdot \nabla \varphi_i) \cdot \varphi_j) \delta \hat{v}_i. -# ``` -# Since it is clear that the values of the Jacobian corresponding to the -# directional derivative of the nonlinear advection term rely on the velocity -# values $v$ and velocity gradient $\nabla v$ that change during the solving -# phase, the finite element contributions to the Jacobian for this term are -# calculated separately. - +# ```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 the integrand. With this +# function 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. # # ## Commented implementation # @@ -494,7 +469,7 @@ apply!(u₀, ch); 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. diff --git a/docs/src/tutorials/index.md b/docs/src/tutorials/index.md index 3fe4bda264..048ac3a3f8 100644 --- a/docs/src/tutorials/index.md +++ b/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 From 5a589f28977a0d5bf2f8689635fcdc8727804d12 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Sat, 29 Mar 2025 11:21:44 +0100 Subject: [PATCH 09/14] removed Whitely practical FEM ref since not used --- docs/src/assets/references.bib | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index 7a69769479..c240c43833 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -190,14 +190,3 @@ @article{CroRav:1973:cnf year={1973}, publisher={EDP Sciences} } -@inbook{Whiteley2017_Ch5NonlinearBoundaryValueProblems, -title = {Nonlinear {Boundary} {Value} {Problems}}, -isbn = {978-3-319-49971-0}, -url = {https://doi.org/10.1007/978-3-319-49971-0_5}, -booktitle = {Finite {Element} {Methods}: {A} {Practical} {Guide}}, -publisher = {Springer International Publishing}, -author = {Whiteley, Jonathan}, -year = {2017}, -pages = {81-101}, -} - From f1301f96b5975db2c6cd307dbc8f673c0431ab31 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Wed, 2 Apr 2025 00:47:41 +0200 Subject: [PATCH 10/14] nonlinear advection Jacobian term needed negative sign in derivation! --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index e62ac9cad1..0ea2bd0dfc 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -184,7 +184,7 @@ nothing #hide # The term in the Jacobian # # ```math -# \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi, +# - \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 @@ -203,7 +203,7 @@ nothing #hide # 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, +# - \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 the integrand. With this # function and $K$, we can fully describe the Jacobian in the manner required From 0f6aeda4d2acb7d38246d77b549a4e19e5157572 Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Sun, 13 Apr 2025 10:41:26 +0200 Subject: [PATCH 11/14] corrected dot(v, gradv) to dot(v, grad)v in Jacobian derivation; note on implementing this operation --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 35 ++++++++++++++++----- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 0ea2bd0dfc..8944102a35 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -132,8 +132,9 @@ nothing #hide # $\frac{\delta N}{\delta u}$ is therefore given by # # ```math -# - \sum_{i=1}^{N} (\int_{\Omega} (\varphi_i \cdot \nabla v + v \cdot \nabla \varphi_i) \cdot \varphi_j) \delta \hat{v}_i, +# - \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 @@ -164,18 +165,18 @@ nothing #hide # 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{\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} \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} \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} @@ -184,7 +185,7 @@ nothing #hide # The term in the Jacobian # # ```math -# - \int_{\Omega} (\delta v \cdot \nabla v + v \cdot \nabla \delta v) \cdot \varphi, +# - \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 @@ -203,12 +204,30 @@ nothing #hide # 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, +# - \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 the integrand. With this -# function and $K$, we can fully describe the Jacobian in the manner required +# 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 # # Now we solve the problem with Ferrite and [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). What follows is a program spliced with comments. From 0376d8f0f5abe11a026dbc60264a6b2974b6c280 Mon Sep 17 00:00:00 2001 From: Jared Frazier <72281389+jfdev001@users.noreply.github.com> Date: Mon, 14 Apr 2025 16:58:15 +0200 Subject: [PATCH 12/14] WIP: Apply suggestions from code review * 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 --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index 8944102a35..e29f8d739a 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -122,17 +122,17 @@ nothing #hide # ``` # 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 +# trial and test functions, ignoring the zero values in the form. 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 +# To derive $\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 +# $\frac{\partial N}{\partial 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, +# \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, # ``` # # on a mesh with characteristic element size $h$ and $N$ trial functions. From c1d225ff5008735cc57730a5b3fbf5194f77d95e Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Mon, 14 Apr 2025 17:37:53 +0200 Subject: [PATCH 13/14] replace N with n for n trial functions --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index e29f8d739a..ecda5f9457 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -132,10 +132,10 @@ nothing #hide # $\frac{\partial N}{\partial u}$ is therefore given by # # ```math -# \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, +# \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, # ``` # -# on a mesh with characteristic element size $h$ and $N$ trial functions. +# for $n$ trial functions. # # With the Jacobian accounting for both the linear and nonlinear contributions # of the semi-discrete weak form, we can easily use @@ -193,18 +193,19 @@ nothing #hide # 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 +# 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, +# \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, +# - \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 From 54fbfbd003775934aec99babc4983bdaf92cfbed Mon Sep 17 00:00:00 2001 From: Jared Frazier Date: Wed, 16 Apr 2025 08:26:36 +0200 Subject: [PATCH 14/14] WIP: Towards use of discrete notation in derivation of Jacobian --- docs/src/literate-tutorials/ns_vs_diffeq.jl | 75 ++++++++++++--------- 1 file changed, 43 insertions(+), 32 deletions(-) diff --git a/docs/src/literate-tutorials/ns_vs_diffeq.jl b/docs/src/literate-tutorials/ns_vs_diffeq.jl index ecda5f9457..43249ca18d 100644 --- a/docs/src/literate-tutorials/ns_vs_diffeq.jl +++ b/docs/src/literate-tutorials/ns_vs_diffeq.jl @@ -97,42 +97,48 @@ nothing #hide # A & B^{\textrm{T}} \\ # B & 0 # \end{bmatrix}}_{:=K} -# \begin{bmatrix} +# \underbrace{\begin{bmatrix} # \hat{v} \\ # \hat{p} -# \end{bmatrix} +# \end{bmatrix}}_{:=\hat{u}} # + # \begin{bmatrix} # N(\hat{v}, \hat{v}, \hat{\varphi}) \\ # 0 # \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. +# Here $M$ is the singular block mass matrix, $K$ is the discretized Stokes +# operator and $N$ the nonlinear advection term, which is also called the trilinear +# form. The solution vector $\hat{u}$ has components $\hat{v}$ and $\hat{p}$ +# that represent the time-dependent vectors of nodal values of the +# discretizations of $v$ and $p$, respectively. The discrete form of the +# test function $\varphi$ is denoted by $\hat{\varphi}$. +# 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), +# $J$ of $f(\hat{u}, t) \equiv f(\hat{v}, \hat{p}, t)$ where $f(\hat{u}, t)$ is +# just the discrete form of the RHS of the mass matrix form given by $f(\hat{u}, t) = K\hat{u} + N(\hat{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}. +# J = \frac{\partial f}{\partial \hat{u}} = \frac{\partial}{\partial \hat{u}}(K \hat{u} + N(\hat{u})) = K + \frac{\partial N}{\partial \hat{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, ignoring the zero values in the form. This implies that once $K$ is assembled, its # sparsity values and pattern can be used to initialize the Jacobian. # -# To derive $\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{\partial N}{\partial u}$ is therefore given by +# To derive $\frac{\partial N}{\partial \hat{u}}$, we take the directional derivative +# of $f$ along the vector $(\delta \hat{v}, \delta \hat{p})$ where +# $\delta \hat{v}$ and $\delta \hat{p}$ belong to finite dimensional subspaces +# of the corresponding function spaces for $\delta v$ and $\delta p$. Naturally, +# the trial functions $\hat{\varphi}$ and $\hat{\psi}$ are used for +# $\delta \hat{v}$ and $\delta \hat{p}$, respectively. The discrete form of +# $\frac{\partial N}{\partial \hat{u}}$ is therefore given by # # ```math -# \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, +# \frac{\partial N}{\partial \hat{u}} = - \sum_{i}^{n} \left(\int_{\Omega} ((\hat{\varphi_i} \cdot \nabla) \hat{v} + (\hat{v} \cdot \nabla) \hat{\varphi_i}) \cdot \hat{\varphi_j} \cdot \delta \hat{v}_i\right), # ``` # # for $n$ trial functions. @@ -144,8 +150,9 @@ nothing #hide # !!! 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)) +# we take the directional derivative of +# $f(\hat{u}, t) \equiv f(\hat{v}, \hat{p}, t)$ along +# $\delta \hat{u} \equiv (\delta \hat{v}, \delta \hat{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 @@ -155,29 +162,29 @@ nothing #hide # # ```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)) \\ +# J &= \nabla f(\hat{v}, \hat{p}) \cdot (\delta \hat{v}, \delta \hat{p}) \\ +# &= \lim_{\epsilon \to 0 } \frac{1}{\epsilon}(f(\hat{v} + \epsilon \delta \hat{v}, \hat{p} + \epsilon \delta \hat{p} ) - f(\hat{v}, \hat{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 +# - \int_\Omega \nu \nabla (\hat{v} + \epsilon \delta \hat{v}) : \nabla \hat{\varphi} - \int_\Omega ((\hat{v} + \epsilon \delta \hat{v}) \cdot \nabla) (\hat{v} + \epsilon \delta \hat{v}) \cdot \hat{\varphi} + \int_\Omega (\hat{p} + \epsilon \delta \hat{p})(\nabla \cdot \hat{\varphi}) \\ +# \int_{\Omega} (\nabla \cdot (\hat{v} + \epsilon \delta \hat{v})) \hat{\psi} # \end{bmatrix} - \begin{bmatrix} -# f_1(v,p) \\ -# f_2(v,p) +# f_1(\hat{v},\hat{p}) \\ +# f_2(\hat{v},\hat{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 +# - \int_{\Omega} \cancel{\nu \nabla \hat{v} : \nabla \hat{\varphi}} + \nu \epsilon \nabla \delta \hat{v} : \nabla \hat{\varphi} - \int_{\Omega} (\cancel{(\hat{v} \cdot \nabla) \hat{v}} + (\epsilon \delta \hat{v} \cdot \nabla) \hat{v} + (\hat{v} \cdot \nabla) \epsilon \delta \hat{v} + (\epsilon \delta \hat{v} \cdot \nabla) \epsilon \delta \hat{v}) \cdot \hat{\varphi} + \int_{\Omega} \cancel{\hat{p} (\nabla \cdot \hat{\varphi})} + \epsilon \delta \hat{p} (\nabla \cdot \hat{\varphi}) \\ +# \int_{\Omega} \cancel{(\nabla \cdot \hat{v})\hat{\psi}} + (\nabla \cdot \epsilon \delta \hat{v})\hat{\psi} # \end{bmatrix} - \begin{bmatrix} -# \cancel{f_1(v,p)} \\ -# \cancel{f_2(v,p)} +# \cancel{f_1(\hat{v},\hat{p})} \\ +# \cancel{f_2(\hat{v},\hat{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 +# - \int_{\Omega} \nu \cancel{\epsilon} \nabla \delta \hat{v} : \nabla \hat{\varphi} - \int_{\Omega} ((\cancel{\epsilon} \delta \hat{v} \cdot \nabla) \hat{v} + (\hat{v} \cdot \nabla) \cancel{\epsilon} \delta \hat{v} + \cancel{\epsilon^2 (\delta \hat{v} \cdot \nabla) \delta \hat{v}}) \cdot \hat{\varphi} + \int_{\Omega} \cancel{\epsilon} \delta \hat{p} (\nabla \cdot \hat{\varphi}) \\ +# \int_{\Omega} (\nabla \cdot \cancel{\epsilon} \delta \hat{v}) \hat{\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 +# - \int_{\Omega} \nu \nabla \delta \hat{v} : \nabla \hat{\varphi} - \int_{\Omega} ((\delta \hat{v} \cdot \nabla) \hat{v} + (\hat{v} \cdot \nabla) \delta \hat{v}) \cdot \hat{\varphi} + \int_{\Omega} \delta \hat{p} (\nabla \cdot \hat{\varphi}) \\ +# \int_{\Omega} (\nabla \cdot \delta \hat{v}) \hat{\psi} # \end{bmatrix}. # \end{aligned} # ``` @@ -185,9 +192,11 @@ nothing #hide # The term in the Jacobian # # ```math -# - \int_{\Omega} ((\delta v \cdot \nabla) v + (v \cdot \nabla) \delta v) \cdot \varphi, +# - \int_{\Omega} ((\delta \hat{v} \cdot \nabla) \hat{v} + (\hat{v} \cdot \nabla) \delta \hat{v}) \cdot \hat{\varphi} # ``` # +# BEGIN TODO: Use discrete notation for finite element stuff in the next section... +# # 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. @@ -212,6 +221,8 @@ nothing #hide # Jacobian in the manner required # by the [ODEFunction](https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction) used by the solver. # +# END TODO +# # Note that the operation # # ```math