Out of memory error in frequency for loop

Errors
1

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

2

It is rarely advicable to create solvers inside loops in Python. I would suggest that you rewrite your code such that k0 and omega is defined as dolfinx.fem.Constant which you in turn can assign values to.
Then you can move the definition of your linear problem and simply assign new values to the constants.

1 Like
3

How can I put the LinearProblem outside of the for loop if its input needs to change? I mean, how defining k0 and omega in with fem.constant allows me to put it outside of the for loop?

4

dolfinx.fem.Constant wraps constant quantities which can change between assemblies of the underlying matrix.

dolfinx.fem.petsc.LinearProblem has a PETSc Mat and Vec. Creating one every time step is extremely expensive, and not their intended use.

...  # Problem setup
k = dolfinx.fem.Constant(mesh, 0.0)
k.value = some_important_value
...  # Weak form definition
problem = LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
...  # Other code
for nf in range(0,len(f_axis)):
    ...  # More code
    k.value = another_important_value
    problem.solve()
    ... # Postprocess
2 Likes
5

Thank you very much, I was thinking that part of the assembly and calculations were done in LinearProblem. I Will deepen my knowledge about dolfinx.

I will share the results and the comparison between fenics and Actran soon.

6

Hi Nate, I have a similar issue but I refine my mesh while iterating.
Is there a way to keep the same solver if the mesh changes? Is it sufficient to manually set the Constant’s mesh?

Edit : I didn’t realise this topic is from may '22 and not from may '23: if it’s an issue, I can create a new one.

7

If you do a multiple refinements and solves within the same code, you might face issues with the python garbage collector (as it collects garbage whenever it feels like it).

Also, without the actual code, it’s very hard to give you any pointers.

8

Thanks for your reply: here is a sketch of my code

... # things before
while error_estimator > tol :
        

            
        msh.topology.create_entities(1)
        
        msh=mesh.refine(msh)
        
        
        a, L, bc = assemble_system( msh) # here I redefine my FE, function space and create the forms
        

        problem = fem.petsc.LinearProblem(a, L,  bcs=[bc], petsc_options=options)
  
        new_solution = problem.solve()

        ... # other things, estimate error and then iterate

I am quite sure that this is not optimal, but I have no better ideas.

9

How many iterations does this do before it fails?
Is the problem 2D or 3D?
Are you running the code in serial or parallel?
How much memory do you have on your system?

10

I do ~7 iterations starting from a coarse mesh and for one run it’s fine, but thing is I would like to run the code above multiple times for different meshes (as for different time steps) : if I do, it takes 4-5 runs with ~7 iterations each.
The problem is 2D.
I’m running serial, on jupyter notebook currently.
I have 32 GB RAM

11

4-5 refinement means 256-1024 factor in number of cells in your mesh. You should carefully construct a representative example that fail for you, that others can reproduce to Get further help.