spectral radius for scalar wave equation

[krober10nd]

• A good way I can estimate the numerically stable timestep and the "goodness" of a 2D/3D simplical mesh is by measuring the spectral radius of the scalar wave operator.
• dt < 2/sqrt(spectral_radius)
• Here I used KMV elements which admit diagonal mass matrices for spatial orders < 5.
• I've confirmed through ad hoc tests that the min. dihedral angle of the mesh reduces the spectral radius, thus increasing the maximum stable timestep.
• The spectral radius can be calculated like this (in Firedrake):

  
  import scipy
  import numpy as np

import firedrake as fd from firedrake import dot, grad import finat

def estimate_timestep(mesh, V, c, estimate_max_eigenvalue=True): """Estimate the maximum stable timestep based on the spectral radius using optionally the Gershgorin Circle Theorem to estimate the maximum generalized eigenvalue. Otherwise computes the maximum generalized eigenvalue exactly Currently only works with KMV elements. """

u, v = fd.TrialFunction(V), fd.TestFunction(V)
quad_rule = finat.quadrature.make_quadrature(
    V.finat_element.cell, V.ufl_element().degree(), "KMV"
)
dxlump = fd.dx(rule=quad_rule)
A = fd.assemble(u * v * dxlump)
ai, aj, av = A.petscmat.getValuesCSR()
av_inv = []
for value in av:
    if value == 0:
        av_inv.append(0.0)
    else:
        av_inv.append(1 / value)
Asp = scipy.sparse.csr_matrix((av, aj, ai))
Asp_inv = scipy.sparse.csr_matrix((av_inv, aj, ai))

K = fd.assemble(c*c*dot(grad(u), grad(v)) * dxlump)
ai, aj, av = K.petscmat.getValuesCSR()
Ksp = scipy.sparse.csr_matrix((av, aj, ai))

# operator
Lsp = Asp_inv.multiply(Ksp)
if estimate_max_eigenvalue:
    # absolute maximum of diagonals
    max_eigval = np.amax(np.abs(Lsp.diagonal()))
else:
    print(
        "Computing exact eigenvalues is extremely computationally demanding!",
        flush=True,
    )
    max_eigval = scipy.sparse.linalg.eigs(
        Ksp, M=Asp, k=1, which="LM", return_eigenvectors=False
    )[0]

# print(max_eigval)
max_dt = np.float(2 / np.sqrt(max_eigval))
print(
    f"Maximum stable timestep should be about: {np.float(2 / np.sqrt(max_eigval))} seconds",
    flush=True,
)
return max_dt

[krober10nd]

Re inputs: mesh is a Firedrake.Mesh, V is a Firedrake.FunctionSpace, and c is a Function(V) with spatially varying scalar wavespeeds.

[nschloe]

Oh man, I'd love to see a simple Python FEM package without C++ (beyond those in numpy/scipy). Anyway, this seems a little too specific with all the input argument. You'd want to pin one single value for each mesh. Looking at the spectral radius is actually a good idea, I'll add that for the Laplacian.

[krober10nd]

Oh man, I'd love to see a simple Python FEM package without C++ (beyond those in numpy/scipy)

Yea, I think scikit-fem strives for some of that. I can't tell you how much time I've spent (wasted?) trying to install Firedrake on different machines.

I agree it's a little too specific but the spectral radius in general is useful as it incorporates both the discretization (order and such) and the mesh quality at the same time.

For my application the difference between a minimum dihedral angle bound of 1 degree and 18 degrees means a variation in maximum stable timestep between 0.0001 and 0.0012 seconds. To get that high of a dihedral angle bound, I typically need to iterate with DistMesh 300-500 iterations (depending on the minimum element size).

[krober10nd]

Here's the result.