Memory use buildup when solving PDE multiple times

Question

Bug introduced in 13.0 or earlier and fixed in 13.1.0

---

Mathematica V 13.0 Mac OS X ARM (64-bit) - MacBook Pro 13" M1 processor macOS Big Sur 11.6.2

I am solving a PDE numerically and have noticed that memory use keeps piling up if I execute the same code multiple times, even when using $HistoryLength = 0. This is becoming a serious limitation for me when carrying out my analysis since I am running parameter sweeps for no apparently good reason (the memory corresponding to the output result is very small)

While debugging I noticed that the following simple example piles up memory:

I initialize the notebook on a fresh kernel with

ClearAll["Global'*"]; Needs["NDSolve`FEM`"]; $HistoryLength = 0; 

Then I run the following example code from the NDEigensystem MM page multiple times:

{vals, funs} = NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y} \[Element] Disk[], 6]; MemoryInUse[] 

When i run the code multiple times (either the whole notebook or the last block) I noted that memory keeps piling up, i.e. MemoryInUse[] outputs 161822960, 162029368, 162202344, ...

The only way I have found to release this memory is to restart the Kernel. I am hoping that understanding this minimal example will help me understand why memory keeps piling up in my application. Can someone help me understand the reason for the memory buildup, whether it is a bug, or if there is a way to prevent the buildup / flush the memory (other than kernel restart)?

EDIT: Based on the discussion in user21's answer below, I have added more details. My process flow is as follows.

I define my own function to assemble the PDE system via the finite element method as described in the MM FEM programming guide (the reason for this is related to a current bug that affects the real PDE I am interested in solving (Inconsistent behaviour between Active and Inactive forms of PDEs (finite element method)). In this case, I am only interested in the eigenvalues:

ClearAll["Global'*"]; Needs["NDSolve`FEM`"]; $HistoryLength = 0; myeigvals[mesh_, u_, pdecoeffs_, bc_, neig_, shift_] := With[ {space = Table[u[[1]][[j]], {j, 1, Length@u[[1]]}], dependvars = Head[#] & /@ u}, vd = NDSolve`VariableData[{"DependentVariables" -> dependvars, "Space" -> space}]; sd = NDSolve`SolutionData[{"Space" -> mesh}]; methodData = InitializePDEMethodData[vd, sd, Method -> {"FiniteElement", "IntegrationOrder" -> 4, "InterpolationOrder" -> Table[u[[j]] -> 2, {j, 1, Length@u}]}]; initCoeffs = InitializePDECoefficients[vd, sd, pdecoeffs]; initBCs = InitializeBoundaryConditions[vd, sd, bc]; discretePDE = DiscretizePDE[initCoeffs, methodData, sd]; discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]; {stiffness = discretePDE["StiffnessMatrix"], load = discretePDE["LoadVector"], mass = discretePDE["MassMatrix"]}; DeployBoundaryConditions[{load, stiffness, mass}, discreteBCs, "ConstraintMethod" -> "Remove"]; Eigenvalues[{stiffness, mass}, neig, Method -> {"Arnoldi", "Shift" -> shift}] ]; 

Then I set up the problem I wish to solve (to illustrate I am using a simpler PDE, but the flow is the same):

mycoeffs[var_, param_] = {"DiffusionCoefficients" -> {{-param*IdentityMatrix[2]}}, "ReactionCoefficients" -> {{var^2}}, "MassCoefficients" -> IdentityMatrix[1]}; mesh = ToElementMesh[Rectangle[{0., 0.}, {50, 10}], "MeshElementType" -> TriangleElement, MaxCellMeasure -> 0.5]; u = {v[x, y]}; bc = Table[{DirichletCondition[u[[j]] == 0., True]}, {j, 1, Length@u}]; 

Now before I describe my user case, I decided to do the following checks. Check 1, repeatedly solve for eigenvalues, but without varying parameters.

shift = 0; neig = 10; Do[myeigvals[mesh, u, mycoeffs[0.01, 1], bc, neig, shift]; If[Mod[n, 100] == 0, Print[MemoryInUse[]]], {n, 1000}]; 

which outputs 151897904,151897904,151897904, ... (no increment)

Check 2, same as above, but vary the variable in the coefficient in each iteration:

shift = 0; neig = 10; Do[myeigvals[mesh, u, mycoeffs[0.0001*n, 1], bc, neig, shift]; If[Mod[n, 100] == 0, Print[MemoryInUse[]]], {n, 1000}]; 

which outputs: 158363160, 163791352, 169475336, ...

I'm not sure if the above is relevant to the following, but the actual application in which I become memory limited (for a larger version of the same problem) corresponds to finding the variable var which minimizes the smallest eigenvalue for a range of different parameters param. This means that I need to repeatedly evaluate the following function:

minfun[x_, p_] := myeigvals[mesh, u, mycoeffs[x, p], bc, 1, 0][[1]]; 

which is minimized w.r.t. x for many different values of p:

mem1 = MemoryInUse[]; mintable = Table[NMinimize[minfun[x, p], x], {p, Range[1, 2, 0.0005]}]; mem2 = MemoryInUse[]; mem2 - mem1 ByteCount[mintable] 

which for me outputs 120718952, 400576 i.e. whatever leak is occuring due to repeated evaluations (or at least I assume so) during NMinimize is building up memory and for my bigger version of the problem this amounts to several gigabytes which is not getting released. I hope this less obscure statement clarifies my issue.

EDIT 2: I may have found another big memory leak during the DiscretizePDE step when imaginary coefficients are involved. See the following example. I set up the same problem as before but now with the clear cache fix suggested by user21 and only go to the discretization step now and don't actually solve for the eigenvalues since I have narrowed down the leak to the discretization step.

ClearAll["Global'*"]; Needs["NDSolve`FEM`"]; $HistoryLength = 0; mydiscretizedpde[mesh_, u_, pdecoeffs_, bc_, neig_, shift_] := With[ {space = Table[u[[1]][[j]], {j, 1, Length@u[[1]]}], dependvars = Head[#] & /@ u}, ClearSystemCache["Symbolic"]; ClearFEMKernelCache[]; vd = NDSolve`VariableData[{"DependentVariables" -> dependvars, "Space" -> space}]; sd = NDSolve`SolutionData[{"Space" -> mesh}]; methodData = InitializePDEMethodData[vd, sd, Method -> {"FiniteElement", "IntegrationOrder" -> 4, "InterpolationOrder" -> Table[u[[j]] -> 2, {j, 1, Length@u}]}]; initCoeffs = InitializePDECoefficients[vd, sd, pdecoeffs]; initBCs = InitializeBoundaryConditions[vd, sd, bc]; discretePDE = DiscretizePDE[initCoeffs, methodData, sd]; ]; mycoeffs[var_, param_] = {"DiffusionCoefficients" -> {{-param*IdentityMatrix[2]}}, "ReactionCoefficients" -> {{var^2}}, "MassCoefficients" -> IdentityMatrix[1]}; mesh = ToElementMesh[Rectangle[{0., 0.}, {50, 10}], "MeshElementType" -> TriangleElement, MaxCellMeasure -> 0.1]; u = {v[x, y]}; bc = Table[{DirichletCondition[u[[j]] == 0., True]}, {j, 1, Length@u}]; 

Now execute the discretization for both imaginary and real values of param. Real case

shift = 0; neig = 1; Do[mydiscretizedpde[mesh, u, mycoeffs[0.1, 1], bc, neig, shift]; If[Mod[n, 100] == 0, Print[MemoryInUse[]]], {n, 1000}]; 

Outputs: 157649888, 157673952, 157697952 ...

Imaginary case

shift = 0; neig = 1; Do[mydiscretizedpde[mesh, u, mycoeffs[0.1, I], bc, neig, shift]; If[Mod[n, 100] == 0, Print[MemoryInUse[]]], {n, 1000}]; 

Outputs: 172341992, 185002408, 197662760 ...

The case with real coefficient leaks something like 200 bytes per execution while the case with the imaginary coefficient is 130 000 bytes per execution. This later leak appears to be important in my application since the memory usage is still high (reduced only by a factor of 2 or so) after the cache clear in my application which has imaginary components.

Answer 1

This is not an answer but an extended comment. I see the memory leak in NDEigensystem. However, I do not see a leak in ToElementMesh.

When I run this in a fresh kernel:

$HistoryLength = 0; Needs["NDSolve`FEM`"]; Do[ToElementMesh[Disk[]]; If[Mod[n, 100] == 0, Print[MemoryInUse[]]], {n, 1000}]; 

I get

115539552

115539280

115539280

115539280

115539552

115539552

115539280

115539280

115539280

115539280

This suggests to me that for this geometry there is no leak in ToElementMesh.

Since your comment suggests that generating the mesh outside of NDEigensystem has an effect on the leak, I assume that your geometry is not a Disk[]. In that case I suggest that you report this to the [email protected] and include an example of the actual code.

I am a ware of a few (small) leaks in the FEM stuff, but those are from functions I do not have any control over. Nevertheless, if I could see your real application I can try to track things down and hopefully have them fixed.

Update:

I finally have a solution for you. You'd need to clear some caches.

ClearAll["Global'*"]; Needs["NDSolve`FEM`"]; $HistoryLength = 0; mesh = ToElementMesh[Disk[]]; Do[ ClearSystemCache["Symbolic"]; ClearFEMKernelCache[]; NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y} \[Element] mesh, 6]; If[Mod[n, 10] == 0, Print[MemoryInUse[]/1024.^2]];, {n, 100}] 119.895 119.897 119.9 119.902 119.904 119.907 119.909 119.911 119.913 119.916 

There is still a tiny leak - but nothing on the scale that we had seen before. This tiny leak is fixed in the development version (13.1)

I'll have to investigate further if this FEM caching can be optimized, otherwise I'll document the clean up function in some way or another.

Update 2

From version 13.1 onward you will be able to use ClearSystemCache[] which will take care of all of this. To specifically clear the FEM cache you can then use ClearSystemCache["FiniteElementMethod"].