DOLFINX Violates Superposition Law for an Eigenvalue Problem with SLEPc

General
1

Hi there,

I am trying to solve the in-homogeneous Helmholtz equation with damping and forcing. The equation is;

c^2\nabla^2p - c^2\omega i d \nabla^2p + \omega^2 p = f

where the d is the damping parameter, the i is complex number the f is the forcing, \omega is the eigenvalue and the p is the eigenvector. In matrix form;

(\mathbf{A} + \omega \mathbf{B} + \omega^2 \mathbf{C})\mathbf{p} = \mathbf{D}\mathbf{p}

where \mathbf{B} is the damping matrix and \mathbf{D} is the forcing matrix. For 1D case, the boundary conditions are Neumann at both ends.

I am aiming to damp the system with \mathbf{B} and did try 4 different cases;

  • \mathbf{A} + \omega^2 \mathbf{C} = 0 (no damping and no forcing)
  • \mathbf{A} + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0 (damping and no forcing)
  • (\mathbf{A}-\mathbf{D}) + \omega^2 \mathbf{C} = 0 (no damping and forcing)
  • (\mathbf{A}-\mathbf{D}) + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0 (damping and forcing)

When I run this cases, I get very weird behaviour. Here is the list of the eigenvalues for the cases;

  • \mathbf{A} + \omega^2 \mathbf{C} = 0 gives \omega=1069.2399-0j
  • \mathbf{A} + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0 gives \omega=1139.6967-0.2463j
  • (\mathbf{A}-\mathbf{D}) + \omega^2 \mathbf{C} = 0 gives \omega=1074.748+3.2431j
  • (\mathbf{A}-\mathbf{D}) + \omega \mathbf{B} + \omega^2 \mathbf{C} = 0 gives \omega=1146.0758+3.5129j

The complex part of the eigenvalue is representing the growth rate of the system. As I impose a damping, the last case’s complex part (+3.5129j) should be lower than the 3. case’s complex part (3.2431j) according to the law of superposition but it is not. Here is the 1D MWE;

from dolfinx.mesh import meshtags,locate_entities,create_unit_interval
from dolfinx.fem import FunctionSpace, form, Constant
from dolfinx.fem.petsc import assemble_matrix
from ufl import Measure, TestFunction, TrialFunction, grad, inner
from petsc4py import PETSc
from slepc4py import SLEPc
from mpi4py import MPI
import numpy as np

degree = 1

mesh = create_unit_interval(MPI.COMM_WORLD, 20)

L, x_h = 0.1, 0.7
def subdomain_func(x, x_f=x_h, a_f=L/2, eps=1e-16):
    x = x[0]
    return np.logical_and(x_f - a_f - eps <= x, x <= x_f + a_f + eps)
marked_cells = locate_entities(mesh, mesh.topology.dim, subdomain_func)
forcing_tag = 2
subdomains = meshtags(mesh, mesh.topology.dim, marked_cells, np.full(len(marked_cells), forcing_tag, dtype=np.int32))

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

V = FunctionSpace(mesh, ("Lagrange", degree))
u, v = TrialFunction(V), TestFunction(V)

c = Constant(mesh, PETSc.ScalarType(340)) #speed of sound

a_form = form(-c**2* inner(grad(u), grad(v))*dx)
A = assemble_matrix(a_form)
A.assemble()

d_damp = 0.25
b_form = form(c**2 * 1j*d_damp*inner(grad(u),grad(v))*dx(2))
B = assemble_matrix(b_form)
B.assemble()

c_form = form(inner(u , v) * dx)
C = assemble_matrix(c_form)
C.assemble()

rows = [ 0,  0,  0,  4,  8, 12, 16, 20, 24, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28]
cols = [2, 3, 5, 6, 2, 3, 5, 6, 2, 3, 5, 6, 2, 3, 5, 6, 2, 3, 5, 6, 2, 3,
       5, 6, 2, 3, 5, 6]
vals = [ 2.11898752e-01+1.22975798e-01j,  1.56520670e+00+9.08370344e-01j,
       -1.56520670e+00-9.08370344e-01j, -2.11898752e-01-1.22975798e-01j,
        8.63179873e+01+5.00947893e+01j,  6.37594561e+02+3.70029077e+02j,
       -6.37594561e+02-3.70029077e+02j, -8.63179873e+01-5.00947893e+01j,
        9.73640596e+02+5.65053959e+02j,  7.19187238e+03+4.17381524e+03j,
       -7.19187238e+03-4.17381524e+03j, -9.73640596e+02-5.65053959e+02j,
        2.69712743e+03+1.56528245e+03j,  1.99225426e+04+1.15620812e+04j,
       -1.99225426e+04-1.15620812e+04j, -2.69712743e+03-1.56528245e+03j,
        9.73640596e+02+5.65053959e+02j,  7.19187238e+03+4.17381524e+03j,
       -7.19187238e+03-4.17381524e+03j, -9.73640596e+02-5.65053959e+02j,
        8.63179873e+01+5.00947893e+01j,  6.37594561e+02+3.70029077e+02j,
       -6.37594561e+02-3.70029077e+02j, -8.63179873e+01-5.00947893e+01j,
        2.11898752e-01+1.22975798e-01j,  1.56520670e+00+9.08370344e-01j,
       -1.56520670e+00-9.08370344e-01j, -2.11898752e-01-1.22975798e-01j]

D = PETSc.Mat().create(PETSc.COMM_WORLD)
D.setSizes([(V.dofmap.index_map.size_local, V.dofmap.index_map.size_global), (V.dofmap.index_map.size_local, V.dofmap.index_map.size_global)])
D.setType('mpiaij')
D.setUp()
D.setValuesCSR(rows, cols, vals, addv=PETSc.InsertMode.ADD_VALUES)
D.assemblyBegin()
D.assemblyEnd()

def eps_solver(A, C, target, nev):
    # This function solves A + (w^2) C = 0
    E = SLEPc.EPS().create(MPI.COMM_WORLD)
    C = - C
    E.setOperators(A, C)
    st = E.getST()
    st.setType('sinvert')
    E.setTarget(target)
    E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
    E.setDimensions(nev, SLEPc.DECIDE)
    E.setTolerances(1e-15)
    E.setFromOptions()
    E.solve()
    return E

def pep_solver(A, B, C, target, nev):
    # This function solves A + wB + w^2 C = 0
    Q = SLEPc.PEP().create(MPI.COMM_WORLD)
    operators = [A, B, C]
    Q.setOperators(operators)
    st = Q.getST()
    st.setType('sinvert')
    Q.setTarget(target)
    Q.setWhichEigenpairs(SLEPc.PEP.Which.TARGET_MAGNITUDE)
    Q.setDimensions(nev, SLEPc.DECIDE)
    Q.setTolerances(1e-15)
    Q.setFromOptions()
    Q.solve()
    return Q

target = 200 * 2 * np.pi

E = eps_solver(A, C, target**2, 2)
vr, vi = A.getVecs()
eig = np.sqrt(E.getEigenvalue(0))
print("Eigenvalue without damping without forcing: ", eig)

E2 = pep_solver(A, B, C, target, 2)
vr, vi = A.getVecs()
eig2 = E2.getEigenpair(0, vr, vi)
print("Eigenvalue with damping: ", eig2)

D_Mat = A -  D
E3 = eps_solver(D_Mat, C, target**2, 2)
vr, vi = A.getVecs()
eig3 = np.sqrt(E3.getEigenvalue(0))
print("Eigenvalue with forcing: ", eig3)

E4 = pep_solver(D_Mat, B, C, target, 2)
vr, vi = A.getVecs()
eig4 = E4.getEigenpair(0, vr, vi)
print("Eigenvalue with damping and forcing: ", eig4)

gives;

Eigenvalue without damping without forcing:  (1069.2399772685112-1.373743027646997e-14j)
Eigenvalue with damping:  (1139.6967684690483-0.24634561557213236j)
Eigenvalue with forcing:  (1074.7488345843487+3.2431003471216204j)
Eigenvalue with damping and forcing:  (1146.075830311377+3.512980597062334j)

This results makes me confused and any suggestions would be warmly welcomed. Thank you in advance!

2

If my problem is not clear, please let me know and I try my best to explain it.

3

Is there any way to define the matrix D in the code without hardcoding those specific entries?

At first glance, I would have wanted to double check if you get the expected results by refining the mesh (i.e., using more intervals than 20 in create_unit_interval), but that is not allowed because D is only defined for the case of mesh size equal to 20.

4

Thank you @francesco-ballarin Here is the full code, which includes the builing of the matrix D;

from dolfinx.mesh import meshtags,locate_entities,create_unit_interval
from dolfinx.fem import FunctionSpace, form, Constant, Function
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from ufl import Measure, TestFunction, TrialFunction, grad, inner, exp, as_vector
from petsc4py import PETSc
from slepc4py import SLEPc
from mpi4py import MPI
import numpy as np

degree = 1

mesh = create_unit_interval(MPI.COMM_WORLD, 100)

L, x_h = 0.1, 0.7
def subdomain_func(x, x_f=x_h, a_f=L/2, eps=1e-16):
    x = x[0]
    return np.logical_and(x_f - a_f - eps <= x, x <= x_f + a_f + eps)
marked_cells = locate_entities(mesh, mesh.topology.dim, subdomain_func)
forcing_tag = 2
subdomains = meshtags(mesh, mesh.topology.dim, marked_cells, np.full(len(marked_cells), forcing_tag, dtype=np.int32))

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

V = FunctionSpace(mesh, ("Lagrange", degree))
u, v = TrialFunction(V), TestFunction(V)

c = Constant(mesh, PETSc.ScalarType(340)) #speed of sound

a_form = form(-c**2* inner(grad(u), grad(v))*dx)
A = assemble_matrix(a_form)
A.assemble()

d_damp = 0.25
b_form = form(c**2 * 1j*d_damp*inner(grad(u),grad(v))*dx(2))
B = assemble_matrix(b_form)
B.assemble()

c_form = form(inner(u , v) * dx)
C = assemble_matrix(c_form)
C.assemble()

def gaussian(mesh, x_ref, sigma):
    def gaussian_func(x, x_ref, sigma):
        return 1/(sigma*1*(2*np.pi)**(1/2))*np.exp(-1*(x[0]-x_ref)**2/(2*sigma**2))
    V = FunctionSpace(mesh, ("CG", 1))
    g = Function(V)
    g.interpolate(lambda x: gaussian_func(x, x_ref, sigma))
    return g

coefficient = 6450
h = gaussian(mesh, 0.25, 0.025)
tau = 0.0032
omega = 1100 
left_form = form(coefficient *  u  * h * exp(1j*omega*tau)  *  dx)

rho = 1.22
w = gaussian(mesh, 0.2, 0.025)
right_form = form(inner(as_vector([1]),grad(v)) / rho * w * dx)        

def indices_and_values(form, V, tol=1E-5):
    temp = assemble_vector(form)
    temp.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    temp.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    packed = temp.array
    packed.real[abs(packed.real) < tol] = 0.0
    packed.imag[abs(packed.imag) < tol] = 0.0
    indices = np.array(np.flatnonzero(packed),dtype=np.int32)
    global_indices = V.dofmap.index_map.local_to_global(indices)
    packed = list(zip(global_indices, packed[indices]))
    return packed

left_vector = indices_and_values(left_form, V)
right_vector = indices_and_values(right_form, V)

row = [item[0] for item in left_vector]
col = [item[0] for item in right_vector]

row_vals = [item[1] for item in left_vector]
col_vals = [item[1] for item in right_vector]

product = np.outer(row_vals,col_vals)
val = product.flatten()

global_size = V.dofmap.index_map.size_global
local_size = V.dofmap.index_map.size_local

D = PETSc.Mat().create(PETSc.COMM_WORLD)
D.setSizes([(local_size, global_size), (local_size, global_size)])
D.setType('mpiaij')
D.setUp()
D.setValues(row, col, val, addv=PETSc.InsertMode.ADD_VALUES)
D.assemblyBegin()
D.assemblyEnd()

def eps_solver(A, C, target, nev):
    # This function solves A + (w^2) C = 0
    E = SLEPc.EPS().create(MPI.COMM_WORLD)
    C = - C
    E.setOperators(A, C)
    st = E.getST()
    st.setType('sinvert')
    E.setTarget(target)
    E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
    E.setDimensions(nev, SLEPc.DECIDE)
    E.setTolerances(1e-15)
    E.setFromOptions()
    E.solve()
    return E

def pep_solver(A, B, C, target, nev):
    # This function solves A + wB + w^2 C = 0
    Q = SLEPc.PEP().create(MPI.COMM_WORLD)
    operators = [A, B, C]
    Q.setOperators(operators)
    st = Q.getST()
    st.setType('sinvert')
    Q.setTarget(target)
    Q.setWhichEigenpairs(SLEPc.PEP.Which.TARGET_MAGNITUDE)
    Q.setDimensions(nev, SLEPc.DECIDE)
    Q.setTolerances(1e-15)
    Q.setFromOptions()
    Q.solve()
    return Q

target = 200 * 2 * np.pi

E = eps_solver(A, C, target**2, 2)
vr, vi = A.getVecs()
eig = np.sqrt(E.getEigenvalue(0))
print("Eigenvalue without damping without forcing: ", eig)

E2 = pep_solver(A, B, C, target, 2)
vr, vi = A.getVecs()
eig2 = E2.getEigenpair(0, vr, vi)
print("Eigenvalue with damping: ", eig2)

D_Mat = A -  D
E3 = eps_solver(D_Mat, C, target**2, 2)
vr, vi = A.getVecs()
eig3 = np.sqrt(E3.getEigenvalue(0))
print("Eigenvalue with forcing: ", eig3)

E4 = pep_solver(D_Mat, B, C, target, 2)
vr, vi = A.getVecs()
eig4 = E4.getEigenpair(0, vr, vi)
print("Eigenvalue with damping and forcing: ", eig4)

The mesh is not affecting the results much.

5

Checking for mesh dependence was the only thing that came to my mind. If that is not the issue, then I do not have any experience with SLEPc.PEP or these type of problems, and thus I am not able to comment any further.

6

@francesco-ballarin Thank you. Actually, I am not sure that if this problem occurs due to the subroutines of dolfinx rather than SLEPc.

7

What I would suggest doing is to do a really small example where you can compute your matrix by hand (say a 5 element problem). Then you can verify if you get similar matrices with DOLFINx, and if it SLEPC is the problem.

1 Like