Hi,
I’m quite new to FEniCS and I’m experiencing a weird memory issue trying to solve a time-dependent system of two PDEs in the form (note: I’m avoiding the complete formulation for brevity; I can provide it in later comments if needed):
on a fairly large mesh (num_vertices = 1030301).
I’m using Linux (Ubuntu Mate 20.04) and I have 16 GB RAM.
Minimal report of the issue
Consider the following variational formulation, which I already succesfully tested in 2d (i.e. the issue is not in the definition of the variational form):
u = fenics.Function(V)
phi, sigma = fenics.split(u)
F = (
(((phi - phi_prec) / dt) * v1 * dx)
+ (lam * dot(grad(phi), grad(v1)) * dx)
+ ((1 / tau) * df_dsigma(phi) * v1 * dx)
+ (- chi * sigma * v1 * dx)
+ (A * phi * v1 * dx)
) + (
(((sigma - sigma_prec) / dt) * v2 * dx)
+ (eps * dot(grad(sigma), grad(v2)) * dx)
+ (- s * v2 * dx)
+ (delta * phi * v2 * dx)
+ (gamma * sigma * v2 * dx)
)
Where df_dsigma is a non-linear term defined as:
def df_dsigma(var_phi):
c1 = fenics.Constant(32)
c2 = fenics.Constant(-96)
c3 = fenics.Constant(64)
return c1 * var_phi + c2 * (var_phi ** 2) + c3 * (var_phi ** 3)
In my script I’m trying to solve it using mumps, as suggested in some topics of the old FEniCS forum:
solve(F == 0, u, form_compiler_parameters={"linear_solver": "mumps"})
And I get an out of Memory error:
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 6.032e+09 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
UMFPACK V5.7.8 (Nov 9, 2018): ERROR: out of memory
Traceback (most recent call last):
File "/home/francop/PycharmProjects/lorenzo-fenics/demonstration-RAM-error.py", line 120, in <module>
fenics.solve(F == 0, u, form_compiler_parameters={"linear_solver": "mumps"})
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 314, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside /build/dolfin-cZGj9S/dolfin-2019.2.0~git20200629.946dbd3/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
So I though that my mesh was too large for my ram capability. However, when I tried to monitor the memory usage using tracemalloc I get that the memory used from the program is not high at all:
Current memory usage is 0.236243MB; Peak was 0.270883MB
Do you know what might be causing the issue? Thank you in advance.
Complete script reproducing the issue
import fenics
import tracemalloc
# monitor malloc
tracemalloc.start()
# macro
DAYSINYEAR = 356 # defined for setting proper dimension to some parameters
# define files to store values
data_folder_name = "data"
vtk_file_phi = fenics.File(data_folder_name + "/phi.pvd")
vtk_file_sigma = fenics.File(data_folder_name + "/sigma.pvd")
# define simulation time
t_final = 1
n_steps = 1000
delta_t = t_final / n_steps
dt = fenics.Constant(delta_t)
# define parameters
lam = fenics.Constant(1.6E5) # um^2 / y
tau = fenics.Constant(0.01) # years
chi = fenics.Constant(600.0) # L / (g * y)
A = fenics.Constant(600.0) # 1 / years
eps = fenics.Constant(5.0E6) # um^2 / y
s_days = 2.70 # (g / L * d)
s = fenics.Constant(DAYSINYEAR * s_days) # g / (L * y)
delta_days = 2.75 # g / (L * d)
delta = fenics.Constant(DAYSINYEAR * delta_days)
gamma = 1000.0 # 1 / year
# define mesh
nx = ny = nz = 100
x_up = 1000 # um
x_down = -1000 # um
y_up = x_up
y_down = x_down
z_up = x_up
z_down = x_down
mesh = fenics.BoxMesh(fenics.Point(x_down, y_down, z_down), fenics.Point(x_up, y_up, z_up), nx, ny, nz)
# define initial condition
semiax_x = 100 # um
semiax_y = 150 # um
semiax_z = semiax_x
phi0_max = 1
phi0_min = 0
sigma0_max = 1
sigma0_min = 0.2
is_inside_ellipsoid_ccode = "((pow(x[0] / semiax_x, 2) + pow(x[1] / semiax_y, 2) + pow(x[2] / semiax_z, 2)) <= 1)"
phi0_ccode = is_inside_ellipsoid_ccode + " ? phi0_max : phi0_min"
sigma0_ccode = is_inside_ellipsoid_ccode + " ? sigma0_min : sigma0_max"
u0 = fenics.Expression((phi0_ccode, sigma0_ccode),
degree=2,
semiax_x=semiax_x, semiax_y=semiax_y, semiax_z=semiax_z,
phi0_max=phi0_max, phi0_min=phi0_min,
sigma0_max=sigma0_max, sigma0_min=sigma0_min)
# define function space
P1 = fenics.FiniteElement("P", fenics.tetrahedron, 1)
element = fenics.MixedElement([P1, P1])
V = fenics.FunctionSpace(mesh, element)
# define test functions
v1, v2 = fenics.TestFunctions(V)
# define trial functions
u = fenics.Function(V)
phi, sigma = fenics.split(u)
# define u_prec (initially assigned to u0)
u_prec = fenics.project(u0, V, mesh=mesh, solver_type="gmres", preconditioner_type="ilu")
phi_prec, sigma_prec = fenics.split(u_prec)
# define variational form
def df_dsigma(var_phi):
c1 = fenics.Constant(32)
c2 = fenics.Constant(-96)
c3 = fenics.Constant(64)
return c1 * var_phi + c2 * (var_phi ** 2) + c3 * (var_phi ** 3)
F = (
(((phi - phi_prec) / dt) * v1 * fenics.dx)
+ (lam * fenics.dot(fenics.grad(phi), fenics.grad(v1)) * fenics.dx)
+ ((1 / tau) * df_dsigma(phi) * v1 * fenics.dx)
+ (- chi * sigma * v1 * fenics.dx)
+ (A * phi * v1 * fenics.dx)
) + (
(((sigma - sigma_prec) / dt) * v2 * fenics.dx)
+ (eps * fenics.dot(fenics.grad(sigma), fenics.grad(v2)) * fenics.dx)
+ (- s * v2 * fenics.dx)
+ (delta * phi * v2 * fenics.dx)
+ (gamma * sigma * v2 * fenics.dx)
)
# store initial values in files
u.assign(u_prec)
phi_curr, sigma_curr = u.split()
vtk_file_phi << (phi_curr, 0)
vtk_file_sigma << (sigma_curr, 0)
# time integration
t = 0
for n in range(n_steps):
# update time
t += delta_t
# solve problem for current time
try:
fenics.solve(F == 0, u, form_compiler_parameters={"linear_solver": "mumps"})
except RuntimeError:
current, peak = tracemalloc.get_traced_memory()
print(f"Current memory usage is {current / 10 ** 6}MB; Peak was {peak / 10 ** 6}MB")
tracemalloc.stop()
raise
# save current solutions to file
phi_curr, sigma_curr = u.split()
vtk_file_phi << (phi_curr, t)
vtk_file_sigma << (sigma_curr, t)
# update u_prec
u_prec.assign(u)
# log
msg = f"Missing steps: {n_steps - (n + 1)}"
print(msg)