SLEPc: non-homogeneous Dirichlet BC not working

dolfinx
1

Hi everyone, I am trying to obtain the eigenvalues / egienfunctions of a simple Laplace problem in a 1D interval (0,1) using SLEPc library.

However I seem to be running into some trouble when defining the BCs. Specifically, I’d like to impose:

u(0) = 1 & u(1) = 0

with u being the solution to the Laplace problem. I tried implementing the following loop

bc_facets = exterior_facet_indices(msh.topology)
bc_dofs = fem.locate_dofs_topological(V, msh.topology.dim - 1, bc_facets)
u_bc = fem.Function(V)
with u_bc.vector.localForm() as loc:
    loc.set(0)

for dof in bc_dofs:
    x = msh.geometry.x
    if np.isclose(x[dof, 0], 0.0)==1:  # Compara si x está cerca de 0
        u_bc.vector[dof] = 1  # Valor en el extremo izquierdo
    elif np.isclose(x[dof, 0], 1.0)==1:  # Compara si x está cerca de 1
        u_bc.vector[dof] = 1  # Valor en el extremo derecho

bc = fem.dirichletbc(u_bc, bc_dofs)

But the solution I get is always that of the Laplace problem with homogeneous Dirichlet BC. How may I make sure that my boundary conditions are being imposed correctly? Thank you in advance. I attach the full code below.


import sys

from mpi4py import MPI

import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type, fem, io, plot
from slepc4py import SLEPc
import petsc4py.PETSc
from dolfinx.fem.petsc import assemble_matrix
import dolfinx.fem as fem
from dolfinx.mesh import (CellType, create_interval, exterior_facet_indices,
                         locate_entities)
from ufl import TestFunction, TrialFunction, dx, grad, inner

try:
   import pyvista
   have_pyvista = True
except ModuleNotFoundError:
   print("pyvista and pyvistaqt are required to visualise the solution")
   have_pyvista = False

# Domain definition 

print("Programa inicializado")

rh = 0
rinf = 1
N = 100

comm = MPI.COMM_WORLD
msh = create_interval(comm,N,[rh, rinf])
# Create mesh and finite element
msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)

lambda0 = 2
lambdaTg = lambda0
# Function Space

V = fem.functionspace(msh, ("Lagrange", 1))

# Weak form
u = TrialFunction(V)
v = TestFunction(V) 

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

print("Problema Variacional definido")

a = fem.form(a)
b = fem.form(b)

# Boundary Conds

bc_facets = exterior_facet_indices(msh.topology)
bc_dofs = fem.locate_dofs_topological(V, msh.topology.dim - 1, bc_facets)
u_bc = fem.Function(V)
with u_bc.vector.localForm() as loc:
   loc.set(0)

for dof in bc_dofs:
   x = msh.geometry.x
   if np.isclose(x[dof, 0], 0.0)==1:  # Compara si x está cerca de 0
       u_bc.vector[dof] = 1  # Valor en el extremo izquierdo
   elif np.isclose(x[dof, 0], 1.0)==1:  # Compara si x está cerca de 1
       u_bc.vector[dof] = 1  # Valor en el extremo derecho

bc = fem.dirichletbc(u_bc, bc_dofs)


# Construction of the solver

A = assemble_matrix(a, bcs=[bc])
A.assemble()
B = assemble_matrix(b, bcs=[bc])
B.assemble()

# Initialization of the solver

eps = SLEPc.EPS().create(msh.comm)
eps.setOperators(A, B)
eps.setProblemType(SLEPc.EPS.ProblemType.GHEP)
tol = 1e-4
eps.setTolerances(tol=tol)

# Choice of method

eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)

# Get ST context from eps
st = eps.getST()

# Set shift-and-invert transformation
st.setType(SLEPc.ST.Type.SINVERT)

# Set which part we want to recover and target value 

eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)
eps.setTarget(lambdaTg)


# How many eigenvalues?

eps.setDimensions(nev=4)
eps.solve()
eps.view()
eps.errorView()

vals = [(i, eps.getEigenvalue(i)) for i in range(eps.getConverged())]
print(vals)
# vals.sort(key=lambda x: x[1].real)

ehr = fem.Function(V)
ehi = fem.Function(V)

lambda_list = []
ehr_list = []

for i, lamb in vals:
   # Save eigenvector in eh
   eps.getEigenpair(i, ehr.vector, ehi.vector)

   # Compute error for i-th eigenvalue
   error = eps.computeError(i, SLEPc.EPS.ErrorType.RELATIVE)

   # Save and visualize 
   ehr_list.append(ehr.vector)
   lambda_list.append(lamb)

print(lambda_list)

eigenone = ehr_list[0]

import pyvista as pv
x = np.linspace(rh, rinf, N+1)
y = eigenone
chart = pv.Chart2D()

_ = chart.scatter(x, y)

_ = chart.line(x, y, 'r')

chart.show()
2

Eigen-solutions by definition relate to homogeneous PDEs. The non-homogeneous parts relate to particular solutions, and these functions can never satisfy the linear relation of the eigenvalue statement.

Indeed, this happens automatically in the linear algebra that you’re executing. At a Dirichlet BC, your assembled A and B have 1s in the diagonal, and the generalized eigenvalue problem for those rows are left to say:
Au = \lambda Bu \rightarrow u|_{BC} = \lambda u|_{BC}. For \lambda not equal to 1, the only solution is u|_{BC} = 0.

2 Likes
3

Thank you very much for your reply! I see, seems I was making a mistake on a theoretical level. My theoretical background is surely rusty, I’ll try to revise a bit :slight_smile: