Adds convergence criteria for mayavi 3D, if timesteps are large enough

PKGBUILD

@@ -1,7 +1,7 @@
 # Maintainer: Daniel Scheiermann  <[email protected]>
 _name=supersolids
 pkgname=python-${_name}
-pkgver=0.1.13
+pkgver=0.1.14
 pkgrel=1
 pkgdesc="Notes and script to supersolids"
 url="https://github.com/Scheiermann/${_name}"

setup.py

@@ -8,7 +8,7 @@
 
 setup(
     name="supersolids",
-    version="0.1.13",
+    version="0.1.14",
     packages=["", "supersolids"],
     package_data={"supersolids": ["results/split_time_imag.mp4",
                                   "results/split_time_real.mp4",

supersolids/MayaviAnimation.py

@@ -48,21 +48,28 @@ def get_image_path(dir_path: Path, dir_name: str = "movie", counting_format: str
 
 
 @mlab.animate(delay=10, ui=True)
-def animate(System: Schroedinger.Schroedinger, x_lim=(-1, 1), y_lim=(-1, 1), z_lim=(-1, 1)):
+def animate(System: Schroedinger.Schroedinger, accuracy=10**-6, x_lim=(-1, 1), y_lim=(-1, 1), z_lim=(-1, 1),
+            slice_x_index=0, slice_y_index=0):
     prob_3d = np.abs(System.psi_val) ** 2
     p = mlab.contour3d(System.x_mesh, System.y_mesh, System.z_mesh, prob_3d,
                        colormap="spectral", opacity=0.5, transparent=True)
 
     slice_x = mlab.volume_slice(System.x_mesh, System.y_mesh, System.z_mesh, prob_3d, colormap="spectral",
                                 plane_orientation="x_axes",
-                                slice_index=System.resolution // 2,
+                                slice_index=slice_x_index,
                                 extent=[*x_lim, *y_lim, *z_lim])
     slice_y = mlab.volume_slice(System.x_mesh, System.y_mesh, System.z_mesh, prob_3d, colormap="spectral",
                                 plane_orientation="y_axes",
-                                slice_index=System.resolution // 2,
+                                slice_index=slice_y_index,
                                 extent=[*x_lim, *y_lim, *z_lim])
     for i in range(0, System.timesteps):
+        s_old = System.s
         System.time_step()
+        s_rel = np.abs((System.s - s_old) / System.s)
+        print(f"s_rel: {s_rel}")
+        if s_rel < accuracy:
+            print(f"accuracy reached: {s_rel}")
+            break
         prob_3d = np.abs(System.psi_val) ** 2
         slice_x.mlab_source.trait_set(scalars=prob_3d)
         slice_y.mlab_source.trait_set(scalars=prob_3d)
@@ -114,7 +121,7 @@ def create_movie(self, input_data_file_pattern="*.png", filename="anim.mp4"):
 # Script runs, if script is run as main script (called by python *.py)
 if __name__ == "__main__":
     Harmonic = Schroedinger.Schroedinger(resolution=2 ** 6, timesteps=100, L=3, dt=1.0, g=1.0, imag_time=True, dim=3,
-                                         s=1,
+                                         s=1.1, E=1.0,
                                          psi_0=functions.psi_gauss_3d,
                                          V=functions.v_harmonic_3d,
                                          psi_sol=functions.thomas_fermi

supersolids/Schroedinger.py

@@ -30,10 +30,11 @@ class Schroedinger(object):
     WARNING: We don't use Baker-Campell-Hausdorff formula, hence the accuracy is small. This is just a draft.
     """
 
-    def __init__(self, resolution, timesteps, L, dt, g=0.0, imag_time=False, dim=1, s=1.0,
+    def __init__(self, resolution, timesteps, L, dt, g=0.0, imag_time=False, dim=1, s=1.1, E=1.0,
                  psi_0=functions.psi_gauss_1d,
                  V=functions.v_harmonic_1d,
-                 psi_sol=functions.thomas_fermi
+                 psi_sol=functions.thomas_fermi,
+                 mu_sol=functions.mu_3d,
                  ):
         """
         Parameters
@@ -49,8 +50,16 @@ def __init__(self, resolution, timesteps, L, dt, g=0.0, imag_time=False, dim=1,
         self.g = float(g)
         self.imag_time = imag_time
         self.dim = dim
+
+        print(f"{self.g}")
+        self.mu_sol = mu_sol(self.g)
+
         # s = - ln(N) / (2 * dtau), where N is the norm of the psi
         self.s = s
+        print(f"init s: {self.s}")
+
+        # E = mu - 0.5 * g * int psi_val ** 2
+        self.E = E
 
         self.psi = psi_0
         self.V = V
@@ -140,13 +149,13 @@ def __init__(self, resolution, timesteps, L, dt, g=0.0, imag_time=False, dim=1,
 
         self.V_z_line = None
 
-    def get_norm(self):
+    def get_norm(self, p=2.0):
         if self.dim == 1:
-            psi_norm = np.sum(np.abs(self.psi_val) ** 2.0) * self.dx
+            psi_norm = np.sum(np.abs(self.psi_val) ** p) * self.dx
         elif self.dim == 2:
-            psi_norm = np.sum(np.abs(self.psi_val) ** 2.0) * self.dx * self.dy
+            psi_norm = np.sum(np.abs(self.psi_val) ** p) * self.dx * self.dy
         elif self.dim == 3:
-            psi_norm = np.sum(np.abs(self.psi_val) ** 2.0) * self.dx * self.dy * self.dz
+            psi_norm = np.sum(np.abs(self.psi_val) ** p) * self.dx * self.dy * self.dz
         else:
             print("Spatial dimension over 3. This is not implemented.", file=sys.stderr)
             sys.exit(1)
@@ -178,8 +187,12 @@ def time_step(self):
 
         # for self.imag_time=False, renormalization should be preserved, but we play safe here (regardless of speedup)
         # if self.imag_time:
-        psi_norm_after_evolution = self.get_norm()
+        psi_norm_after_evolution = self.get_norm(p=2.0)
+        psi_quadratic_integral = self.get_norm(p=4.0)
 
-        self.psi_val /= np.sqrt(psi_norm_after_evolution)
+        self.s = - np.log(psi_norm_after_evolution) / (2.0 * self.dt)
+        self.E = self.s - 0.5 * self.g * psi_quadratic_integral
 
-        # self.s = - np.log(self.get_norm()) / (2.0 * self.dt)
+        self.psi_val /= np.sqrt(psi_norm_after_evolution)
+        # print(f"s: {self.s}")
+        print(f"E: {self.E}, E_sol: {self.mu_sol - 0.5 * self.g * psi_quadratic_integral}")

supersolids/functions.py

@@ -221,6 +221,13 @@ def thomas_fermi(x, g):
         return None
 
 
+def mu_3d(g):
+    # mu is the chemical potential
+    mu = ((15 * g) / (16 * np.sqrt(2) * np.pi)) ** (2 / 5)
+
+    return mu
+
+
 def v_harmonic_1d(x):
     return 0.5 * x ** 2
 

supersolids/main.py

@@ -24,22 +24,26 @@
 from supersolids import Schroedinger
 
 
-def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, dim=1, s=1,
+def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, dim=1, s=1.1, E=1.0, accuracy=10**-6,
                   psi_0=functions.psi_gauss_1d,
                   V=functions.v_harmonic_1d,
                   psi_sol=functions.thomas_fermi,
+                  mu_sol=functions.mu_3d,
                   file_name="split.mp4",
                   x_lim=(-1, 1),
                   y_lim=(-1, 1),
                   z_lim=(0, 1.0),
+                  slice_x_index=0, slice_y_index=0,
                   view_height=20.0,
                   view_angle=45.0,
                   view_distance=10.0
                   ):
     with run_time.run_time():
-        Harmonic = Schroedinger.Schroedinger(resolution, timesteps, L, dt, g=g, imag_time=imag_time, dim=dim, s=s,
+        Harmonic = Schroedinger.Schroedinger(resolution, timesteps, L, dt, g=g, imag_time=imag_time, dim=dim,
+                                             s=s, E=E,
                                              psi_0=psi_0, V=V,
                                              psi_sol=psi_sol,
+                                             mu_sol=mu_sol,
                                              )
 
     if dim < 3:
@@ -64,8 +68,13 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, dim=1, s=1,
         may = MayaviAnimation.MayaviAnimation(dim=dim)
         with run_time.run_time():
             # may.animate(Harmonic)
-            may.animate(Harmonic, x_lim=x_lim, y_lim=y_lim, z_lim=z_lim)
+            may.animate(Harmonic, accuracy=accuracy, x_lim=x_lim, y_lim=y_lim, z_lim=z_lim,
+                        slice_x_index=slice_x_index, slice_y_index=slice_y_index)
         mlab.show()
+        # TODO: close window after last frame
+        # print(f"{Harmonic.t}, {Harmonic.dt * Harmonic.timesteps}")
+        # if Harmonic.t >= Harmonic.dt * Harmonic.timesteps:
+        #     mlab.close()
         may.create_movie(input_data_file_pattern="*.png", filename=file_name)
 
 
@@ -81,15 +90,16 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, dim=1, s=1,
     # constants needed for the Schroedinger equation
     g = 10.0
     g_step = 10
-    dt = 1.0
+    dt = 0.1
 
     # box length [-L,L]
     # generators for L, g, dt to compute for different parameters
-    L_generator = (10,)
-    G = (i for i in np.arange(g, g + g_step, g_step))
+    L_generator = (4,)
+    g_generator = (i for i in np.arange(g, g + g_step, g_step))
+    mu_sol_list = [functions.mu_3d for g in g_generator]
     factors = np.linspace(0.2, 0.3, max_workers)
-    DT = (i * dt for i in factors)
-    cases = itertools.product(L_generator, G, DT)
+    dt_generator = (i * dt for i in factors)
+    cases = itertools.product(L_generator, g_generator, dt_generator, mu_sol_list)
 
     # functions needed for the Schroedinger equation (e.g. potential: V, initial wave function: psi_0)
     V_1d = functions.v_harmonic_1d
@@ -102,27 +112,38 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, dim=1, s=1,
     psi_0_2d = functools.partial(functions.psi_gauss_2d_pdf, mu=[0.0, 0.0], var=np.array([[1.0, 0.0], [0.0, 1.0]]))
     psi_0_3d = functools.partial(functions.psi_gauss_3d, a=1, x_0=0, y_0=0, z_0=0, k_0=0)
 
+    # TODO: get mayavi lim to work
     # 3D works in single core mode
-    simulate_case(resolution, timesteps=30, L=L_generator[0], g=g, dt=dt, imag_time=True, dim=3, s=1,
-                  psi_0=psi_0_3d, V=V_3d, psi_sol=functions.thomas_fermi, file_name="anim.mp4",
-                  x_lim=(-8, 8), y_lim=(-5, 5), z_lim=(0, 0.4),
+    simulate_case(resolution, timesteps=500, L=L_generator[0], g=g, dt=dt, imag_time=True, dim=3,
+                  s=1.1, E=1.0,
+                  psi_0=psi_0_3d, V=V_3d,
+                  psi_sol=functions.thomas_fermi, mu_sol=mu_sol_list[0],
+                  accuracy=10 ** -6,
+                  file_name="anim.mp4",
+                  x_lim=(-2, 2), y_lim=(-2, 2), z_lim=(0, 0.4),
+                  slice_x_index=resolution//3, slice_y_index=resolution//3,
                   view_height=15.0, view_angle=75.0, view_distance=10.0
                   )
     print("Single core done")
 
+    # TODO: get mayavi concurrent to work
     i: int = 0
-    with futures.ProcessPoolExecutor(max_workers=max_workers) as e:
-        for L, g, dt in cases:
+    with futures.ProcessPoolExecutor(max_workers=max_workers) as executor:
+        for L, g, dt, mu_sol in cases:
             i = i + 1
-            print(f"i={i}, L={L}, g={g}, dt={dt}")
+            print(f"i={i}, L={L}, g={g}, dt={dt}, mu={mu_sol}")
             file_name = f"split_{i:03}.mp4"
             psi_sol = functools.partial(functions.thomas_fermi, g=g)
-            e.submit(simulate_case, resolution, timesteps=30, L=L, g=g, dt=dt, imag_time=True, dim=3, s=1,
-                     psi_0=psi_0_3d, V=V_3d, psi_sol=psi_sol, file_name=file_name,
-                     x_lim=(-8, 8),
-                     y_lim=(-5, 5),
-                     z_lim=(0, 0.4),
-                     view_height=15.0,
-                     view_angle=75.0,
-                     view_distance=10.0
-                     )
+            executor.submit(simulate_case, resolution, timesteps=30, L=L, g=g, dt=dt, imag_time=True,
+                            dim=3, s=1.1, accuracy=10**-6,
+                            psi_0=psi_0_3d, V=V_3d,
+                            psi_sol=psi_sol, mu_sol=mu_sol,
+                            file_name=file_name,
+                            x_lim=(-8, 8),
+                            y_lim=(-5, 5),
+                            z_lim=(0, 0.4),
+                            slice_x_index=0, slice_y_index=0,
+                            view_height=15.0,
+                            view_angle=75.0,
+                            view_distance=10.0
+                            )