Helmholtz equation on a curved manifold

[Mathieu]

Hi everyone,

I work with legacy FEniCS.
I’m trying to solve the Helmholtz equation (generalized eigenproblem Ax = lambda Bx) on a 2D surface which is not a plane (like a half-sphere or the lateral surface of a cylinder). An example of such a mesh is available here: Dropbox (.h5 file). The associated .xdmf file is here: Dropbox.

What i’am doing for now is the following code:

from dolfin import *

# Read the mesh
mesh = Mesh()
    with XDMFFile(mesh.mpi_comm(), mesh_walls.xdmf) as file:
        file.read(mesh)
fs = FunctionSpace(mesh, "CG", 1)
gdim = mesh.topology().dim()
global_num_vertices = mesh.num_entities_global(0)
global_num_cells = mesh.num_entities_global(gdim)
print(f"      Number of vertices: {global_num_vertices}") #  27866
print(f"      Number of cells:  {global_num_cells}") #  55428

# Variational form
u = TrialFunction(fs)
v = TestFunction(fs)

a = inner(grad(u), grad(v)) * dx
b = inner(u, v) * dx

# Assembling
a_petsc = PETScMatrix(fs.mesh().mpi_comm())
assemble(a, tensor=a_petsc)
b_petsc = PETScMatrix(fs.mesh().mpi_comm())
assemble(b, tensor=b_petsc)

# Solving
eps = SLEPc.EPS().create(mesh.mpi_comm())
eps.setOperators(a_petsc.mat(), b_petsc.mat())
eps.setProblemType(eps.ProblemType.GHEP)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
eps.setTolerances(1e-6, 100000)
eps.setConvergenceTest(SLEPc.EPS.Conv.REL)
eps.setWhichEigenpairs(eps.Which.SMALLEST_MAGNITUDE)
eps.setDimensions(20, 40, SLEPc.DECIDE)
eps.view()
eps.solve()
eps.errorView()
n_modes_converged = eps.getConverged()

# Create the results vectors
vr, vi = a_petsc.mat().getVecs()
for i in range(n_modes_converged):
    k = eps.getEigenpair(i, vr, vi)
    error = eps.computeError(i)
    if k.real > 0:
        print(k.real)

This code works well for 3D meshes or 2D planes meshes. However if I use it for 2D curved manifolds (like the mesh of the example), I don’t get any error but the results are false because the curvature of the surface is not taken into account. How can I account for the curvature?

[dokken]

Have you considered the paper:

gmd-6-2099-2013.pdf

Or in DOLFINx:

Curvature by finite elements

I could solve the problem. It seems that project does not work well as indicated here so I used interpolate instead. In addition I defined the identity as an expression-although I am not sure if that was also affecting the code-as: Identity = Expression(("x[0]","x[1]"), degree=1) Idt = interpolate(Identity, V) Finally, here Rogens indicates that the approximation only works with linear elements, so I changed the degree of the vector space to 1 and it worked. I am going to run some extra simula…

[Mathieu]

Yes I read the paper of Rognes, but it seems like they only deal with linear problem of the form Ax=b. It is a bit different when you deal with Eigenproblem, particularly with the syntax. But do we agree on that: The only thing I need to change is the following variational form:

a = inner(grad(u), grad(v)) * dx
b = inner(u, v) * dx

that must be adapted for manifold, isn’t it?

[dokken]

I don’t think you should need to adapt it.
I’ve made a version that runs in DOLFINx (which also supports second order geometries)
Generate the geometry with:

// Gmsh project created on Tue Jul 30 10:39:54 2024
SetFactory("OpenCASCADE");
Mesh.CharacteristicLengthMin=0.05;
Mesh.CharacteristicLengthMax=0.05;
Mesh.RecombinationAlgorithm = 0;
Mesh.RecombineAll= 0;
Mesh.ElementOrder = 1;
Circle(1) = {0, 0, 0, 0.5, 0, 2*Pi};

Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 1} {Surface{1};Layers {100};}

Physical Surface("cylinder", 5) = {2};

and the finite element code with

from mpi4py import MPI
from slepc4py import SLEPc
import dolfinx.fem.petsc
import ufl
mesh, _, _ = dolfinx.io.gmshio.read_from_msh(
    "mesh.msh", MPI.COMM_WORLD, 0, gdim=3)
fs = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
gdim = mesh.geometry.dim
global_num_vertices = mesh.topology.index_map(0).size_global
global_num_cells = mesh.topology.index_map(mesh.topology.dim).size_global
print(f"      Number of vertices: {global_num_vertices}")  # 27866
print(f"      Number of cells:  {global_num_cells}")  # 55428
print(mesh.geometry.cmap.degree)

# Variational form
u = ufl.TrialFunction(fs)
v = ufl.TestFunction(fs)

a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
b = ufl.inner(u, v) * ufl.dx

A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a))
b_petsc = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(b))
A.assemble()
b_petsc.assemble()
# Solving
eps = SLEPc.EPS().create(mesh.comm)
eps.setOperators(A, b_petsc)
eps.setProblemType(eps.ProblemType.GHEP)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
eps.setTolerances(1e-6, 100000)
eps.setConvergenceTest(SLEPc.EPS.Conv.REL)
eps.setWhichEigenpairs(eps.Which.SMALLEST_MAGNITUDE)
eps.setDimensions(20, 40, SLEPc.DECIDE)
eps.view()
eps.solve()
eps.errorView()
n_modes_converged = eps.getConverged()

# Create the results vectors
vr, vi = A.getVecs()
for i in range(n_modes_converged):
    k = eps.getEigenpair(i, vr, vi)
    error = eps.computeError(i)
    if k.real > 0:
        print(k.real)

where ElementOrder=1 yields

4.003317762773205
4.003317762798583
9.869604414452821
13.872925058217891
13.872925058262837
16.013303961584196
16.013303961637565
25.882947090313447
25.882947090321526
36.03025378613448
36.030253786146886
39.478418458417735
43.48175501280492
43.48175501290896
45.90004669141309
45.90004669148018
55.49190611901173
55.49190611911046
64.05530840216181

and element order 2 yields

4.000000824105506
4.0000008241849585
9.869604414447382
13.869608070341062
13.869608070357602
16.000023042694973
16.000023042708833
25.869665994250983
25.869665994309962
36.000322402222615
36.000322402245196
39.4784184584009
43.47843787667883
43.478437876924744
45.870114983222244
45.87011498328191
55.47862446266824
55.47862446275751
64.0019744922168

What are the eigenvalues you would expect?

[Mathieu]

Thanks Dokken for your quick reply. What I expect is the following: when your geometry has a symmetry for instance a circular symmetry (as it is the case in your exemple if I understand well), all the eigenvalues related to azimuthal modes (i.e: depends on theta) must be degenerated (double eigenvalue). But, if your break the symmetry, (for instance the lateral surface of an oval cylinder), the eigenvalues shouldn’t be degenerated. A simple test case can be just to compute the eigenvalues on 1D circle and on 1D ellipse. Eigenvalues on the 1D circle must be degenerated whereas ones on the 1D ellipse shouldn’t be.

[Mathieu]

Hello Dokken, what happens if you run the same calculation on an ellipse instead of a circle?

[dokken]

Running it on:

SetFactory("OpenCASCADE");
Mesh.CharacteristicLengthMin=0.05;
Mesh.CharacteristicLengthMax=0.05;
Mesh.RecombinationAlgorithm = 0;
Mesh.RecombineAll= 0;
Mesh.ElementOrder = 1;
Ellipse(1) = {0, 0, 0, 1, 0.5, 0, 2*Pi};

Curve Loop(1) = {1};
Plane Surface(1) = {1};
Extrude {0, 0, 1} {Surface{1};Layers {100};}

Physical Surface("cylinder", 5) = {2};

yields

6.418284128717607e-12
1.6832879501991458
1.683287950227771
6.733154266059622
6.733154266681479
9.869604414428064
11.552893025617017
11.552893025957019
15.149621092820258
15.149621118787046
16.60276616991171
16.60276617386619
25.019260505949525
25.019260506070154
26.93277437127246
26.93277446713447
36.80248598132148
36.802486006530714
39.47841845837415
41.1617121244406

for first order geometry and

1.6823310834368876
1.6823315666473102
6.729326136740639
6.7293261896099
9.8696044143546
11.551936153872989
11.551936617939162
15.141000027678729
15.141003706105918
16.59893800710186
16.5989380178191
25.010639265248866
25.010643034341008
26.917431026459415
26.91744174824543
36.78714236767646
36.78715318440688
39.478418458190006
41.160755561183024

with second order geometry.

[Mathieu]

Thanks a lot Dokken for your interest in my issue !
I see that you have degenerated frequencies here on the Ellipse (twice the same eigenfreq). This is not correct as the symmetry is broken when using an Ellipse instead of a Circle, the eigenfrequencies must not be degenerated. I think that this issue is related to the fact that the Ellipse has a non constant curvature K(theta) and the solver don’t handle this. The solver see the Circle and the Ellipse as the same 1-dimensional line, and the transformation x → R * theta for the circle is ok (I mean the gradient operator is the same d/dx = 1/R*d/dtheta), but it is false for the Ellipse.

[dokken]

I think we need to split between solver and assembler.

From an assembly point of view:
assemble_matrix does account for the “curvature” of the cell, in the sense that the Jacobian is non-constant on the manifold in the case of a second order mesh.

This has for instance been shown in:

Normal vector in manifold

Note that when you use linear elements (straight elements), the divergence of the normal is 0 on every element (since the normal is constant on every line segment. If you run your code for a while, you will see that the normal is updated, but H will always remain 0, as the elements are not curved. If you use a second order mesh, H will be non-zero. import os from mpi4py import MPI from petsc4py import PETSc import numpy as np import math import ufl from basix.ufl import element, mixed_elem…

If SLEPc handles it or not, I do not know (or if there is something missing from the variational form). This is outside my area of expertise.