apply boundary condition on a generalized eigenvalue solver

Question

Is there a simple way to apply boundary condition on a generalized eigenvalue solver?

Here is what I did to apply a zero dirichlet boundary condition on a generalized eigenvalue solver:

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
u = TrialFunction(V)
v = TestFunction(V)

c1s=0.5*inner(nabla_grad(u),nabla_grad(v))*dx+u*v2s*v*dx-u*v_ext_corr*v*dx
m=u*v*dx
b=v*dx

C1s, _ = assemble_system(c1s, b, bc)
M, _ = assemble_system(m, b, bc)
bc.zero(M)
C1s = down_cast(C1s)
M = down_cast(M)
eigensolver = SLEPcEigenSolver(C1s, M)

cevito · Answer 1 · May 15, 2014

Don't use assemble_system in this case, but just assemble and apply the boundary conditions yourself, i.e.

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
u = TrialFunction(V)
v = TestFunction(V)

c1s=0.5*inner(nabla_grad(u),nabla_grad(v))*dx+u*v2s*v*dx-u*v_ext_corr*v*dx
m=u*v*dx

C1s = PETScMatrix()
assemble(c1s, tensor=C1s)
bc.apply(C1s)

M = PETScMatrix()
assemble(m, tensor=M)
bc.zero(M)

eigensolver = SLEPcEigenSolver(C1s, M)

Also have a look at the eigensolver tutorial shipping with fenics:

answered May 15, 2014 by cevito FEniCS User (5,550 points)
edited May 27, 2014 by cevito

By the method you suggested, the boundary value of the eigenvector calculated is not zero strictly. By my method, the boundary value of the eigenvector calculated is zero strictly. The eigenvalues by these two methods has a little bit difference. How do you think about it?

commented May 27, 2014 by vincehouhou FEniCS Novice (540 points)

And if the eigenvalue is 1, your method will not give zero boundary condition.

commented May 27, 2014 by vincehouhou FEniCS Novice (540 points)

In order to make the boundary condition as zero exactly.

Based on your method, I use "bc.zero(M)" instead of "bc.apply(M)".

Then I get exactly same results as my code.

commented May 27, 2014 by vincehouhou FEniCS Novice (540 points)

If I want to apply non-zero dirichlet boundary condition, what should I do?

commented May 27, 2014 by vincehouhou FEniCS Novice (540 points)

You're correct, this should be bc.zero(M). Note, that, in general, you can't set arbitrary non-zero dirichlet boundary condition in case of an eigenvalue problem as this problem is overdetermined.

commented May 27, 2014 by cevito FEniCS User (5,550 points)
edited May 27, 2014 by cevito

Thanks for your reply.

Why will this problem be overdetermined if we set arbitrary non-zero dirichlet boundary condition?
After I set the zero-boundary condition, sometimes I will meet such error:

Computing eigenvalues. This can take a minute.
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Detected zero pivot in LU factorization
see http://www.mcs.anl.gov/petsc/petsc-as/documentation/faq.html#ZeroPivot!
[0]PETSC ERROR: Zero pivot row 0 value 0 tolerance 1e-12!
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Release Version 3.2.0, Patch 5, Sat Oct 29 13:45:54 CDT 2011
[0]PETSC ERROR: See docs/changes/index.html for recent updates.
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[0]PETSC ERROR: See docs/index.html for manual pages.
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Unknown Name on a linux-gnu named vincehouhou-ThinkPad-X230 by vincehouhou Tue May 27 10:44:06 2014
[0]PETSC ERROR: Libraries linked from /build/buildd/petsc-3.2.dfsg/linux-gnu-c-opt/lib
[0]PETSC ERROR: Configure run at Tue Jun 12 21:14:38 2012
[0]PETSC ERROR: Configure options --with-shared-libraries --with-debugging=0 --useThreads 0 --with-clanguage=C++ --with-c-support --with-fortran-interfaces=1 --with-mpi-dir=/usr/lib/openmpi --with-mpi-shared=1 --with-blas-lib=-lblas --with-lapack-lib=-llapack --with-blacs=1 --with-blacs-include=/usr/include --with-blacs-lib="[/usr/lib/libblacsCinit-openmpi.so,/usr/lib/libblacs-openmpi.so]" --with-scalapack=1 --with-scalapack-include=/usr/include --with-scalapack-lib=/usr/lib/libscalapack-openmpi.so --with-mumps=1 --with-mumps-include=/usr/include --with-mumps-lib="[/usr/lib/libdmumps.so,/usr/lib/libzmumps.so,/usr/lib/libsmumps.so,/usr/lib/libcmumps.so,/usr/lib/libmumps_common.so,/usr/lib/libpord.so]" --with-umfpack=1 --with-umfpack-include=/usr/include/suitesparse --with-umfpack-lib="[/usr/lib/libumfpack.so,/usr/lib/libamd.so]" --with-cholmod=1 --with-cholmod-include=/usr/include/suitesparse --with-cholmod-lib=/usr/lib/libcholmod.so --with-spooles=1 --with-spooles-include=/usr/include/spooles --with-spooles-lib=/usr/lib/libspooles.so --with-hypre=1 --with-hypre-dir=/usr --with-ptscotch=1 --with-ptscotch-include=/usr/include/scotch --with-ptscotch-lib="[/usr/lib/libptesmumps.so,/usr/lib/libptscotch.so,/usr/lib/libptscotcherr.so]" --with-hdf5=1 --with-hdf5-dir=/usr
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: MatPivotCheck_none() line 499 in src/mat/impls/aij/seq//build/buildd/petsc-3.2.dfsg/include/private/matimpl.h
[0]PETSC ERROR: MatPivotCheck() line 518 in src/mat/impls/aij/seq//build/buildd/petsc-3.2.dfsg/include/private/matimpl.h
[0]PETSC ERROR: MatLUFactorNumeric_SeqAIJ() line 565 in src/mat/impls/aij/seq/aijfact.c
[0]PETSC ERROR: MatLUFactorNumeric() line 2876 in src/mat/interface/matrix.c
[0]PETSC ERROR: PCSetUp_LU() line 160 in src/ksp/pc/impls/factor/lu/lu.c
[0]PETSC ERROR: PCSetUp() line 819 in src/ksp/pc/interface/precon.c
[0]PETSC ERROR: KSPSetUp() line 260 in src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: PCSetUp_Redundant() line 176 in src/ksp/pc/impls/redundant/redundant.c
[0]PETSC ERROR: PCSetUp() line 819 in src/ksp/pc/interface/precon.c
[0]PETSC ERROR: KSPSetUp() line 260 in src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: STSetUp_Shift() line 94 in src/st/impls/shift/shift.c
[0]PETSC ERROR: STSetUp() line 269 in src/st/interface/stsolve.c
[0]PETSC ERROR: EPSSetUp() line 155 in src/eps/interface/setup.c
[0]PETSC ERROR: EPSSolve() line 122 in src/eps/interface/solve.c
Traceback (most recent call last):
File "har1D_test.py", line 98, in
r, c, rx, cx = eigensolver.get_eigenpair(0)
File "/usr/lib/python2.7/dist-packages/dolfin/cpp.py", line 5436, in get_eigenpair
lr, lc = self._get_eigenpair(r_vec, c_vec, i)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** https://answers.launchpad.net/dolfin

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to extract eigenpair from SLEPc eigenvalue solver.
*** Reason: Requested eigenpair (0) has not been computed.
*** Where: This error was encountered inside SLEPcEigenSolver.cpp.
*** -------------------------------------------------------------------------

commented May 27, 2014 by vincehouhou FEniCS Novice (540 points)

If we set "bc.zero(M)", M should be always singular.

Should "Detected zero pivot in LU factorization" error be usual? How should we handle this.

commented May 27, 2014 by vincehouhou FEniCS Novice (540 points)

I face the same problem of zero pivot. Did you find a solution ?

commented Dec 16, 2014 by praveen FEniCS User (2,760 points)

Try to use shift-invert, i.e.,

eigensolver.parameters["spectrum"] = "target magnitude"
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectral_shift"] = <insert_value_for_shift>

commented Dec 17, 2014 by cevito FEniCS User (5,550 points)

Hello cevito, even with these options I get zero pivot error. If I dont apply bc then it runs but of course solution is wrong. Applying the bc gives me zero pivot error.

commented Dec 17, 2014 by praveen FEniCS User (2,760 points)

apply boundary condition on a generalized eigenvalue solver

1 Answer