Non periodic corner - eigenvalue problem with Slepc

dolfinx
,
1 
import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, plot
import ufl
import pyvista as pv
from petsc4py import PETSc
from slepc4py import SLEPc
from dolfinx_mpc import MultiPointConstraint, assemble_matrix

def slave_boundary_x(x):
    return np.isclose(x[0], 1.0)

def slave_boundary_y(x):
    return np.logical_and(np.isclose(x[1], 1.0), np.logical_not(np.isclose(x[0], 1.0)))
    # return np.isclose(x[1], 1.0) # first try

# Second try
# def slave_corner(x):
#     return np.logical_and(np.isclose(x[0],1.0),
#                           np.isclose(x[1],1.0))

# mapping from slave -> master
def periodic_relation_x(x):
    out = np.zeros_like(x)
    out[0] = x[0] - 1.0
    out[1] = x[1]
    return out

# mapping from slave -> master
def periodic_relation_y(x):
    out = np.zeros_like(x)
    out[0] = x[0]
    out[1] = x[1] - 1.0
    return out

# Second try
# def periodic_relation_corner(x):
#     out = np.zeros_like(x)
#     out[0] = x[0] - 1.0
#     out[1] = x[1] - 1.0
#     return out

# Mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

# Function space
V = fem.functionspace(domain, ("CG", 1))
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, slave_boundary_x, periodic_relation_x, [])
mpc.create_periodic_constraint_geometrical(V, slave_boundary_y, periodic_relation_y, [])
# mpc.create_periodic_constraint_geometrical(V, slave_corner, periodic_relation_corner, []) # second try
mpc.finalize()

W = mpc.function_space

# Variational forms
u = ufl.TrialFunction(W)
v = ufl.TestFunction(W)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
m = ufl.inner(u, v) * ufl.dx

# Assemble matrices
A = assemble_matrix(fem.form(a), mpc)
A.assemble()
B = assemble_matrix(fem.form(m), mpc)
B.assemble()

# Slepc settings
eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setDimensions(10)   # number of eigenvalues
eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)

# Solving and recovering eigenvectors
eps.solve()
nconv = eps.getConverged()
eigvec_list = []
for i in range(nconv):
    eigval = eps.getEigenvalue(i)
    xr = A.createVecRight()
    eps.getEigenvector(i, xr)
    print("Eigenvalue:", eigval)
    eigvh = fem.Function(W, name = "eigvec{0}".format(i))
    eigvh.x.array[:] = xr[:]
    mpc.backsubstitution(eigvh)
    eigvec_list.append(eigvh)


# Plot pyvista
topology, cell_types, geometry = plot.vtk_mesh(V)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)

plotter = pv.Plotter()

grid.point_data["eigv"] = eigvec_list[8].x.array # also in 12, ...  
plotter.add_mesh(grid, scalars="eigv")

plotter.view_xy()
plotter.enable_parallel_projection()
plotter.add_scalar_bar("eigv")
plotter.show()

Dear all,

Here you have MWE of an eigenvalue problem solved by Slepc. The periodic boundary conditions are imposed by dolfinx_mpc. I’ve noticed some weird behaviour on the upper-right corner for some eigenmodes. I’m using dolfinx 0.10, serial mode. Do you any idea how to fix that? I have left in comments some wrong fixes that I tried.

2

What about enforce the constraint in one go, i.e.

def slave_boundary(x):
    return np.isclose(x[0], 1.0) | np.isclose(x[1], 1.0)


# mapping from slave -> master
def periodic_relation(x):
    out = np.zeros_like(x)
    out[0] = x[0] - np.isclose(x[0], 1.0)
    out[1] = x[1] - np.isclose(x[1], 1.0)
    return out


# Mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

# Function space
V = fem.functionspace(domain, ("CG", 1))
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, slave_boundary, periodic_relation, [])
mpc.finalize()

resulting in the complete code:

import numpy as np
from mpi4py import MPI
from dolfinx import mesh, fem, plot
import ufl
import pyvista as pv
from slepc4py import SLEPc
from dolfinx_mpc import MultiPointConstraint, assemble_matrix

def slave_boundary(x):
    return np.isclose(x[0], 1.0) | np.isclose(x[1], 1.0)


# mapping from slave -> master
def periodic_relation(x):
    out = np.zeros_like(x)
    out[0] = x[0] - np.isclose(x[0], 1.0)
    out[1] = x[1] - np.isclose(x[1], 1.0)
    return out


# Mesh
domain = mesh.create_unit_square(MPI.COMM_WORLD, 20, 20)

# Function space
V = fem.functionspace(domain, ("CG", 1))
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, slave_boundary, periodic_relation, [])
mpc.finalize()

W = mpc.function_space

# Variational forms
u = ufl.TrialFunction(W)
v = ufl.TestFunction(W)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
m = ufl.inner(u, v) * ufl.dx

# Assemble matrices
A = assemble_matrix(fem.form(a), mpc)
A.assemble()
B = assemble_matrix(fem.form(m), mpc)
B.assemble()

# Slepc settings
eps = SLEPc.EPS().create(MPI.COMM_WORLD)
eps.setOperators(A, B)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
eps.setDimensions(10)   # number of eigenvalues
eps.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)

# Solving and recovering eigenvectors
eps.solve()
nconv = eps.getConverged()
eigvec_list = []
for i in range(nconv):
    eigval = eps.getEigenvalue(i)
    eigvh = fem.Function(W, name = "eigvec{0}".format(i))
    eps.getEigenvector(i, eigvh.x.petsc_vec)
    eigvh.x.scatter_forward()
    print("Eigenvalue:", eigval)
    mpc.backsubstitution(eigvh)
    eigvec_list.append(eigvh)


# Plot pyvista
topology, cell_types, geometry = plot.vtk_mesh(W)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)

plotter = pv.Plotter()

grid.point_data["eigv"] = eigvec_list[8].x.array # also in 12, ...  
plotter.add_mesh(grid, scalars="eigv")

plotter.view_xy()
plotter.enable_parallel_projection()
plotter.add_scalar_bar("eigv")
plotter.show()

which to me yields:

image
image874×711 55.5 KB

2 Likes
3

Thanks! It has solved the issue.