Hi everyone,

I am solving Helmholtz equation with a point source in the center of a square domain. Nothing was wrong until I decided to solve it for a certain number of frequencies. It solves correctly up to about the 300th frequency of the list I provide and then it stops with the following error:

```
Traceback (most recent call last):
File "/media/Discone/Fenics/Docker/MWE.py", line 63, in <module>
problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
File "/usr/local/dolfinx-complex/lib/python3.8/dist-packages/dolfinx/fem/petsc.py", line 545, in __init__
opts[k] = v
File "PETSc/Options.pyx", line 23, in petsc4py.PETSc.Options.__setitem__
File "PETSc/Options.pyx", line 91, in petsc4py.PETSc.Options.setValue
petsc4py.PETSc.Error: error code 55
[0] PetscOptionsSetValue() at /usr/local/petsc/src/sys/objects/options.c:1165
[0] PetscOptionsSetValue_Private() at /usr/local/petsc/src/sys/objects/options.c:1215
[0] Out of memory. Allocated: 0, Used by process: 104472576
[0] Number of options 512 < max number of options 512, can not allocate enough space
```

What is happening? Why do I get this error if my system has quite a lot of free memory (about 20 GB)

this is th MWE:

```
import numpy as np
import ufl
from dolfinx import geometry
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form, Constant
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square
from ufl import dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
f_axis = np.arange(50, 2000, 5) #frequencies at which the Helmholtz equation is solved
# approximation space polynomial degree
deg = 1
# number of elements in each direction of msh
n_elem = 32
msh = create_unit_square(MPI.COMM_SELF, n_elem, n_elem)
n = ufl.FacetNormal(msh)
# Test and trial function space
V = FunctionSpace(msh, ("Lagrange", deg))
# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = Function(V)
x = ufl.SpatialCoordinate(msh)
#Source location
Sx = 0.5
Sy = 0.5
# Source amplitude
Q = 1
#Dirac
a = 0.001
delta_tmp = Function(V)
delta_tmp.interpolate(lambda x : 1/(np.abs(a)*np.sqrt(np.pi))*np.exp(-(((x[0]-Sx)**2+(x[1]-Sy)**2)/(a**2))))
int_delta_tmp = assemble_scalar(form(delta_tmp*ufl.dx))
delta = delta_tmp/int_delta_tmp
int_delta = assemble_scalar(form(delta*ufl.dx))
for nf in range(0,len(f_axis)):
freq = f_axis[nf]
print("Computing frequency: " + str(freq) + "...")
omega = freq*2*np.pi
c0 = 343
rho_0 = 1.21
k0 = 2*np.pi*freq/c0
f = 1j*rho_0*omega*Q*delta
a = inner(grad(u), grad(v)) * dx - k0**2 * inner(u, v) * dx
L = inner(f, v) * dx
# Compute solution
uh = Function(V)
uh.name = "u"
problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()
```

As you can see the list goes from 50 to 2000 Hz. The for loop stops at about 1400 Hz.

As always, thank you for the help