Eigensolver fails when I apply boundary conditions


119
views
0
3 months ago by

Hi,

I am following an eigensolver tutorial for FEniCS. The script below works if I comment out the lines "bc2D.apply(A) bc2D.zero(B)" but I wish to apply a boundary condition of 0 on the border. Is this possible? The script is below

from __future__ import print_function
from dolfin import *

# Define mesh, function space
n = 25
mesh = RectangleMesh(Point(0,0),Point(2.0,2.0),n,n,"left/right")
V = FunctionSpace(mesh, "Lagrange", 1)

# Define basis and bilinear forms
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
b = u*v*dx

def u0_boundary(x,on_boundary):
    return on_boundary
bc2D = DirichletBC(V,Constant(0.0),u0_boundary)

# Assemble stiffness forms
A = PETScMatrix()
B = PETScMatrix()
assemble(a,tensor=A)
#bc2D.apply(A)
assemble(b,tensor=B)
#bc2D.zero(B)

# Create eigensolver
eigensolver = SLEPcEigenSolver(A, B)
eigensolver.parameters["spectrum"]="smallest magnitude"

# Compute all eigenvalues of A x = \lambda x
print("Computing eigenvalues. This can take a minute.")
eigensolver.solve(2)

# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(1)
print("eigenvalue = ", r)

# Initialize function and assign eigenvector
u = Function(V)
u.vector()[:] = rx

# Plot eigenfunction
File('vec.pvd') << u
Community: FEniCS Project
If you leave out the DirichletBC in the code above, does that imply Neumann BC by default?
written 9 weeks ago by Singularity  

1 Answer


2
3 months ago by
I've had some problems with this in the past as well, that no amount of tinkering with the SLEPcEigenSolver parameters could resolve.  It seems that setting options directly through the slepc4py API can be more effective (although I'm not sure why).  Try something like the following:

from __future__ import print_function
from dolfin import *

# Define mesh, function space
n = 25
mesh = RectangleMesh(Point(0,0),Point(2.0,2.0),n,n,"left/right")
V = FunctionSpace(mesh, "Lagrange", 1)

# Define basis and bilinear forms
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
b = u*v*dx

def u0_boundary(x,on_boundary):
    return on_boundary
bc2D = DirichletBC(V,Constant(0.0),u0_boundary)

# Assemble stiffness forms
A = PETScMatrix()
B = PETScMatrix()

assemble(a,tensor=A)

# Trick to shift eigenvalues from BCs to the high end of the spectrum.
# (I'm not sure whether this is considered the best way to skip over the BC
# eigenvalues; setting the diagonal values too large produces weird results.)
bc2D.zero(A)
bc2D.zero_columns(A,Function(V).vector(),1.0e10)

assemble(b,tensor=B)
bc2D.apply(B)


# Adapted from MiroK's answer to this question in the old forum:
#
#   https://fenicsproject.org/qa/10134/slepc-no-longer-works-in-parallel/
#
from petsc4py import PETSc
from slepc4py import SLEPc

solver = SLEPc.EPS().create()
solver.setOperators(as_backend_type(A).mat(), as_backend_type(B).mat())
solver.setType(solver.Type.KRYLOVSCHUR)
#solver.setType(solver.Type.LAPACK)
solver.setDimensions(10, PETSc.DECIDE)
solver.setWhichEigenpairs(solver.Which.SMALLEST_MAGNITUDE)
solver.setTolerances(1E-8, 1000)

solver.solve()

nconv = solver.getConverged()
for i in range(nconv):
    fr = Function(V)
    fi = Function(V)
    k = solver.getEigenpair(i,
                            as_backend_type(fr.vector()).vec(),
                            as_backend_type(fi.vector()).vec())
    File("evecs/evec"+str(i)+".pvd") << fr
    print(k.real)
I've also seen some weird behavior where eigensolver convergence depends on the SLEPc version.  The above runs correctly for me in the 2017.2 Docker container, though.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »