# Improving NDSolve computation speed for 3D stationary PDE

Continuing my last question, I am trying to solve 3D Schrodinger equation (with mass equals to 1) and with no potential (I will add a potential later). The analytical solution for this equation is simply the spherical harmonics times the spherical Bessel function.

As in my last question, I defined:

rCoor[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2]
θCoor[x_, y_, z_] := ArcTan[z, Sqrt[x^2 + y^2]]
ϕCoor[x_, y_, z_] := ArcTan[x, y]
YlmCart[l_, m_, x_, y_, z_] := SphericalHarmonicY[l, m, θCoor[x, y, z], ϕCoor[x, y, z]]

Now I'm trying to solve the equation, with l = m = 0 (the spherical harmonic paramters) and for energy parameter of Ef = 1:

Needs["NDSolveFEM"];
l = 0;
m = 0;
Ef = 1;
kf = Sqrt[2 Ef];
rmax = 15 kf^-1;
mesh = ToElementMesh[Ball[{0, 0, 0}, rmax], MaxCellMeasure -> {"Length" -> 1.2}];
boundary[l_, m_, k_, x_, y_, z_] :=
SphericalBesselJ[l, k rCoor[x, y, z]] YlmCart[l, m, x, y, z]
op = -(1/2) Laplacian[u[x, y, z], {x, y, z}] - Ef u[x, y, z];
sol = NDSolveValue[{op == 0, DirichletCondition[u[x, y, z] ==
boundary[l, m, kf, x, y, z] , True]}, u, {x, y, z} ∈ mesh];

This is working just fine for small range, but I'm trying to solve it for much larger range, like rmax = 1000 kf^-1 or even more. After rmax = 30 kf^-1 the computation becomes very very slow and does not finish.

I broke the NDSolve to parts and found out that the longest part is in the LinearSolve function, which for rmax = 30 kf^-1 has to solve about 180,000 equations.

What can I do to speed things up?

• I'm pretty sure the size of the LinearSolve call is directly dependent on the size of the mesh. You could increase the "Length" of your mesh, but you have to balance that against the loss of accuracy. I expect LinearSolve does its work as fast as is possible. Perhaps someone else will know more. – Michael E2 Aug 16 '15 at 12:18
• Yes, it is directly dependent on the Length of the mesh. The question is how to overcome this obstacle. Is there a way to change the mesh resolution and still get accurate results? Is there a way to speed up LinearSolve? Maybe change some settings? – Amit Abir Aug 16 '15 at 12:23
• Have a look at the section Solving Memory-Intensive PDEs in the documentation. – user21 Aug 16 '15 at 14:54
• @AmitAbir Why not solve the equation in spherical coordinates rather than Cartesian coordinates? Also, if the number of modes in theta and phi is small, why not expand the differential equation in spherical harmonics? – bbgodfrey Aug 17 '15 at 2:38
• @bbgodfrey When I'm using spherical coordinates, I get the following message: The PDE coefficient {{-(1/2),0,0},{0,-(1/(2 r^2)),0},{0,0,-(Csc[\[Theta]]^2/(2 r^2))}} does not evaluate to a numeric matrix of dimensions {3,3}. I thought that FEM does not work with non constant coefficients. – Amit Abir Aug 17 '15 at 19:04

For the specific parameters in the question,

rmax = 15 kf^-1;
mesh = ToElementMesh[Ball[{0, 0, 0}, rmax], MaxCellMeasure -> {"Length" -> 1.2}]
(* ElementMesh[ ..., {TetrahedronElement[<10795>]}] *)
MaxMemoryUsed[sol = NDSolveValue[{op == 0,
DirichletCondition[u[x, y, z] == boundary[l, m, kf, x, y, z] , True]},
u, {x, y, z} ∈ mesh]] // AbsoluteTiming
(*  {0.974145, 315694792} *)

In other words, for these parameters NDSolve uses a mesh to almost 10000 tetrahedrons and 300MB of memory. Computational time is negligible on my four-processor PC. (Typically, these computations use about 60% of my total CPU capacity.) Both memory and run-time can be reduced by using the NDSolve option

Method -> {"PDEDiscretization" -> {"FiniteElement", "LinearSolveMethod" -> {"Pardiso"}}

as described in Solving Memory-Intensive PDEs. (My thanks to User21 for suggesting this reference.) The new timing and memory usage are

(* {0.780953,161811160} *)

Thus, the answer to the specific question posed is to use the "Pardiso" LinearSolve option to reduce memory requirements by between half and two-thirds.

Of course, using this option only delays the point at which the problem becomes too large for a PC. For instance, I run out of memory on an 8GB PC just beyond

rmax = 30 kf^-1;
mesh = ToElementMesh[Ball[{0, 0, 0}, rmax], MaxCellMeasure -> {"Length" -> 1.1}]
(* ElementMesh[ ..., {TetrahedronElement[<121663>]}] *)
MaxMemoryUsed[sol = NDSolveValue[{op == 0,
DirichletCondition[u[x, y, z] == boundary[l, m, kf, x, y, z] , True]},
u, {x, y, z} ∈ meshMethod -> {"PDEDiscretization" -> {"FiniteElement",
"LinearSolveMethod" -> {"Pardiso"}}}]] // AbsoluteTiming
(*  {91.0104, 4829279056} *)

To save more memory, one could design a custom mesh with high resolution in the radial direction and low resolution in the angular directions, appropriate for small l and m. Unfortunately, the problem described by the OP in comments above may require high resolution in all directions.