# Speeding up solving nonlinear Schroedinger equation in 3D with NDSolve with periodic boundary conditions

I have a question on speeding up solving nonlinear Schroedinger equation in 3D with NDSolve with periodic boundary conditions.

I build the ODE system with NDSolveFiniteDifferenceDerivative[] and NDSolveProcessEquations[]. After that I iterate with NDSolveProcessSolutions[InState, "Forward"] .

For 20x20x20 grid ProcessEquations works for about 15 minutes. After that the iterations take from 20 to 80 seconds. I do not understand why just building the system of ODE takes so long. For 50x50x50 grid NDSolveProcessEquations could not finish overnight.

Please tell me if I am doing somthing wrong or if there are ways to speed up.

My code is given below.

np = 10;
nb = 5;
wp = 15;
n = 20;
npoints = 5000;

xgrid = ygrid = zgrid = Range[1, n + 1];

sd2x = NDSolveFiniteDifferenceDerivative[{2, 0, 0},{xgrid,ygrid,zgrid},PeriodicInterpolation -> {True, True,True},"DifferenceOrder" -> 2];
sd2y = NDSolveFiniteDifferenceDerivative[{0, 2, 0},{xgrid,ygrid,zgrid}, PeriodicInterpolation -> {True, True, True},"DifferenceOrder" -> 2];
sd2z = NDSolveFiniteDifferenceDerivative[{0, 0, 2},{xgrid,ygrid,zgrid}, PeriodicInterpolation -> {True, True, True},"DifferenceOrder" -> 2];
termlap = Drop[sd2x[Outer[f[t, #1, #2, #3] &,xgrid, ygrid,zgrid]],{-1},{-1},{-1}] + Drop[sd2y[Outer[f[t, #1, #2, #3] &,xgrid,ygrid,zgrid]],{-1},{-1}, {-1}] + Drop[sd2z[Outer[f[t, #1, #2, #3] &,xgrid,ygrid,zgrid]],{-1},{-1}, {-1}];
term1 = Drop[Outer[(f[t, #1, #2, #3] Abs[f[t, #1, #2, #3]]^2) &, xgrid,ygrid,zgrid], {-1}, {-1}, {-1}];
term2 = Drop[Outer[(f[t, #1, #2, #3] Abs[f[t, #1, #2, #3]]^4) &, xgrid,ygrid,zgrid], {-1}, {-1}, {-1}];
termf = Drop[Outer[(f[t, #1, #2, #3]) &, xgrid, ygrid, zgrid], {-1}, {-1}, {-1}];
termf1 = Flatten@termf;


The initial conditions are more complicated but they do not really matter.

incd = Thread[(termf1 /. t -> 0) == Flatten[1 + RandomReal[.1, {n, n, n}]]];
eqns = Thread[Flatten[I \!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]termf\) + 0.25 termlap + term1 - 0.1875 term2] == 0];
{InState} =  NDSolveProcessEquations[Join[eqns, incd], termf1, {t}]


The last line above is the the longest to evaluate in the program.

Next I iterate with

Ani3D[\[CapitalDelta]T_, frames_] := Block[{OutState, time1, time2, mem1, mem2},
time1 = TimeUsed[]; mem1 = MemoryInUse[];
OutState = NDSolveProcessSolutions[InState, "Forward"];
Print[Plotic3[OutState, 0.]];
Print[{0, TimeUsed[] - time1, MemoryInUse[] - mem1}];
time2 = TimeUsed[]; mem2 = MemoryInUse[];
Do[
NDSolveIterate[InState, j \[CapitalDelta]T];
OutState = NDSolveProcessSolutions[InState, "Forward"];
Print[Plotic3[OutState, j \[CapitalDelta]T]];
time1 = TimeUsed[]; mem1 = MemoryInUse[];
Print[{j \[CapitalDelta]T, time1 - time2, mem1 - mem2}];
time2 = time1; mem2 = mem1;
, {j, frames}]]


The function plotic3 visualizes the output plotting "npoints" random points with the probability density given by the absolute value of the solution squared and calculates the correlation function using Corr[], and plots it.

Plotic3[fu_, tt_] := Block[{dots, height, length, fabs, ffi, npi, pic, pic1, li},
fabs = Abs[termf /. t -> tt /. fu]^2;
height = Max[fabs];
ffi = ListInterpolation[fabs, InterpolationOrder -> 0];
npi = Round[npoints Max[fabs]];
dots = {Table[Random[Real, {1, n}], {npi}],
Table[Random[Real, {1, n}], {npi}],
Table[Random[Real, {1, n}], {npi}],
Table[Random[Real, {0, height}], {npi}]}\[Transpose];
pic = Graphics3D[Point[Transpose[Select[dots, #[[4]] < ffi[#[[1]], #[[2]], #[[3]]] &]][[{1, 2, 3}]]\[Transpose]]];
li = Reap[Corr[fu, tt]][[2]][[1]];
length = N[Sqrt[3] n/2 1/nb];
pic1 = ListPlot[
Reap[Do[Sow[{length (i - 1/2), Mean[Select[li, (#[[1]] < i length) && (#[[1]] >= (i-1) length) &].{0,1}]}], {i, 1, nb, 1}]][[2]], Filling -> Axis];
Sow[Show[GraphicsRow[{pic, pic1}]]]
]


Corr starts from the point with the coordinates (1,1,1), chooses another point randomly and calculates the distance between the points and the product of the abs. values of the solutions in these points. It repeats this task for (number of bins) * (number of points in the bin rougthly) times. And the distance between the points is physical when it is smaller than the half of the box diagonal.

Corr[fj_, tt_] := Block[{len, x1, x2, y1, y2, z1, z2},
x1 = y1 = z1 = 1;
Do[len = Round[RandomReal[{0, Sqrt[3] n/2}] {Sin[#1]Cos[#2],Sin[#1]Sin[#2],Cos[#2]}&[RandomReal[2 \[Pi]],RandomReal[{-\[Pi]/2, \[Pi]/2}]]];
x2 = Mod[x1 + len[[1]], n] /. 0 -> n;
y2 = Mod[y1 + len[[2]], n] /. 0 -> n;
z2 = Mod[z1 + len[[3]], n] /. 0 -> n;
Sow[{N[Sqrt[Min[Abs[x1 - x2], n - Abs[x1 - x2]]^2 + Min[Abs[y1 - y2], n - Abs[y1 - y2]]^2 + Min[Abs[z1 - z2], n - Abs[z1 - z2]]^2]], N[Abs[f[t, x1, y1, z1] f[t, x2, y2, z2] /. t -> tt /. fj]^2]}];
x1 = x2;
y1 = y2;
z1 = z2;
, {i, 0, np nb, 1}]];


Finally I get the pictures with e.g.

res = Reap[Ani3D[10., 5]]
`

In fact the main bottleneck is the construction of the InState, the visualization as it is here is just to see that it gives a reasonable solution. I will need a much larger grid for the real analysis. When I talk to people doing Mathcad, they say that it takes them seconds to solve such equations on such a small grid. And I do not understand the problem.

Thank you.