Hi everyone,
I’m using macOS Mojave (10.14.6), miniconda, and FEniCS 2018. I want to compute the eigenfunctions/vectors of
- \nabla \cdot (\mu(u) \nabla \phi_m) = \lambda_m \phi_m
for m = 1, \ldots, K with zero Dirichlet boundary condition, where \mu is defined as
\mu(u) := \frac{1}{\sqrt{|\nabla u|^2 + \epsilon^2}}
for an \epsilon > 0 (normally set as \epsilon = 10^{-6} ). I want u to be 0 on the boundary and piecewise constant else.
I chose u as
u1 = Expression('x[0] >= 0.4 && x[0] <= 0.6 && x[1] >= 0.3 && x[1] <= 0.4 ? 1.5 : 0', degree = 2)
u2 = Expression('x[0] >= 0.65 && x[0] <= 0.8 && x[1] >= 0.65 && x[1] <= 0.8 ? 1.5 : 0', degree = 2)
u = u1 + u2
u = project(u, V)
The error I get is:
*** -------------------------------------------------------------------------
*** 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.
*** Process: 0
***
*** DOLFIN version: 2018.1.0
*** Git changeset:
*** -------------------------------------------------------------------------
And here comes the funny part: If I only take u1 as u my code works. If I set \mu = 1 (then I have the eigenfunctions/vector of the Laplacian), it also works with u = u1 + u2.
Here is my code:
from fenics import *
import numpy as np
N = 50
epsilon = 1e-6
K = 2
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, 'P', 1)
u1 = Expression('x[0] >= 0.4 && x[0] <= 0.6 && x[1] >= 0.3 && x[1] <= 0.4 ? 1.5 : 0', degree = 2)
u2 = Expression('x[0] >= 0.65 && x[0] <= 0.8 && x[1] >= 0.65 && x[1] <= 0.8 ? 1.5 : 0', degree = 2)
u = u1 + u2
u = project(u, V)
u_D = Expression('0*x[0]*x[1]', degree = 2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# variable coefficient
def mu(z, epsilon):
return 1/(sqrt(dot(grad(z), grad(z))**2 + epsilon**2))
# define variational problem
phi = TrialFunction(V)
w = TestFunction(V)
a = mu(u, epsilon)*dot(grad(phi), grad(w))*dx
m = phi*w*dx
# assemble matrix and apply boundary conditions
A = PETScMatrix()
assemble(a, tensor = A)
bc.apply(A)
M = PETScMatrix()
assemble(m, tensor = M)
bc.apply(M)
# specify eigensolver properties (smalles magnitued till K-th eigenvalue/vector)
eigensolver = SLEPcEigenSolver(A, M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.solve(K)
# extract first eigenvector/pair
r, c, rx, cx = eigensolver.get_eigenpair(0)
What I tried is to calculate mu directly with mu = ... and then mu = project(mu, V) and used this in the bilinear form a = mu*dot(grad(phi), grad(w))*dx. So the problem should not lie in u and also not in mu(u, epsilon) (?) and the problem will possibly be with the computation of the eigenvectors? But with the Laplacian, it works.
I hope I formulated everything as clear as possible. Perhaps I only have overseen some small detail in the code which crashes the SLEPcEigenSolver? I’m only a few days in learning FEniCS/Python
I’m happy for any hints/tips I can get.
edit: spelling; the last code did not point out the problem