Implementing Periodic BC for eigenvectors

dolfinx
1

Dear all,

I am trying to solve Helmholtz Equation for 3D geometries. My case includes degenerate modes (Same eigenvalue, shifted eigenvectors). So I am looking for a BC which enables me to calculate the eigenvector in a much more computationally efficient way. Here is my MWE for 2D disk;

import gmsh
gmsh.initialize()
r = 1
membrane = gmsh.model.occ.addDisk(0, 0, 0, r, r)
gmsh.model.occ.synchronize()

gdim = 2
status = gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
lc = 0.05
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",lc)
gmsh.model.mesh.generate(gdim)

from dolfinx.io import extract_gmsh_geometry, extract_gmsh_topology_and_markers, ufl_mesh_from_gmsh
from mpi4py import MPI
if MPI.COMM_WORLD.rank == 0:
    # Get mesh geometry
    geometry_data = extract_gmsh_geometry(gmsh.model)
    # Get mesh topology for each element
    topology_data = extract_gmsh_topology_and_markers(gmsh.model)

import numpy as np
if MPI.COMM_WORLD.rank == 0:
    # Extract the cell type and number of nodes per cell and broadcast
    # it to the other processors 
    gmsh_cell_type = list(topology_data.keys())[0]    
    properties = gmsh.model.mesh.getElementProperties(gmsh_cell_type)
    name, dim, order, num_nodes, local_coords, _ = properties
    cells = topology_data[gmsh_cell_type]["topology"]
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([gmsh_cell_type, num_nodes], root=0)
else:        
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([None, None], root=0)
    cells, geometry_data = np.empty([0, num_nodes]), np.empty([0, gdim])

from dolfinx.io import cell_perm_gmsh
from dolfinx.cpp.mesh import to_type
from dolfinx.mesh import create_mesh

# Permute topology data from MSH-ordering to dolfinx-ordering
ufl_domain = ufl_mesh_from_gmsh(cell_id, gdim)
gmsh_cell_perm = cell_perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]

# Create distributed mesh
mesh = create_mesh(MPI.COMM_WORLD, cells, geometry_data[:, :gdim], ufl_domain)

from ufl import TestFunction, TrialFunction, dx, grad, inner
from slepc4py import SLEPc
from dolfinx.fem import FunctionSpace, Function, assemble_matrix

V = FunctionSpace(mesh, ("CG", 1))

u = TrialFunction(V); v = TestFunction(V)

a = -inner(grad(u), grad(v))*dx
A = assemble_matrix(a); A.assemble()

c = inner(u , v) * dx
C = assemble_matrix(c); C.assemble()

target = np.pi**2
nev = 2
C = - C 

solver = SLEPc.EPS().create(MPI.COMM_WORLD) # A p = \lambda C p
solver.setOperators(A, C)
st = solver.getST()
st.setType('sinvert')
solver.setTarget(target)
solver.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)  
solver.setDimensions(nev, SLEPc.DECIDE)
solver.setTolerances(1e-15)
solver.setFromOptions()
solver.solve()

A = solver.getOperators()[0]
vr, vi = A.createVecs()

eig = solver.getEigenvalue(0)
omega = np.sqrt(eig)
solver.getEigenvector(0, vr, vi)
print(omega, np.sqrt(omega))

p = Function(V)
p.vector.setArray(vr.array)
p.x.scatter_forward()

import dolfinx.io
p.name = "Acoustic_Wave"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "p.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(p)

which gives me;

Screenshot from 2022-02-08 14-54-20

I want to implement a BC which should be able to result in same eigenvector with 1/8th of the full circle;

Screenshot from 2022-02-08 15-12-50

Is it possible?

Thanks for your interest in advance.

2

Consider: https://github.com/jorgensd/dolfinx_mpc
and in particular:
https://github.com/jorgensd/dolfinx_mpc/blob/master/python/demos/demo_periodic_gep.py

1 Like