Eigenvalues of Laplacian with a inhomogenes boundary condition

mssh · December 12, 2022, 10:07am

Hi,

I am interested in the eigen-Problem for the Laplacian with boundary conditions which are not the homogene Dirichlet condition. In particular I tried to modify the Poisson example and got it to work with the Dirichlet boundary condition:

import numpy as np
import ufl
from dolfinx import fem, mesh, plot
from ufl import grad, dx, dot, inner
from mpi4py import MPI
from slepc4py import SLEPc
from petsc4py.PETSc import ScalarType
import matplotlib.pyplot as plt
import pyvista

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                            cell_type=mesh.CellType.triangle,)
V = fem.FunctionSpace(msh, ("Lagrange", 1))

facets = mesh.locate_entities_boundary(msh, dim=1,
                                       marker=lambda x: np.logical_or(np.logical_or(np.isclose(x[0], 0.0),
                                                                      np.isclose(x[0], 2.0)),np.logical_or(np.isclose(x[1], 0.0),
                                                                      np.isclose(x[1], 1.0))))


dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)


bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = inner(u, v) * dx


A = fem.petsc.assemble_matrix(fem.form(a), bcs=[bc])
A.assemble()

B = fem.petsc.assemble_matrix(fem.form(L), bcs=[bc])
B.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A,B)
eigensolver.setType(SLEPc.EPS.Type.JD)
eigensolver.setWhichEigenpairs(eigensolver.Which.SMALLEST_REAL)
eigensolver.setDimensions(10)
eigensolver.solve()
vr, vi = A.getVecs()
uh = fem.Function(V)
cells, types, x = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
for i in range(4,eigensolver.getConverged()):
    lmbda = eigensolver.getEigenpair(i, vr, vi)
    print(lmbda/np.pi)
    uh.vector.setArray(vr.array)
    grid.point_data["u"] = uh.x.array.real
    grid.set_active_scalars("u")
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=True)
    warped = grid.warp_by_scalar()
    plotter.add_mesh(warped)
    plotter.show()

Now if I change my boundary bc to fem.dirichletbc(value=ScalarType(1), dofs=dofs, V=V) I get the same result, the boundary condition seems to be ignored. I read up a bit, and saw that it might make sense to write u=u_0+u_g, where u_0 is equal to zero on the boundary and \forall x\in\partial \Omega : u_g(x)=g(x).

So my first question would be what do I have to do, to get a u which respects my boundary condition?

My second question would be quite similar: I wanted to solve the same problem, but with the same boundary condition as in the Poisson example. My naive approach would have been to start from the Poisson code and modifiy L=inner(f, v) * dx + inner(g, v) * ds and then add the code to solve the eigen problem right after:

import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from ufl import ds, dx, grad, inner

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from slepc4py import SLEPc
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                            cell_type=mesh.CellType.triangle,)
V = fem.FunctionSpace(msh, ("Lagrange", 1))


facets = mesh.locate_entities_boundary(msh, dim=1,
                                       marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                                                      np.isclose(x[0], 2.0)))


dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)


bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
a = inner(grad(u), grad(v)) * dx
L = inner(u, v) * dx + inner(g, v) * ds

A = fem.petsc.assemble_matrix(fem.form(a), bcs=[bc])
A.assemble()

B = fem.petsc.assemble_matrix(fem.form(L), bcs=[bc])
B.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A,B)
eigensolver.setType(SLEPc.EPS.Type.JD)
eigensolver.setWhichEigenpairs(eigensolver.Which.SMALLEST_REAL)
eigensolver.setDimensions(10)
eigensolver.solve()
vr, vi = A.getVecs()
uh = fem.Function(V)
cells, types, x = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
for i in range(9,eigensolver.getConverged()):
    lmbda = eigensolver.getEigenpair(i, vr, vi)
    print(lmbda/(np.pi))
    uh.vector.setArray(vr.array)
    grid.point_data["u"] = uh.x.array.real
    grid.set_active_scalars("u")
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=True)
    warped = grid.warp_by_scalar()
    plotter.add_mesh(warped)
    plotter.show()

But I get an error while creating B:

 File "[...]/src/eigenfunctions_laplace2.py", line 38, in <module>
    B = fem.petsc.assemble_matrix(fem.form(L), bcs=[bc])
  File "[...]/anaconda3/envs/fenicsx050/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 166, in form
    return _create_form(form)
...
ufl.algorithms.check_arities.ArityMismatch: Integrand arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None),) differ from form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 0, None), Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Lagrange', triangle, 1)), 1, None)).

How can I apply the same / similar boundary conditions to my eigen problem?

dokken · December 12, 2022, 10:17am

This term is bilinear and has to be moved to the LHS of your equation.
In general, you can write your form on residual from, i.e.

F = inner(grad(u), grad(v)) * dx - inner(u, v) * dx - inner(g, v) * ds

and use ufl.lhs and ufl.rhs (or ufl.system) to split the form into a bilinear and linear form

mssh · December 14, 2022, 4:32pm

Thank you for your answer. Do I need to change other things too?

The following (modified) code

import numpy as np

from dolfinx import fem, io, mesh, plot
import ufl
from ufl import ds, dx, grad, inner, lhs,rhs

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from slepc4py import SLEPc
msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (2.0, 1.0)), n=(32, 16),
                            cell_type=mesh.CellType.triangle,)
V = fem.FunctionSpace(msh, ("Lagrange", 1))


facets = mesh.locate_entities_boundary(msh, dim=1,
                                       marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                                                      np.isclose(x[0], 2.0)))


dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)


bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

# +
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
F = inner(grad(u), grad(v)) * dx-inner(u, v) * dx -inner(g, v) * ds
r=rhs(F)
l=lhs(F)

A = fem.petsc.assemble_matrix(fem.form(l), bcs=[bc])
A.assemble()

B = fem.petsc.assemble_matrix(fem.form(r), bcs=[bc])
B.assemble()

eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A,B)
eigensolver.setType(SLEPc.EPS.Type.JD)
eigensolver.setWhichEigenpairs(eigensolver.Which.SMALLEST_REAL)
eigensolver.setDimensions(10)
eigensolver.solve()
vr, vi = A.getVecs()
uh = fem.Function(V)
cells, types, x = plot.create_vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
for i in range(9,eigensolver.getConverged()):
    lmbda = eigensolver.getEigenpair(i, vr, vi)
    print(lmbda/(np.pi))
    uh.vector.setArray(vr.array)
    grid.point_data["u"] = uh.x.array.real
    grid.set_active_scalars("u")
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=True)
    warped = grid.warp_by_scalar()
    plotter.add_mesh(warped)
    plotter.show()

throws this error:

Traceback (most recent call last):
  File "", line 39, in <module>
    B = fem.petsc.assemble_matrix(fem.form(r), bcs=[bc])
  File "...fenicsx050/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File ".../fenicsx050/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 348, in _assemble_matrix_form
    A = _cpp.fem.petsc.create_matrix(a)
RuntimeError: Cannot create sparsity pattern. Form is not a bilinear form

So I`m not sure how I can use those commands? What do I have to do after creating lhs and rhs?

dokken · December 14, 2022, 5:09pm

r is a linear form, i.e. a vector, so Im not sure how it would make sende to add a Dirichlet condition weakly (Which in fact is an inhomogenuous Neumann condition).

mssh · December 14, 2022, 11:34pm

My thought process was that the poisson example has two boundary conditions: one neumann and one dirichlet which apply to two different sets of opposite edges of our domain. But I can see how this does not make sense for our r.

Sadly I fail to see how I can use this r for my Eigenproblem, could give me an example / some suggestions what I can do after splitting my equation to get Eigenvalues and functions?

dokken · December 15, 2022, 6:52am

A different input value to the DirichletBC would never alter the structure of the matrix you would like to solve for.

In legacy dolfin, if a dof has a Dirichlet bc, the row is set to the identity row (if assembly is used), and to the identity row/column (if assemble_system is used).
The boundary condition is then applied directly to the RHS in the first case, and using lifting in the second case (computing L-=Ag)

In DOLFINx the behavior of assembly is the same as assemble_system, and thus the matrix does not change when you change the value of g in your Dirichlet bc

mssh · December 16, 2022, 9:06am

Ok, i think I understand that. How do I continue from here? Say I drop the Dirichlet-Boundary and just use one inhomogenuous Neumann condition. How do I get from

F = inner(grad(u), grad(v)) * dx - inner(u, v) * dx - inner(g, v) * ds

and the lhs / rhs of F to solutions of the generalized Eigenproblem?

B = fem.petsc.assemble_matrix(fem.form(r))
B.assemble()

still gives the same error.

Topic		Replies	Views
Boundary conditions for a generalised eigenvalue problem Linear Algebra	5	1424	November 26, 2021
Eigenvalues Laplacian on the 2D half sphere	0	433	August 27, 2021
User defined Boundary Conditions in multidimensional problems General	4	445	December 9, 2022
Non-trivial boundary conditions for linear elasticity General	3	55	November 5, 2024
The Problem of Dirichlet Boundary Conditions for DG General	4	536	October 30, 2023

Eigenvalues of Laplacian with a inhomogenes boundary condition

Related topics