Linear stability analysis functionality (#176)

* Helmholtz eigenvalue problem example

* [WIP] testing direct stability on cylinder

* Link to old code with working SLEPc call

* Derive linear control term

* Simple test problems for KS algorithm

* Arnoldi example in Firedrake with heat equation

* Debugging KS in Firedrake

* Working Krylov-Schur on heat equation

* Add linearized NS solver

* Add base flow forcing terms

* Fixed Arnoldi iteration with Navier-Stokes mass matrix

* Tweaking Arnoldi config

* Testing shift-invert on heat eqution

* Working cylinder stability with spectral transformation

* Add spectral transformation script

* Sorting ritz vals

* Partial cleanup of cyl_st script

* Cleaned up cyl_st.  Save before making the analysis generic

* Output formatting

* Non-crashing stability analysis for cavity

* Working shift-invert

* Add general-purpose 'residual' method

* Filtering duplicate eigenvalues

* Working general functions with filtered eigenvalues

* Shift-invert Arnoldi test on cavity

* [WIP] adding CLI args to cyl script

* Clean up examples/cylinder/stability/

* Combined Arnoldi and K-S into one function

* Move  operatorto __matmul__

* Pre-assemble linear operator

* Separate code for  linear stability and linear operators

* Cavity stability analysis

* Fix render method for cavity and step

* Move common functionality to utils.stability

* Black

* yapf

* Fix codespell

* Undo re-running gmsh:

* Remove examples/cavity/stability/

* Final cleanup of examples/

* codespell complaint

examples/cavity/stability.py

@@ -0,0 +1,179 @@
+"""
+To get the leading mode, run with:
+
+```
+mpiexec -np 20 python stability.py --output-dir eig_output --sigma 1.0+11j
+```
+"""
+
+import argparse
+import os
+
+# from cyl_common import base_checkpoint, evec_checkpoint, flow
+from typing import NamedTuple
+
+import firedrake as fd
+import numpy as np
+
+import hydrogym.firedrake as hgym
+
+parser = argparse.ArgumentParser(
+    description="Stability analysis of the open cavity flow.")
+parser.add_argument(
+    "--mesh",
+    default="fine",
+    type=str,
+    help="Identifier for the mesh resolution",
+)
+parser.add_argument(
+    "--reynolds",
+    default=7500.0,
+    type=float,
+    help="Reynolds number of the flow",
+)
+parser.add_argument(
+    "--krylov-dim",
+    default=100,
+    type=int,
+    help="Dimension of the Krylov subspace (number of Arnoldi vectors)",
+)
+parser.add_argument(
+    "--tol",
+    default=1e-10,
+    type=float,
+    help="Tolerance to use for determining converged eigenvalues.",
+)
+parser.add_argument(
+    "--schur",
+    action="store_true",
+    dest="schur",
+    default=False,
+    help="Use Krylov-Schur iteration to restart the Arnoldi process.",
+)
+parser.add_argument(
+    "--no-adjoint",
+    action="store_true",
+    default=False,
+    help="Skip computing the adjoint modes.",
+)
+parser.add_argument(
+    "--sigma",
+    default=0.0,
+    type=complex,
+    help="Shift for the shift-invert Arnoldi method.",
+)
+parser.add_argument(
+    "--base-flow",
+    type=str,
+    help="Path to the HDF5 checkpoint containing the base flow.",
+)
+parser.add_argument(
+    "--output-dir",
+    type=str,
+    default="eig_output",
+    help="Directory in which output files will be stored.",
+)
+if __name__ == "__main__":
+  args = parser.parse_args()
+
+  mesh = args.mesh
+  m = args.krylov_dim
+  sigma = args.sigma
+  tol = args.tol
+  Re = args.reynolds
+
+  output_dir = args.output_dir
+  if not os.path.exists(output_dir):
+    os.makedirs(output_dir)
+
+  velocity_order = 2
+  stabilization = "none"
+
+  flow = hgym.Cavity(
+      Re=Re,
+      mesh=mesh,
+      velocity_order=velocity_order,
+  )
+
+  dof = flow.mixed_space.dim()
+
+  hgym.print("|--------------------------------------------------|")
+  hgym.print("| Linear stability analysis of the open cavity flow |")
+  hgym.print("|--------------------------------------------------|")
+  hgym.print(f"Reynolds number:       {Re}")
+  hgym.print(f"Krylov dimension:      {m}")
+  hgym.print(f"Spectral shift:        {sigma}")
+  hgym.print(f"Include adjoint modes: {not args.no_adjoint}")
+  hgym.print(f"Krylov-Schur restart:  {args.schur}")
+  hgym.print(f"Total dof: {dof} --- dof/rank: {int(dof/fd.COMM_WORLD.size)}")
+  hgym.print("")
+
+  if args.base_flow:
+    hgym.print(f"Loading base flow from checkpoint {args.base_flow}...")
+    flow.load_checkpoint(args.base_flow)
+
+  else:
+
+    # Since this flow is at high Reynolds number we have to
+    #    ramp to get the steady state
+    hgym.print("Solving the steady-state problem...")
+
+    steady_solver = hgym.NewtonSolver(
+        flow,
+        stabilization=stabilization,
+        solver_parameters={"snes_monitor": None})
+    Re_init = [500, 1000, 2000, 4000, Re]
+
+    for i, Re in enumerate(Re_init):
+      flow.Re.assign(Re)
+      hgym.print(f"Steady solve at Re={Re_init[i]}")
+      qB = steady_solver.solve()
+
+    flow.save_checkpoint(f"{args.output_dir}/base.h5")
+
+  hgym.print("Computing direct modes...")
+  dir_results = hgym.utils.stability_analysis(
+      flow, sigma, m, tol, schur_restart=args.schur, adjoint=False)
+  np.save(f"{output_dir}/raw_evals", dir_results.raw_evals)
+
+  # Save checkpoints
+  evals = dir_results.evals
+  eval_data = np.zeros((len(evals), 3), dtype=np.float64)
+  eval_data[:, 0] = evals.real
+  eval_data[:, 1] = evals.imag
+  eval_data[:, 2] = dir_results.residuals
+  np.save(f"{output_dir}/evals", eval_data)
+
+  with fd.CheckpointFile(f"{output_dir}/evecs.h5", "w") as chk:
+    for i in range(len(evals)):
+      chk.save_mesh(flow.mesh)
+      dir_results.evecs_real[i].rename(f"evec_{i}_re")
+      chk.save_function(dir_results.evecs_real[i])
+      dir_results.evecs_imag[i].rename(f"evec_{i}_im")
+      chk.save_function(dir_results.evecs_imag[i])
+
+  if not args.no_adjoint:
+    hgym.print("Computing adjoint modes...")
+    adj_results = hgym.utils.stability_analysis(
+        flow, sigma, m, tol, adjoint=True)
+
+    # Save checkpoints
+    evals = adj_results.evals
+    eval_data = np.zeros((len(evals), 3), dtype=np.float64)
+    eval_data[:, 0] = evals.real
+    eval_data[:, 1] = evals.imag
+    eval_data[:, 2] = adj_results.residuals
+    np.save(f"{output_dir}/adj_evals", eval_data)
+
+    with fd.CheckpointFile(f"{output_dir}/adj_evecs.h5", "w") as chk:
+      for i in range(len(evals)):
+        chk.save_mesh(flow.mesh)
+        adj_results.evecs_real[i].rename(f"evec_{i}_re")
+        chk.save_function(adj_results.evecs_real[i])
+        adj_results.evecs_imag[i].rename(f"evec_{i}_im")
+        chk.save_function(adj_results.evecs_imag[i])
+
+  hgym.print(
+      "NOTE: If there is a warning following this, ignore it.  It is raised by PETSc "
+      "CLI argument handling and not argparse. It does not indicate that any CLI "
+      "arguments are ignored by this script.")

examples/cylinder/lti_system.py

@@ -0,0 +1,65 @@
+"""Illustrate constructing an LTI system from a flow configuration.
+
+This is a work in progress - so far it's just the base flow and the control
+vector. Ultimately it will need LinearOperator functionality to have the form
+
+````
+x' = A*x + B*u
+y  = C*x + D*u
+```
+
+where `x` is a firedrake.Function and `u` and `y` are numpy arrays.
+"""
+import os
+
+import firedrake as fd
+import matplotlib.pyplot as plt
+from firedrake.pyplot import tripcolor
+from ufl import dx, inner
+
+import hydrogym.firedrake as hgym
+
+show_plots = False
+output_dir = "output"
+if not os.path.exists(output_dir):
+  os.makedirs(output_dir)
+
+flow = hgym.RotaryCylinder(Re=100, mesh="medium")
+
+# 1. Compute base flow
+steady_solver = hgym.NewtonSolver(flow)
+qB = steady_solver.solve()
+qB.rename("qB")
+
+# Check lift/drag
+print(flow.compute_forces(qB))
+
+if show_plots:
+  vort = flow.vorticity(qB.subfunctions[0])
+  fig, ax = plt.subplots(1, 1, figsize=(6, 3))
+  tripcolor(vort, axes=ax, cmap="RdBu", vmin=-2, vmax=2)
+
+# 2. Derive flow field associated with actuation BC
+# See Barbagallo et al. (2009) for details on the "lifting" procedure
+F = steady_solver.steady_form(qB)  # Nonlinear variational form
+J = fd.derivative(F, qB)  # Jacobian with automatic differentiation
+
+flow.linearize_bcs(mixed=True)
+flow.set_control([1.0])
+bcs = flow.collect_bcs()
+
+# Solve steady, inhomogeneous problem
+qC = fd.Function(flow.mixed_space, name="qC")
+v, s = fd.TestFunctions(flow.mixed_space)
+zero = inner(fd.Constant((0.0, 0.0)), v) * dx
+fd.solve(J == zero, qC, bcs=bcs)
+
+if show_plots:
+  vort = flow.vorticity(qC.subfunctions[0])
+  fig, ax = plt.subplots(1, 1, figsize=(6, 3))
+  tripcolor(vort, axes=ax, cmap="RdBu", vmin=-2, vmax=2)
+
+with fd.CheckpointFile(f"{output_dir}/lin_fields.h5", "w") as chk:
+  chk.save_mesh(flow.mesh)
+  chk.save_function(qB)
+  chk.save_function(qC)

examples/cylinder/pressure-probes.py

@@ -1,6 +1,7 @@
 """
 Simulate the flow around the cylinder with evenly spaced surface pressure probes.
 """
+
 import os
 import matplotlib.pyplot as plt
 import numpy as np