Merge upstream main into perf/tinyasm-dofmap-copy

.github/workflows/core.yml

@@ -536,7 +536,7 @@ jobs:
       - name: Setup flake8 annotations
         uses: rbialon/flake8-annotations@v1
       - name: Install linting packages
-        run: pip install flake8 pylint
+        run: pip install flake8 pylint cython-lint
       - name: Lint codebase
         run: make lint GITHUB_ACTIONS_FORMATTING=1
 

Makefile

@@ -34,6 +34,8 @@ srclint:
 	@python -m flake8 $(FLAKE8_FORMAT) pyop2/scripts --filename=*
 	@echo "    Linting TSFC"
 	@python -m flake8 $(FLAKE8_FORMAT) tsfc
+	@echo "    Linting Cython code"
+	@cython-lint firedrake/cython
 
 .PHONY: actionlint
 actionlint:

demos/demo_references.bib

@@ -515,3 +515,89 @@ @article{John:2017
 	journal = {{SIAM} Review}
 }
 
+@article{CarstensenGedicke:2014,
+  title={Guaranteed lower bounds for eigenvalues},
+  author={Carstensen, Carsten and Gedicke, Joscha},
+  journal={Mathematics of Computation},
+  volume={83},
+  number={290},
+  pages={2605--2629},
+  year={2014},
+  doi={10.1090/S0025-5718-2014-02833-0}
+}
+
+@article{Boffi:2010,
+    title={Finite element approximation of eigenvalue problems},
+    volume={19},
+    DOI={10.1017/S0962492910000012},
+    journal={Acta Numerica},
+    author={Boffi, Daniele},
+    year={2010},
+    pages={1--120}
+}
+
+@article{Boffi:2014,
+title = {Some remarks on finite element approximation of multiple eigenvalues},
+journal = {Applied Numerical Mathematics},
+volume = {79},
+pages = {18-28},
+year = {2014},
+doi = {10.1016/j.apnum.2012.08.006},
+author = {Daniele Boffi and Lucia Gastaldi},}
+}
+
+@article{Gander:2012,
+    author = {Martin J. Gander and Gerhard Wanner},
+    title = {From {Euler, Ritz, and Galerkin} to Modern Computing},
+    journal = {SIAM Review},
+    volume = {54},
+    number = {4},
+    pages = {627--666},
+    year = {2012},
+    doi = {10.1137/100804036},
+}
+
+@article{Crouzeix:1973,
+    doi = {10.1051/m2an/197307r300331},
+    year = 1973,
+    volume = {7},
+    number = {R3},
+    pages = {33--75},
+    author = {Michel Crouzeix and Pierre-Arnaud Raviart},
+    title = {Conforming and nonconforming finite element methods for solving the stationary {Stokes} equations {I}},
+    journal = {Revue fran{\c{c}}aise d'automatique informatique recherche op{\'{e}}rationnelle}
+}
+
+@article{Duran:2003,
+  title={A posteriori error estimates for the finite element approximation of eigenvalue problems},
+  author={Durán, Ricardo G. and Padra, Claudio and Rodríguez, Rodolfo},
+  journal={Mathematical Models and Methods in Applied Sciences},
+  volume={13},
+  number={8},
+  pages={1219--1229},
+  year={2003},
+  doi={10.1142/S0218202503002878}
+}
+
+@article{Larson:2000,
+  title={A posteriori and a priori error analysis for finite element approximations of self-adjoint elliptic eigenvalue problems},
+  author={Larson, Mats G.},
+  journal={SIAM Journal on Numerical Analysis},
+  volume={38},
+  number={2},
+  pages={608--625},
+  year={2000},
+  doi={10.1137/S0036142997320164}
+}
+
+@article{Hu:2014,
+  title={Lower bounds for eigenvalues of elliptic operators: By nonconforming finite element methods},
+  author={Jun Hu and Yunqing Huang and Qun Lin},
+  journal={Journal of Scientific Computing},
+  volume={61},
+  number={1},
+  pages={196--221},
+  year={2014},
+  doi={10.1007/s10915-014-9821-5}
+}
+

demos/poisson_adaptive_eigenvalue/poisson_adaptive_eigenvalue.py.rst

@@ -0,0 +1,218 @@
+Adaptive eigenvalue problem on an L-shaped domain
+=================================================
+
+.. rst-class:: emphasis
+
+    This demo computes the lowest eigenvalue of the Dirichlet Laplacian on an L-shaped domain to within a guaranteed
+    error tolerance using adaptive mesh refinement. The error estimate in the eigenvalue is based on solving the problem
+    with two complementary discretisations: a conforming method naturally yields an upper bound, while the
+    Crouzeix-Raviart non-conforming method yields a lower bound after postprocessing proposed by Carstensen & Gedicke.
+    The demo also demonstrates that adaptive mesh refinement is strictly more efficient than uniform refinement,
+    due to the singularity of the eigenfunction at the re-entrant corner.
+
+    The demo was contributed by `Patrick Farrell
+    <mailto:patrick.farrell@maths.ox.ac.uk>`__.
+
+We consider the Dirichlet eigenvalue problem for the Poisson equation on an L-shaped domain :math:`\Omega`:
+
+.. math::
+   -\Delta u &= \lambda u \quad \text{in } \Omega, \\
+           u &= 0         \quad \text{on } \partial \Omega.
+
+We aim to compute the lowest eigenvalue. Since the domain has a re-entrant corner, the eigenfunction has a singularity. Thus, uniform mesh refinement leads to suboptimal convergence. We demonstrate an adaptive strategy driven by a residual-based a posteriori error estimator. To establish rigorous bounds on the true eigenvalue :math:`\lambda`, we use conforming elements (CG) for an upper bound :cite:`Boffi:2010`, and nonconforming elements (CR) with Carstensen-Gedicke postprocessing :cite:`CarstensenGedicke:2014` for a guaranteed lower bound.
+
+We start by importing the necessary libraries: ::
+
+  from firedrake import *
+  from netgen.occ import *
+
+We then define the L-shaped domain using Netgen's Open CASCADE technology (OCC) interface, and generate an initial mesh: ::
+
+  rect1 = WorkPlane(Axes((0,0,0), n=Z, h=X)).Rectangle(1,2).Face()
+  rect2 = WorkPlane(Axes((0,1,0), n=Z, h=X)).Rectangle(2,1).Face()
+  L = rect1 + rect2
+
+  geo = OCCGeometry(L, dim=2)
+  ngmesh = geo.GenerateMesh(maxh=0.5)
+  mesh = Mesh(ngmesh)
+
+We create a function to solve the eigenvalue problem for both continuous (CG) and nonconforming Crouzeix-Raviart (CR) :cite:`Crouzeix:1973` elements. By the Rayleigh-Ritz principle :cite:`Boffi:2010,Gander:2012`, the CG solution provides an upper bound :math:`\lambda_{\text{ub}}`. The CR solution gives a discrete eigenvalue :math:`\lambda_{\text{CR}}`, which can be postprocessed to yield a guaranteed lower bound:
+
+.. math::
+    \lambda_{\text{lb}} = \frac{\lambda_{\text{CR}}}{1 + \kappa_{\text{CR}}^2 h_{\max}^2 \lambda_{\text{CR}}}
+
+where :math:`\kappa_{\text{CR}} \approx 0.1893` is the constant established by Carstensen and Gedicke :cite:`CarstensenGedicke:2014`. A general theory for deriving lower bounds for eigenvalues with nonconforming methods has been developed by Hu et al. :cite:`Hu:2014`. We will use the postprocessed lower bound to terminate the adaptive iteration, while for technical reasons we will plot the Galerkin gap :math:`\lambda_{\text{CG}} - \lambda_{\text{CR}}` to demonstrate optimal convergence.
+
+To efficiently compute the smallest eigenvalue, we configure SLEPc using a solver parameters dictionary. We specify a Krylov-Schur eigensolver (``eps_type``) and a shift-and-invert spectral transformation (``st_type``) with a target of zero (``eps_target``). We also flag the generalized eigenvalue problem as Hermitian (``eps_gen_hermitian``) and request the smallest real eigenvalue (``eps_smallest_real``). ::
+
+  def solve_eigenproblem(mesh):
+      h_symbolic = CellDiameter(mesh)
+      DG0 = FunctionSpace(mesh, "DG", 0)
+      h_cell = Function(DG0).interpolate(h_symbolic)
+      with h_cell.dat.vec_ro as hvec:
+          h_max = hvec.max()[1]
+      eigenfunction = None
+
+      for space in ["CG", "CR"]:
+          V = FunctionSpace(mesh, space, 1)
+
+          u = TrialFunction(V)
+          v = TestFunction(V)
+
+          a = inner(grad(u), grad(v))*dx
+          b = inner(u, v)*dx
+          bc = DirichletBC(V, 0, "on_boundary")
+          eigenproblem = LinearEigenproblem(a, b, bc)
+
+          sp = {
+                "eps_gen_hermitian": None,
+                "eps_smallest_real": None,
+                "eps_monitor_cancel": None,
+                "eps_type": "krylovschur",
+                "eps_target": 0,
+                "st_type": "sinvert",
+                }
+
+          eigensolver = LinearEigensolver(eigenproblem, 1, solver_parameters=sp)
+          eigensolver.solve()
+
+          if space == "CR":
+              lambda_CR = eigensolver.eigenvalue(0).real
+              kappa_CR = 0.1893
+              lambda_lb = lambda_CR / (1 + kappa_CR**2 * h_max**2 * lambda_CR)
+          if space == "CG":
+              lambda_ub = eigensolver.eigenvalue(0).real
+              eigenfunction = eigensolver.eigenfunction(0)[0]
+
+      eigenfunction.rename("Eigenfunction")
+      return (lambda_lb, lambda_ub, lambda_CR, eigenfunction)
+
+These bounds do not describe where the mesh should be refined so as to reduce the error. For this purpose we employ a standard residual-based a posteriori error estimator :cite:`Duran:2003,Larson:2000`. Note that this assumes there is a single eigenfunction associated with the lowest eigenvalue; if the eigenvalue were of higher multiplicity the estimator would need to consider the entire eigenspace :cite:`Boffi:2014`. ::
+
+  def estimate_error(mesh, uh, lam):
+      W = FunctionSpace(mesh, "DG", 0)
+      eta_sq = Function(W)
+      w = TestFunction(W)
+      h = CellDiameter(mesh)
+      n = FacetNormal(mesh)
+      v = CellVolume(mesh)
+
+      G = (
+            inner(eta_sq / v, w)*dx
+          - inner(h**2 * (lam*uh + div(grad(uh)))**2, w) * dx
+          - inner(h('+')/2 * jump(grad(uh), n)**2, w('+')) * dS
+          - inner(h('-')/2 * jump(grad(uh), n)**2, w('-')) * dS
+          )
+
+      sp = {"mat_type": "matfree", "ksp_type": "preonly", "pc_type": "jacobi"}
+      solve(G == 0, eta_sq, solver_parameters=sp)
+      eta = Function(W).interpolate(sqrt(eta_sq))
+
+      with eta.dat.vec_ro as eta_:
+          error_est = eta_.norm()
+      return (eta, error_est)
+
+We define a function to adapt the mesh by refining elements with large error indicators, using the maximum Dörfler-like marking strategy with :math:`\theta = 0.5`: ::
+
+  def adapt(mesh, eta):
+      W = FunctionSpace(mesh, "DG", 0)
+      markers = Function(W)
+
+      with eta.dat.vec_ro as eta_:
+          eta_max = eta_.max()[1]
+
+      theta = 0.5
+      should_refine = conditional(gt(eta, theta*eta_max), 1, 0)
+      markers.interpolate(should_refine)
+
+      return mesh.refine_marked_elements(markers)
+
+Finally, we run the adaptive loop until the upper and lower bounds agree to within a tolerance. ::
+
+  max_iterations = 20
+  error_estimators = []
+  dofs = []
+  err = 1
+
+  for i in range(max_iterations):
+      lam_lb, lam_ub, lam_CR, uh = solve_eigenproblem(mesh)
+      err = lam_ub - lam_lb
+      gap = lam_ub - lam_CR
+      error_estimators.append(gap)
+      dofs.append(uh.function_space().dim())
+      print(f"Level {i}: Upper bound {lam_ub:.5f}, Lower bound {lam_lb:.5f}, Bound gap {err:.5e}, Galerkin gap {gap:.5e}")
+
+      VTKFile(f"l_eigenfunction_{i}.pvd").write(uh)
+
+      if err < 1e-2 or dofs[-1] > 1000000:
+          break
+
+      eta, _ = estimate_error(mesh, uh, lam_ub)
+      mesh = adapt(mesh, eta)
+
+To demonstrate that adaptivity is necessary to achieve the optimal convergence rate, we can run the same script with :math:`\theta = 0`, which forces uniform refinement (all cells are marked for refinement at every step). We make this optional by guarding it behind the Boolean ``run_uniform``. ::
+
+  run_uniform = True
+
+  if run_uniform:
+      mesh_uniform = Mesh(ngmesh)
+      uniform_error_estimators = []
+      uniform_dofs = []
+
+      def adapt_uniform(mesh, eta):
+          markers = Function(FunctionSpace(mesh, "DG", 0)).assign(1.0)
+          return mesh.refine_marked_elements(markers)
+
+      for i in range(max_iterations):
+          lam_lb, lam_ub, lam_CR, uh = solve_eigenproblem(mesh_uniform)
+          err = lam_ub - lam_lb
+          gap = lam_ub - lam_CR
+          uniform_error_estimators.append(gap)
+          uniform_dofs.append(uh.function_space().dim())
+          if err < 5e-3 or uniform_dofs[-1] > 1000000:
+              break
+          eta, _ = estimate_error(mesh_uniform, uh, lam_ub)
+          mesh_uniform = adapt_uniform(mesh_uniform, eta)
+
+We can plot the convergence of the Galerkin gap :math:`\lambda_{\text{ub}} - \lambda_{\text{CR}}` against the number of degrees of freedom. With adaptivity, we achieve the optimal :math:`O(N^{-1})` convergence rate. For uniform refinement, the error is initially dominated by the smooth part of the solution (yielding a pre-asymptotic :math:`O(N^{-1})` rate), but as the mesh is refined, the singularity inevitably dominates and limits the asymptotic convergence to the suboptimal rate of :math:`O(N^{-2/3})`. ::
+
+  try:
+      import matplotlib.pyplot as plt
+      import numpy as np
+
+      plt.grid()
+      plt.loglog(dofs, error_estimators, '-ok', label=r"Adaptive refinement ($\theta = 0.5$)")
+      scaling = error_estimators[0] / dofs[0]**-1
+      plt.loglog(dofs, np.array(dofs)**(-1.0) * scaling, '--', label="Optimal convergence $N^{-1}$")
+
+      if run_uniform:
+          plt.loglog(uniform_dofs, uniform_error_estimators, '-or', label=r"Uniform refinement ($\theta = 0$)")
+          scaling_uniform = uniform_error_estimators[-1] / uniform_dofs[-1]**(-2.0/3.0)
+          plt.loglog(uniform_dofs, np.array(uniform_dofs)**(-2.0/3.0) * scaling_uniform, ':r', label=r"Suboptimal convergence $N^{-2/3}$")
+      plt.xlabel("Number of degrees of freedom $N$")
+      plt.ylabel(r"Galerkin gap $\lambda_{\text{ub}} - \lambda_{\text{CR}}$")
+      plt.legend()
+      plt.savefig("adaptive_eigenvalue_convergence.png")
+  except ImportError:
+      warning("Matplotlib not imported")
+
+.. figure:: adaptive_eigenvalue_convergence.png
+   :align: center
+   :figwidth: 80%
+
+   Convergence of the Galerkin gap :math:`\lambda_{\text{ub}} - \lambda_{\text{CR}}`. Note that the adaptive scheme achieves the optimal convergence rate of :math:`O(N^{-1})`, whereas uniform refinement is limited to the suboptimal rate of :math:`O(N^{-2/3})`.
+
+To visualize how the adaptive algorithm resolves the singularity, the sequence of images below shows the mesh (at refinement level 15), zooming into the re-entrant corner at 10x, 100x, and 1000x magnification.
+
+.. image:: mesh_zoom_10x.png
+   :width: 32%
+.. image:: mesh_zoom_100x.png
+   :width: 32%
+.. image:: mesh_zoom_1000x.png
+   :width: 32%
+
+.. rubric:: References
+
+.. bibliography:: demo_references.bib
+   :filter: docname in docnames
+

docs/source/advanced_tut.rst

@@ -38,3 +38,4 @@ element systems.
    Coupled volume-surface reaction-diffusion on a torus using Submesh and geometric multigrid.<demos/submesh_reaction_diffusion.py>
    Nonlinear preconditioning using an auxiliary SNES for the Allen-Cahn equation.<demos/nonlinear_pc_allen_cahn.py>
    Reynolds-robust preconditioning of the stationary Navier-Stokes equations.<demos/reynolds_robust_navier_stokes_hdiv.py>
+   Eigenvalue problem with guaranteed bounds and adaptive refinement.<demos/poisson_adaptive_eigenvalue.py>

docs/source/conf.py

@@ -193,6 +193,7 @@
 
 # Dodgy links
 linkcheck_ignore = [
+    # journals and institutions
     r'https://zenodo.org/.*',
     r'https://doi\.org/.*',
     r'https://epubs\.siam\.org/doi/.*',
@@ -201,17 +202,23 @@
     r'https://ieeexplore\.ieee\.org/document/1634311/',
     r'http://www.cs.virginia.edu/stream/',
     r'https://www.sciencedirect.com',
-    r'https://.*\.baylor\.edu.*',
-    r'https://www.crosscountrytrains.co.uk/',
     r'https://www.siam.org/',
+    r'https://.*\.baylor\.edu.*',
     r'https://aims.ac.rw',
     r'https://mpecdt.ac.uk',
+    r'https://www.ox.ac.uk',
+    r'https://le.ac.uk/people/',
+
+    # event links
+    r'https://www.crosscountrytrains.co.uk/',
     r'https://www.hilton.com/en/hotels/leehnhn-hilton-leeds-city/',
     r'https://www.radissonhotels.com/*',
     r'https://all.accor.com/hotel/*',
     r'https://fluids.leeds.ac.uk/',
-    r'https://www.ox.ac.uk',
     r'https://buy.crosscountrytrains.co.uk',
+
+    # other
+    r'https://firedrakeproject.slack.com',
     r'https://join.slack.com/t/firedrakeproject/*',
 ]
 linkcheck_timeout = 30

firedrake/bcs.py

@@ -176,8 +176,8 @@ def hermite_stride(bcnodes):
         with temp_internal_comm(self._function_space.mesh().comm) as icomm:
             num_global_nodes = icomm.reduce(len(bcnodes), MPI.SUM, root=0)
             if num_global_nodes == 0 and icomm.rank == 0:
-                logger.warn(f"Subdomain {self.sub_domain} is empty. This is likely an error. "
-                            "Did you choose the right label?")
+                logger.warning(f"Subdomain {self.sub_domain} is empty. This is likely an error. "
+                               "Did you choose the right label?")
 
         return bcnodes
 

firedrake/cofunction.py

@@ -519,6 +519,7 @@ def __call__(self, value):
             else:
                 solve, rhs, soln = self._solver
                 rhs.assign(value)
+                soln.zero()
                 solve()
                 output = Function(self._function_space)
                 output.assign(soln)