I have a simple code which solves an equation by an explicit method (FTCS). It takes mathematica several minutes (mathematica 10.0.2) to finish the calculation while the same code in Fortran runs less than a second (on the same machine, no parallelization). I expected Mathematica to be slower, but that slower? I think something weird is going on. Can somebody help?
ClearAll["Global`*"];
(* define viscosity function *)
fGetNu := Compile[{r},
1.34*10^14*(r/AU2cm[1])^(3/2)
];
(* define conversion functions *)
yr2sec = Compile[{x},
3.155*10^7*x
];
AU2cm = Compile[{x},
1.496*10^13*x
];
calc = Compile[{},
(* define global constants *)
ngrid = 102;
rin = AU2cm[0.1];
rout = AU2cm[100];
tmax = yr2sec[2*10^6];
tout = yr2sec[1*10^5];
tcounter = 0;
c0 = 0.5;
(* arrays *)
x = ConstantArray[0, ngrid];
v = ConstantArray[0, ngrid];
r = ConstantArray[0, ngrid];
Σ = ConstantArray[0, ngrid];
vNew = ConstantArray[0, ngrid];
ν = ConstantArray[0, ngrid];
(* main part *)
(* setup the grid *)
dX = 2.0*(Sqrt[rout] - Sqrt[rin])/(ngrid - 2);
For[i = 1, i <= ngrid, i++,
x[[i]] = 2 Sqrt[rin] + (i - 1.5)*dX;
r[[i]] = 0.25*x[[i]]^2;
];
(* Define zero-torque (Sigma=0) boundary conditions *)
Σ[[1]] = 0; Σ[[ngrid]] = 0;
(* Define initial conditions *)
For[i = 2, i <= ngrid - 1, i++,
If[TrueQ[r[[i]] <= AU2cm[30]],
Σ[[i]] = 1000*(r[[i]]/AU2cm[1])^(-3/2);,
Σ[[i]] = 1*10^-4;
];
ν[[i]] = fGetNu[r[[i]]];
v[[i]] = 3/2 ν[[i]]*Σ[[i]]*Sqrt[r[[i]]];
];
(* initial time *)
t = 0; tnext = 0;
(* main loop *)
While[TrueQ[t <= tmax],
dtmin = 1*10^20;
For[i = 2, i <= ngrid - 1, i++,
(* this is what is in Armitage code, but NB.
he uses explicit FTCS and this breaks von Neumann condition for \
stability, which has strict inequality,
also the further from the limit, the better *)
dtzone = c0*(dX*x[[i]])^2/(24*ν[[i]]);
dtmin = Min[dtmin, dtzone];
];
dt = dtmin;
(* Evolve one timestep,store result in array Vnew *)
(* This is FTCS time step *)
For[i = 2, i <= ngrid - 1, i++,
vNew[[i]] =
v[[i]] +
12*ν[[i]]*
dt (v[[i - 1]] - 2*v[[i]] + v[[i + 1]])/(x[[i]]*dX)^2;
];
For[i = 2, i <= ngrid - 1, i++,
v[[i]] = vNew[[i]];
Σ[[i]] = 2/3 v[[i]]/(ν[[i]]*Sqrt[r[[i]]]);
];
t = t + dt;
If[TrueQ[t >= tnext],
Print["Writing a slice at t = ", t/yr2sec[1]];
(*For[i=2,i≤ngrid-1,i++,
Print[{t,
r〚i〛,Σ\
〚i〛,ν〚i\
〛}];
];*)
tcounter = tcounter + 1;
tnext = tnext + tout;
];
];
Print["Wrote ", ngrid - 2, " radial results per time step."];
Print["Wrote ", tcounter - 1, " time steps."];
];
gives
In[103]:= calc[] // AbsoluteTiming
During evaluation of In[103]:= Writing a slice at t = 16.6278
During evaluation of In[103]:= Writing a slice at t = 100016.
During evaluation of In[103]:= Writing a slice at t = 200016.
During evaluation of In[103]:= Writing a slice at t = 300015.
During evaluation of In[103]:= Writing a slice at t = 400015.
During evaluation of In[103]:= Writing a slice at t = 500015.
During evaluation of In[103]:= Writing a slice at t = 600014.
During evaluation of In[103]:= Writing a slice at t = 700014.
During evaluation of In[103]:= Writing a slice at t = 800013.
During evaluation of In[103]:= Writing a slice at t = 900013.
During evaluation of In[103]:= Writing a slice at t = 1.00001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.10001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.20001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.30001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.40001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.50001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.60001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.70001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.80001*10^6
During evaluation of In[103]:= Writing a slice at t = 1.90001*10^6
During evaluation of In[103]:= Writing a slice at t = 2.00001*10^6
During evaluation of In[103]:= Wrote 100 radial results per time step.
During evaluation of In[103]:= Wrote 20 time steps.
Out[103]= {233.787534, Null}
@Jens has provided awesome answer, thank you so much. Documentation is quite scarce in this regard. I have three questions though:
- You have the following twice in: vNew = Table[0., {ngrid}]; [Nu] = Table[0., {ngrid}]; But if I delete one occurrence the code slows down (with the duplicity the code runs: 0.769044, without: 0.782045). I understand that the time difference is negligible, yet it exists and it confuses me why.
- I do not really understand the interpretation of input by Mathematica, is 2 an integer by default? Because if so, I can understand the need to use 2. as you have in the code. On the other hand in some places, you do not use 2.0 but only 2.
- What is the advantage of using Table instead of ConstantArray at the beginning?
- Next is a more general question, in one of the comments above, @belisarius mentioned to avoid procedural loops, so I read some text about it where Mathematica authors suggest to use Map etc. But I fail to understand how that helps; in essence they are also loops, even when they can be vectorized or parallelized I doubt they would be able to make 2 orders of magnitude of processing time speed up.
When I enable the line with Print["Writing a slice at t = ", t]; I get uppon run CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. So I dug a bit deeper and right after the line where I initialize the t variable I put:
t = 0.0; Print[t];
And now,
In[33]:= AbsoluteTiming[calc[]]
During evaluation of In[33]:= t
During evaluation of In[33]:= CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. >>
During evaluation of In[33]:= 0.
During evaluation of In[33]:= Writing a slice at t = 5.24607*10^8
Out[33]= $Aborted
My conclusion is that for whatever reason there is some sort of "dry run" when Mathematica refuses to do what it is told - set t to 0.0 (and this is not the first time, I came across this strange behaviour with LogPlot), only after that the code runs as it should. This would also explain the necessity to use TrueQ - see below.
@shrx the TrueQ was necessary because Mathematica refused to compile if I had While[t<=tmax,...], it kept complaining. I do not understand, though, how Mathematica can tell during compilation time whether the expression is or is not evaluated.
fGetNugets re-compiled every single time it's invoked because it's defined withSetDelayed. Probably if you state the problem you want to solve, someone can come up with a better method using built-in functionality. $\endgroup$