# Memory use buildup when solving PDE multiple times

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["NDSolveFEM"];$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["NDSolveFEM"];
$HistoryLength = 0; myeigvals[mesh_, u_, pdecoeffs_, bc_, neig_, shift_] := With[ {space = Table[u[[1]][[j]], {j, 1, Length@u[[1]]}], dependvars = Head[#] & /@ u}, vd = NDSolveVariableData[{"DependentVariables" -> dependvars, "Space" -> space}]; sd = NDSolveSolutionData[{"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["NDSolveFEM"];$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 = NDSolveVariableData[{"DependentVariables" -> dependvars,
"Space" -> space}];
sd = NDSolveSolutionData[{"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.

• This is certainly something for @user21. ;) Feb 13, 2022 at 16:12
• I am not at a computer right now but try to generate the mesh first mesh = ToElementMesh[Disk[]] and then use the mesh in stead of the Disk. See if that reduces the memory leak. Feb 13, 2022 at 17:25
• I tried what you suggested @user21 It seems that this may indeed be reducing the memory increment slightly. The reduction seems to be slightly less than a factor of two, with typical values changing from approx. 170 000 bytes to around 100 000 (the exact value of the increment appears to vary a bit from run to run). Feb 13, 2022 at 23:05
• Please see the update, that should help. Feb 17, 2022 at 8:11
• Yes, that is very helpful @user21. Thank you! Feb 18, 2022 at 22:27

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["NDSolveFEM"]; 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 support@wolfram.com 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["NDSolveFEM"];$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"].

• Sorry, I was being unclear. The testing I reported in my previous comment was referring to the Disk[] example in my question. I will test the example in your extended comment later this evening, maybe I missed something. Feb 14, 2022 at 8:45
• Ok, I ran your example and get a similar output. Yet, when solving the PDE it seems to matter whether mesh is assigned to a variable or not. Consider ex1: $HistoryLength = 0; Needs["NDSolveFEM"]; Do[NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y} [Element] Disk[], 6]; If[Mod[n, 100] == 0, Print[MemoryInUse[]]], {n, 1000}]; which outputs 143139888 156944704 170440928 ... with an increment of about 13*10^6 bytes Feb 14, 2022 at 20:25 • and ex2:$HistoryLength = 0; Needs["NDSolveFEM"]; mesh = ToElementMesh[Disk[]]; Do[NDEigensystem[{-Laplacian[u[x, y], {x, y}], DirichletCondition[u[x, y] == 0, True]}, u[x, y], {x, y} [Element] mesh, 6]; If[Mod[n, 100] == 0, Print[MemoryInUse[]]], {n, 1000}]; which outputs 267670608 275109456 281566592 with an increment of about 7*10^6 bytes. Maybe this behavior is more obvious to you though. Either way, I will edit my question with a more complete code in a few hours. Feb 14, 2022 at 20:28
• I have edited my question to provide a more clear example of my user case. I hope this is enough for some debugging to be possible but let me know if you need more information. Feb 14, 2022 at 23:33
• I may have found another big leak when imaginary components are present in the PDE coefficients. I described this in the second edit of my original post and would greatly appreciate it if you could take a look. Feb 20, 2022 at 18:26