Update time-dependent problems. (#280)

* Fix nitsche

* Update further tutorials. Add a small section on petsc LU factorization

Author: jorgensd

chapter1/fundamentals_code.ipynb

@@ -387,6 +387,7 @@
     "Using {py:meth}`problem.solve()<dolfinx.fem.petsc.LinearProblem.solve>` we solve the linear system of equations and\n",
     "return a {py:class}`Function<dolfinx.fem.Function>` containing the solution.\n",
     "\n",
+    "(error-norm)=\n",
     "## Computing the error\n",
     "Finally, we want to compute the error to check the accuracy of the solution.\n",
     "We do this by comparing the finite element solution `u` with the exact solution.\n",
@@ -417,7 +418,7 @@
     "only assembles over the cells on the local process.\n",
     "This means that if we use 2 processes to solve our problem,\n",
     "we need to accumulate the local contributions to get the global error (on one or all processes).\n",
-    "We can do this with the {func}`Comm.allreduce<mpi4py.MPI.Comm.allreduce>` function."
+    "We can do this with the {py:meth}`Comm.allreduce<mpi4py.MPI.Comm.allreduce>` function."
    ]
   },
   {
@@ -438,9 +439,10 @@
    "metadata": {},
    "source": [
     "Secondly, we compute the maximum error at any degree of freedom.\n",
-    "As the finite element function $u$ can be expressed as a linear combination of basis functions $\\phi_j$, spanning the space $V$:\n",
-    "$ u = \\sum_{j=1}^N U_j\\phi_j.$\n",
-    "By writing `problem.solve()` we compute all the coefficients $U_1,\\dots, U_N$.\n",
+    "As the finite element function $u$ can be expressed as a linear combination of basis functions $\\phi_j$,\n",
+    "spanning the space $V$: $ u = \\sum_{j=1}^N U_j\\phi_j.$\n",
+    "By writing {py:meth}`problem.solve()<dolfinx.fem.petsc.LinearProblem.sovle>`\n",
+    "we compute all the coefficients $U_1,\\dots, U_N$.\n",
     "These values are known as the _degrees of freedom_ (dofs).\n",
     "We can access the degrees of freedom by accessing the underlying vector in `uh`.\n",
     "However, as a second order function space has more dofs than a linear function space,\n",
@@ -457,8 +459,7 @@
    "outputs": [],
    "source": [
     "error_max = numpy.max(numpy.abs(uD.x.array - uh.x.array))\n",
-    "# Only print the error on one process\n",
-    "if domain.comm.rank == 0:\n",
+    "if domain.comm.rank == 0:  # Only print the error on one process\n",
     "    print(f\"Error_L2 : {error_L2:.2e}\")\n",
     "    print(f\"Error_max : {error_max:.2e}\")"
    ]

chapter1/fundamentals_code.py

@@ -265,6 +265,7 @@
 # Using {py:meth}`problem.solve()<dolfinx.fem.petsc.LinearProblem.solve>` we solve the linear system of equations and
 # return a {py:class}`Function<dolfinx.fem.Function>` containing the solution.
 #
+# (error-norm)=
 # ## Computing the error
 # Finally, we want to compute the error to check the accuracy of the solution.
 # We do this by comparing the finite element solution `u` with the exact solution.
@@ -281,16 +282,17 @@
 # only assembles over the cells on the local process.
 # This means that if we use 2 processes to solve our problem,
 # we need to accumulate the local contributions to get the global error (on one or all processes).
-# We can do this with the {func}`Comm.allreduce<mpi4py.MPI.Comm.allreduce>` function.
+# We can do this with the {py:meth}`Comm.allreduce<mpi4py.MPI.Comm.allreduce>` function.
 
 L2_error = fem.form(ufl.inner(uh - uex, uh - uex) * ufl.dx)
 error_local = fem.assemble_scalar(L2_error)
 error_L2 = numpy.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))
 
 # Secondly, we compute the maximum error at any degree of freedom.
-# As the finite element function $u$ can be expressed as a linear combination of basis functions $\phi_j$, spanning the space $V$:
-# $ u = \sum_{j=1}^N U_j\phi_j.$
-# By writing `problem.solve()` we compute all the coefficients $U_1,\dots, U_N$.
+# As the finite element function $u$ can be expressed as a linear combination of basis functions $\phi_j$,
+# spanning the space $V$: $ u = \sum_{j=1}^N U_j\phi_j.$
+# By writing {py:meth}`problem.solve()<dolfinx.fem.petsc.LinearProblem.sovle>`
+# we compute all the coefficients $U_1,\dots, U_N$.
 # These values are known as the _degrees of freedom_ (dofs).
 # We can access the degrees of freedom by accessing the underlying vector in `uh`.
 # However, as a second order function space has more dofs than a linear function space,
@@ -299,8 +301,7 @@
 # we can compare the maximum values at any degree of freedom of the approximation space.
 
 error_max = numpy.max(numpy.abs(uD.x.array - uh.x.array))
-# Only print the error on one process
-if domain.comm.rank == 0:
+if domain.comm.rank == 0:  # Only print the error on one process
     print(f"Error_L2 : {error_L2:.2e}")
     print(f"Error_max : {error_max:.2e}")
 

chapter1/membrane_code.py

@@ -246,7 +246,7 @@ def on_boundary(x):
 # This function returns a list of cells whose bounding box collide for each input point.
 # As different points might have different number of cells, the data is stored in
 # {py:class}`dolfinx.graph.AdjacencyList`, where one can access the cells for the
-# `i`th point by calling `links(i)<dolfinx.graph.AdjacencyList.links>`.
+# `i`th point by calling {py:meth}`links(i)<dolfinx.graph.AdjacencyList.links>`.
 # However, as the bounding box of a cell spans more of $\mathbb{R}^n$ than the actual cell,
 # we check that the actual cell collides with the input point using
 # {py:func}`dolfinx.geometry.compute_colliding_cells`,

chapter1/nitsche.ipynb

@@ -8,8 +8,10 @@
     "# Weak imposition of Dirichlet conditions for the Poisson problem\n",
     "Author: Jørgen S. Dokken\n",
     "\n",
-    "In this section, we will go through how to solve the Poisson problem from the [Fundamentals](./fundamentals_code.ipynb) tutorial using Nitsche's method {cite}`Nitsche1971`.\n",
-    "The idea of weak imposition is that we add additional terms to the variational formulation to impose the boundary condition, instead of modifying the matrix system using strong imposition (lifting).\n",
+    "In this section, we will go through how to solve the Poisson problem from the\n",
+    "[Fundamentals](./fundamentals_code.ipynb) tutorial using Nitsche's method {cite}`nitsche-Nitsche1971`.\n",
+    "The idea of weak imposition is that we add additional terms to the variational formulation to\n",
+    "impose the boundary condition, instead of modifying the matrix system using strong imposition (lifting).\n",
     "\n",
     "We start by importing the required modules and creating the mesh and function space for our solution"
    ]
@@ -48,7 +50,12 @@
    "id": "2",
    "metadata": {},
    "source": [
-    "Next, we create a function containing the exact solution (which will also be used in the Dirichlet boundary condition) and the corresponding source function for the right hand side. Note that we use `ufl.SpatialCoordinate` to define the exact solution, which in turn is interpolated into `uD` and used to create the source function `f`."
+    "Next, we create a function containing the exact solution (which will also be used in the Dirichlet boundary condition)\n",
+    "and the corresponding source function for the right hand side. Note that we use {py:class}`ufl.SpatialCoordinate`\n",
+    "to define the exact solution, which in turn is interpolated into {py:class}`uD<dolfinx.fem.Function>`\n",
+    "by wrapping the {py:class}`ufl-expression<ufl.core.expr.Expr>` as a {py:class}`dolfinx.fem.Expression`.\n",
+    "For the source function, we use the symbolic differentiation capabilities of UFL to\n",
+    "compute the negative Laplacian of the exact solution, which we can use directly in the variational formulation."
    ]
   },
   {
@@ -72,21 +79,34 @@
    "source": [
     "As opposed to the first tutorial, we now have to have another look at the variational form.\n",
     "We start by integrating the problem by parts, to obtain\n",
+    "\n",
+    "$$\n",
     "\\begin{align}\n",
-    "    \\int_{\\Omega} \\nabla u \\cdot \\nabla v~\\mathrm{d}x - \\int_{\\partial\\Omega}\\nabla u \\cdot n v~\\mathrm{d}s = \\int_{\\Omega} f v~\\mathrm{d}x.\n",
+    "  \\int_{\\Omega} \\nabla u \\cdot \\nabla v~\\mathrm{d}x\n",
+    "  - \\int_{\\partial\\Omega}\\nabla u \\cdot n v~\\mathrm{d}s = \\int_{\\Omega} f v~\\mathrm{d}x.\n",
     "\\end{align}\n",
+    "$$\n",
+    "\n",
     "As we are not using strong enforcement, we do not set the trace of the test function to $0$ on the outer boundary.\n",
     "Instead, we add the following two terms to the variational formulation\n",
+    "\n",
+    "$$\n",
     "\\begin{align}\n",
-    "    -\\int_{\\partial\\Omega} \\nabla  v \\cdot n (u-u_D)~\\mathrm{d}s + \\frac{\\alpha}{h} \\int_{\\partial\\Omega} (u-u_D)v~\\mathrm{d}s.\n",
+    "  -\\int_{\\partial\\Omega} \\nabla  v \\cdot n (u-u_D)~\\mathrm{d}s\n",
+    "  + \\frac{\\alpha}{h} \\int_{\\partial\\Omega} (u-u_D)v~\\mathrm{d}s.\n",
     "\\end{align}\n",
+    "$$\n",
+    "\n",
     "where the first term enforces symmetry to the bilinear form, while the latter term enforces coercivity.\n",
     "$u_D$ is the known Dirichlet condition, and $h$ is the diameter of the circumscribed sphere of the mesh element.\n",
     "We create bilinear and linear form, $a$ and $L$\n",
+    "\n",
+    "$$\n",
     "\\begin{align}\n",
     "    a(u, v) &= \\int_{\\Omega} \\nabla u \\cdot \\nabla v~\\mathrm{d}x + \\int_{\\partial\\Omega}-(n \\cdot\\nabla u) v - (n \\cdot \\nabla v) u + \\frac{\\alpha}{h} uv~\\mathrm{d}s,\\\\\n",
     "    L(v) &= \\int_{\\Omega} fv~\\mathrm{d}x + \\int_{\\partial\\Omega} -(n \\cdot \\nabla v) u_D + \\frac{\\alpha}{h} u_Dv~\\mathrm{d}s\n",
-    "\\end{align}"
+    "\\end{align}\n",
+    "$$"
    ]
   },
   {
@@ -112,7 +132,8 @@
    "id": "6",
    "metadata": {},
    "source": [
-    "As we now have the variational form, we can solve the linear problem"
+    "As we now have the variational form, we can solve the arising\n",
+    "{py:class}`linear problem<dolfinx.fem.petsc.LinearProblem>`"
    ]
   },
   {
@@ -176,7 +197,9 @@
    "id": "12",
    "metadata": {},
    "source": [
-    "We observe that as we weakly impose the boundary condition, we no longer fullfill the equation to machine precision at the mesh vertices. We also plot the solution using `pyvista`"
+    "We observe that as we weakly impose the boundary condition,\n",
+    "we no longer fullfill the equation to machine precision at the mesh vertices.\n",
+    "We also plot the solution using {py:mod}`pyvista`"
    ]
   },
   {
@@ -206,7 +229,9 @@
    "metadata": {},
    "source": [
     "```{bibliography}\n",
-    "   :filter: cited and ({\"chapter1/nitsche\"} >= docnames)\n",
+    "   :filter: cited\n",
+    "   :labelprefix:\n",
+    "   :keyprefix: nitsche-\n",
     "```"
    ]
   }

chapter1/nitsche.py

@@ -16,8 +16,10 @@
 # # Weak imposition of Dirichlet conditions for the Poisson problem
 # Author: Jørgen S. Dokken
 #
-# In this section, we will go through how to solve the Poisson problem from the [Fundamentals](./fundamentals_code.ipynb) tutorial using Nitsche's method {cite}`Nitsche1971`.
-# The idea of weak imposition is that we add additional terms to the variational formulation to impose the boundary condition, instead of modifying the matrix system using strong imposition (lifting).
+# In this section, we will go through how to solve the Poisson problem from the
+# [Fundamentals](./fundamentals_code.ipynb) tutorial using Nitsche's method {cite}`nitsche-Nitsche1971`.
+# The idea of weak imposition is that we add additional terms to the variational formulation to
+# impose the boundary condition, instead of modifying the matrix system using strong imposition (lifting).
 #
 # We start by importing the required modules and creating the mesh and function space for our solution
 
@@ -44,7 +46,12 @@
 V = fem.functionspace(domain, ("Lagrange", 1))
 # -
 
-# Next, we create a function containing the exact solution (which will also be used in the Dirichlet boundary condition) and the corresponding source function for the right hand side. Note that we use `ufl.SpatialCoordinate` to define the exact solution, which in turn is interpolated into `uD` and used to create the source function `f`.
+# Next, we create a function containing the exact solution (which will also be used in the Dirichlet boundary condition)
+# and the corresponding source function for the right hand side. Note that we use {py:class}`ufl.SpatialCoordinate`
+# to define the exact solution, which in turn is interpolated into {py:class}`uD<dolfinx.fem.Function>`
+# by wrapping the {py:class}`ufl-expression<ufl.core.expr.Expr>` as a {py:class}`dolfinx.fem.Expression`.
+# For the source function, we use the symbolic differentiation capabilities of UFL to
+# compute the negative Laplacian of the exact solution, which we can use directly in the variational formulation.
 
 uD = fem.Function(V)
 x = SpatialCoordinate(domain)
@@ -54,21 +61,34 @@
 
 # As opposed to the first tutorial, we now have to have another look at the variational form.
 # We start by integrating the problem by parts, to obtain
+#
+# $$
 # \begin{align}
-#     \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d}x - \int_{\partial\Omega}\nabla u \cdot n v~\mathrm{d}s = \int_{\Omega} f v~\mathrm{d}x.
+#   \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d}x
+#   - \int_{\partial\Omega}\nabla u \cdot n v~\mathrm{d}s = \int_{\Omega} f v~\mathrm{d}x.
 # \end{align}
+# $$
+#
 # As we are not using strong enforcement, we do not set the trace of the test function to $0$ on the outer boundary.
 # Instead, we add the following two terms to the variational formulation
+#
+# $$
 # \begin{align}
-#     -\int_{\partial\Omega} \nabla  v \cdot n (u-u_D)~\mathrm{d}s + \frac{\alpha}{h} \int_{\partial\Omega} (u-u_D)v~\mathrm{d}s.
+#   -\int_{\partial\Omega} \nabla  v \cdot n (u-u_D)~\mathrm{d}s
+#   + \frac{\alpha}{h} \int_{\partial\Omega} (u-u_D)v~\mathrm{d}s.
 # \end{align}
+# $$
+#
 # where the first term enforces symmetry to the bilinear form, while the latter term enforces coercivity.
 # $u_D$ is the known Dirichlet condition, and $h$ is the diameter of the circumscribed sphere of the mesh element.
 # We create bilinear and linear form, $a$ and $L$
+#
+# $$
 # \begin{align}
 #     a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d}x + \int_{\partial\Omega}-(n \cdot\nabla u) v - (n \cdot \nabla v) u + \frac{\alpha}{h} uv~\mathrm{d}s,\\
 #     L(v) &= \int_{\Omega} fv~\mathrm{d}x + \int_{\partial\Omega} -(n \cdot \nabla v) u_D + \frac{\alpha}{h} u_Dv~\mathrm{d}s
 # \end{align}
+# $$
 
 u = TrialFunction(V)
 v = TestFunction(V)
@@ -80,7 +100,8 @@
 L = inner(f, v) * dx
 L += -inner(n, grad(v)) * uD * ds + alpha / h * inner(uD, v) * ds
 
-# As we now have the variational form, we can solve the linear problem
+# As we now have the variational form, we can solve the arising
+# {py:class}`linear problem<dolfinx.fem.petsc.LinearProblem>`
 
 problem = LinearProblem(a, L, petsc_options_prefix="nitsche_poisson")
 uh = problem.solve()
@@ -102,7 +123,9 @@
 if domain.comm.rank == 0:
     print(f"Error_max : {error_max:.2e}")
 
-# We observe that as we weakly impose the boundary condition, we no longer fullfill the equation to machine precision at the mesh vertices. We also plot the solution using `pyvista`
+# We observe that as we weakly impose the boundary condition,
+# we no longer fullfill the equation to machine precision at the mesh vertices.
+# We also plot the solution using {py:mod}`pyvista`
 
 # +
 import pyvista
@@ -120,5 +143,7 @@
 # -
 
 # ```{bibliography}
-#    :filter: cited and ({"chapter1/nitsche"} >= docnames)
+#    :filter: cited
+#    :labelprefix:
+#    :keyprefix: nitsche-
 # ```