Memory Problem with Modules and NDSolve

Question

I'm using NDSolve inside a module, and I appear to have a memory leak. The relevant code is:

  mpl = 1/Sqrt[6.70837*10^-39];
  gsT = 106.75;
  Sup[ΛI_?NumericQ, ΓI_?NumericQ] := 
  Module[{a, ρr, Trad, tf, s, t},
     ClearAll[tf]; ClearAll[Trad]; ClearAll[s];
     tf = 10/ΓI;
     ClearAll[a]; ClearAll[ρr];
     s = NDSolve[{a'[t] == a[t]*Sqrt[(8 π)/(3 mpl^2) (ρr[t] + ΛI^4/a[t]^3 Exp[-ΓI t])], ρr'[t] + 4*Sqrt[(8 π)/(3 mpl^2) (ρr[t] + ΛI^4/a[t]^3 Exp[-ΓI t])] ρr[t] == ΓI ΛI^4/a[t]^3 Exp[-ΓI t], a[0] == 1, ρr[0] == 0}, {a, ρr}, {t, 0, tf}, MaxStepFraction -> 10^-5, MaxSteps -> 10^6];
    {a = a /. s[[1]], ρr = ρr /. s[[1]]};
    Trad[t_] := (30/(π^2 gsT) ρr[t])^(1/4);
    Trad[tf]*a[tf]]

(I then evaluate Sup[x,y] for a variety of x and y values; actually, my goal is to make a contour plot of it.)

As you can see, I'm explicitly clearing a[t] and Rho[t] before I call NDSolve; however, these don't appear to get deleted after the module evaluates. Looking at the global variables, I have a$1058, a$1186, and so on. This is causing my memory usage to spiral out of control.

The problem seems to be that DownValues are being assigned to these functions and so they're not getting cleared. (I have set $HistoryLength to zero.)

How can I get the module to really, truly clear these functions?

Answer 1

It is a little known fact and probably not well documented, but since version 9 one can use just strings as variables (dependent and independent) in NDSolve, which in this case helps to solve the memory problem in a rather elegant way:

mpl=1/Sqrt[6.70837*10^-39];
gsT=106.75;
Sup[LamdaI_?NumericQ,GammaI_?NumericQ]:=Module[{
        a,rhor,Trad,tf,s,t
    },
    tf=10/GammaI;
    s=NDSolve[{
            "a"'["t"]=="a"["t"]*Sqrt[(8 Pi)/(3 mpl^2) ("rhor"["t"]+LamdaI^4/"a"["t"]^3 Exp[-GammaI "t"])],
            "rhor"'["t"]+4*Sqrt[(8 Pi)/(3 mpl^2) ("rhor"["t"]+LamdaI^4/"a"["t"]^3 Exp[-GammaI "t"])] "rhor"["t"]==GammaI LamdaI^4/"a"["t"]^3 Exp[-GammaI "t"],"a"[0]==1,"rhor"[0]==0
        },{"a","rhor"},{"t",0,tf},
        MaxStepFraction->10^-5,MaxSteps->10^6
    ];
    a="a"/.s[[1]];
    rhor="rhor"/.s[[1]];
    Trad=Function[t,(30/(Pi^2 gsT) rhor[t])^(1/4)];
    Trad[tf]*a[tf]
]

note that I used a pure function instead of DownValues to define Trad which is also a trick which helps to prevent the known Module leakages. For me this function does not leak any variables anymore and I think should be free of memory leaks as well. Of course one could do the same thing in a programmatical way, but that needs some care to prevent leakage as for some reason even the evaluation of just the derivatives seems to be enough to make the local variables survive as Szabolcs has found out...

Note this seems to be an unwanted feature which just happend to work by accident. Newer versions (I think starting from 11.2) strings are not accepted as variable names anymore, of course making my answer useless for these newer versions...

Answer 2

From the comments that Szabolcs gave, Clear and ClearAll are ineffective, but using Remove works. So now the module reads:

    Sup[\[CapitalLambda]I_?NumericQ, \[CapitalGamma]I_?NumericQ] := 
    Module[{a, \[Rho]r, Trad, tf, s, t, result},
      tf = 10/\[CapitalGamma]I;
      s = NDSolve[{a'[t] ==a[t]*Sqrt[(8 \[Pi])/(3 mpl^2) (\[Rho]r[t] + \[CapitalLambda]I^4/a[t]^3 Exp[-\[CapitalGamma]I t])], \[Rho]r'[t] +4*Sqrt[(8 \[Pi])/(3 mpl^2) (\[Rho]r[t] + \[CapitalLambda]I^4/a[t]^3 Exp[-\[CapitalGamma]I t])] \[Rho]r[t] == \[CapitalGamma]I \[CapitalLambda]I^4/a[t]^3 Exp[-\[CapitalGamma]I t],a[0] == 1, \[Rho]r[0] == 0}, {a, \[Rho]r}, {t, 0, tf},MaxStepFraction -> 10^-5, MaxSteps -> 10^6];
     {a = a /. s[[1]], \[Rho]r = \[Rho]r /. s[[1]]};
     Trad[t_] := (30/(\[Pi]^2 gsT) \[Rho]r[t])^(1/4);
     result = Trad[tf]*a[tf];
     Remove[tf]; Remove[Trad]; Remove[s]; Remove[a]; Remove[\[Rho]r]; Remove[t];
     result
     ]

I'm still not entirely sure why Remove works when the others didn't, but it does.

Answer 3

The reason Clear did not work is because the variables are referenced from the system cache. You need to ClearSystemCache[] as well. (Or wait for the references to get cleared from the cache themselves.)

This thus works:

Sup[ΛI_?NumericQ, ΓI_?NumericQ] := 
 Module[{a, ρr, Trad, tf, s, t, result}, tf = 10/ΓI;
  s = NDSolve[{a'[t] == 
  a[t]*Sqrt[(8 π)/(3 mpl^2) (ρr[
        t] + ΛI^4/
         a[t]^3 Exp[-ΓI t])], ρr'[t] + 
   4*Sqrt[(8 π)/(3 mpl^2) (ρr[
         t] + ΛI^4/
          a[t]^3 Exp[-ΓI t])] ρr[
     t] == ΓI ΛI^4/
    a[t]^3 Exp[-ΓI t], 
 a[0] == 1, ρr[0] == 0}, {a, ρr}, {t, 0, tf}, 
MaxStepFraction -> 10^-5, MaxSteps -> 10^6];
{a = a /. s[[1]], ρr = ρr /. s[[1]]};
Trad[t_] := (30/(π^2 gsT) ρr[t])^(1/4);
result = Trad[tf]*a[tf];
Clear[tf, Trad, s, a, ρr, t];
ClearSystemCache[];
result]

Edit: In fact, one can (at least in 11.1.1) forgo the Clear entirely and just ClearSystemCache. This makes me wonder why it is necessary to do this at all, since in principle the system cache should clear old entries by itself.

Answer 4

The string as a function name method by @Albert Retey no longer work in MMA 11. But Block[] seems to work in this specific case.

  mpl = 1/Sqrt[6.70837*10^-39];
  gsT = 106.75;
  Sup[ΛI_?NumericQ, ΓI_?NumericQ] := 
  Block[{a, ρr, Trad, tf, s, t},
     tf = 10/ΓI;
     s = NDSolve[{a'[t] == a[t]*Sqrt[(8 π)/(3 mpl^2) (ρr[t] + ΛI^4/a[t]^3 Exp[-ΓI t])], ρr'[t] + 4*Sqrt[(8 π)/(3 mpl^2) (ρr[t] + ΛI^4/a[t]^3 Exp[-ΓI t])] ρr[t] == ΓI ΛI^4/a[t]^3 Exp[-ΓI t], a[0] == 1, ρr[0] == 0}, {a, ρr}, {t, 0, tf}, MaxStepFraction -> 10^-5, MaxSteps -> 10^6];
    {a = a /. s[[1]], ρr = ρr /. s[[1]]};
    Trad[t_] := (30/(π^2 gsT) ρr[t])^(1/4);
    Trad[tf]*a[tf]]

Try to run it by

Table[{Sup[10^16, 10^8], Information["a$*"]}, {3}]

You should not see any new a$XXXXX generated.