How can I solve a three-dimensional Gross-Pitaevskii equation?

Question

I want to solve a three-dimensional Gross-Pitaevskii equation with an uniform potential. The related papers are Emergence of a Turbulent Cascade in a Quantum Gas and Synthetic dissipation and cascade fluxes in a turbulent quantum gas.

The code I tried was

hbar = 1;

kB = 1.380649*10^-23;
\[Mu] = 1;(*2 10^-9 kB;*)
\[Xi] = 1;(*1.2 10^-6;*)

Udepth = 30 \[Mu];
Uz = 27/1.2 \[Xi];
UR = 16/1.2 \[Xi];

N0 = 1. 10^5;
n0 = N0/(Pi UR^2 Uz);
g = \[Mu]/n0;
\[CapitalDelta]U = \[Mu]; (*shaking amplitude*)
\[Omega]shaking =
 2 Pi 8;

Lx = 40 \[Xi]; Ly = 40 \[Xi]; Lz = 40 \[Xi];
x0 = Lx;
M = 1;
Vtrap[x_, y_, z_, t_] :=
 If[x^2 + y^2 < UR && -Uz/2 < z < Uz/2, 0,
   Udepth] + \[CapitalDelta]U Sin[\[Omega]shaking t] z/Uz
eqs = { D[\[Psi]1[t, x, y, z],
     t] == ((-g n0 Laplacian[\[Psi]1[t, x, y, z], {x, y,
         z}]) + (Vtrap[x, y, z, t] +
        g (Abs[\[Psi]1[t, x, y, z]]^2)) \[Psi]1[t, x, y, z])/(
    I g n0 hbar)};
bc = {\[Psi]1[0, x, y, z] ==
    Exp[-x^2/(10 \[Xi])^2] Exp[-y^2/(10 \[Xi])^2] Exp[-z^2/(10 \
\[Xi])^2](*1/Sqrt[rcx rcy Pi]Exp[-x^2/rcx^2/2]Exp[-y^2/rcy^2/2]*),
   DirichletCondition[\[Psi]1[t, x, y, z] == 0, True]};
solv = NDSolve[{eqs, bc}, \[Psi]1, {t, 0, 0.5}, {x, -Lx/2,
   Lx/2}, {y, -Ly/2, Ly/2}, {z, -Lz/2, Lz/2},
    Method -> {"PDEDiscretization" -> {"MethodOfLines", \
{"SpatialDiscretization" -> "FiniteElement"}}}]

But this costs too much time and can not output. What should I do?

Answer 1

After proper normalization we have

hbar = 1;

kB = 1.380649*10^-23;
\[Mu] = 1;(*2 10^-9 kB;*)
\[Xi] = 1;(*1.2 10^-6;*)

Udepth = 30 \[Mu];
Uz = 27/1.2 \[Xi];
UR = 16/1.2 \[Xi];

N0 = 1. 10^5;
n0 = N0/(Pi UR^2 Uz);
g = \[Mu]/n0;
\[CapitalDelta]U = \[Mu]; (*shaking amplitude*)
\[Omega]shaking = 2 Pi 8;

Lx = 40 \[Xi]; Ly = 40 \[Xi]; Lz = 40 \[Xi];
x0 = Lx; g1 = g/(g n0 );
M = 1;
Vtrap[x_, y_, z_, t_] := 
 If[x^2 + y^2 < UR && -Uz/2 < z < Uz/2, 0, 
   Udepth] + \[CapitalDelta]U Sin[\[Omega]shaking t] z/(Uz g n0)
eqs = {I hbar D[\[Psi]1[t, x, y, z], 
      t] == ((-Laplacian[\[Psi]1[t, x, y, z], {x, y, z}]) + (Vtrap[x, 
          y, z, t] + g1 (Abs[\[Psi]1[t, x, y, z]]^2)) \[Psi]1[t, x, y,
         z])};
bc = {\[Psi]1[0, x, y, z] == 
    Exp[-x^2/(10 \[Xi])^2] Exp[-y^2/(10 \[Xi])^2] Exp[-z^2/(10 \
\[Xi])^2](*1/Sqrt[rcx rcy Pi]Exp[-x^2/rcx^2/2]Exp[-y^2/rcy^2/2]*), 
   DirichletCondition[\[Psi]1[t, x, y, z] == 0, True]};

solv = NDSolve[{eqs, bc}, \[Psi]1, {t, 0, 0.5}, {x, -Lx/2, 
   Lx/2}, {y, -Ly/2, Ly/2}, {z, -Lz/2, Lz/2}]

Visualization

Table[Plot3D[
  Evaluate[Abs[\[Psi]1[t, 0, y, z] /. solv[[1]]]], {y, -Ly/2, 
   Ly/2}, {z, -Lz/2, Lz/2}, PlotLabel -> t, 
  ColorFunction -> "ArmyColors", Mesh -> None, PlotRange -> All], {t, 
  0, .5, .1}]