Simulations#

Here we illustrate how to run simulations, using the atomsmltr.simulation.simulator subpackage

We only provide minimal examples here. Please have a look at the examples in the atomic physic section for more examples.

Minimal example of a simulation#

Here is a very simple example, where we only include a pair of lasers along z (modelled as plane waves), and launch one ytterbium atom with an initial speed of 100 m/s along z.

We use the simulator based on Scipy’s solve_IVP method, namely ScipyIVP_3D. When we only need to perform one simulation, the easiest is to use the integrate() method.

import numpy as np
import matplotlib.pyplot as plt

from atomsmltr.environment import PlaneWaveLaserBeam
from atomsmltr.atoms import Ytterbium
from atomsmltr.simulation import Configuration, ScipyIVP_3D

# - setup atom
atom = Ytterbium()
main = atom.trans["main"] # get transition, to help setting up lasers

# - setup laser
laser_1 = PlaneWaveLaserBeam()
laser_1.direction = (0, 0, 1)
laser_1.set_power_from_I(main.Isat) # set power to reach Isat
laser_1.tag = "las1"

laser_2 = laser_1.copy() # create a copy
laser_2.direction = (0, 0, -1)  # propagating in opposite direction
laser_2.tag = "las2"

# - config
config = Configuration()
config.atom = atom
config += laser_1, laser_2
config.add_atomlight_coupling("las1", "main", -0.5 * main.Gamma)
config.add_atomlight_coupling("las2", "main", -0.5 * main.Gamma)

# - simulation
sim = ScipyIVP_3D(config=config)
t = np.linspace(0, 0.1, 1000)  # timesteps for integration
u0 = (0, 0, 0, 0, 0, 100)  # atom starts with vz=100m/s
res = sim.integrate(u0, t)

# plot
fix, axes = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
axes[0].plot(res.t * 1e3, res.y[2])
axes[0].set_ylabel("z (m)")
axes[1].plot(res.t * 1e3, res.y[5])
axes[1].set_ylabel("vz (m/s)")
for ax in axes:
    ax.set_xlabel("t (ms)")
    ax.grid()
../../_images/b5a38ffb807575ae610172ff9be927cc9844bbe0cd84968704bd29e9945d9efe.png

Parallel computation#

The simulator makes it easy to perform several simulations in parallel. For that, one has to provide a list of initial conditions u0_list and use the run() method :

# - simulation
# using the same configuration as before
sim = ScipyIVP_3D(config=config)
t = np.linspace(0, 0.1, 1000)  # timesteps for integration

# prepare initial conditions
vz_list = np.linspace(1, 130, 100)
u0_list = [(0,0,0,0,0,vz) for vz in vz_list]

# run
sim.u0_list = u0_list
res_coll = sim.run(t, npools=5, verbose=True)
Hide code cell output
100%|β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ| 100/100 [00:01<00:00, 85.50it/s]
# plot
plt.figure(figsize=(8, 4), tight_layout=True)
for res in res_coll[::2]:
    plt.plot(res.y[2], res.y[5], color="C0")
plt.xlabel("z (m)")
plt.ylabel("vz (m/z)")
plt.xlim(0, 10)
plt.ylim(0, 120)
plt.grid()
plt.show()
../../_images/2110579c219f38d4be6cc550e1abf91c8d00cd8022de4e3abb8ac042048585b9.png

Use zones to stop the simulation#

One can use zones to make the simulation stop.

To test that, let’s make a stupid simulation with atoms flying around freely

import numpy as np
import matplotlib.pyplot as plt

from atomsmltr.simulation import Configuration, ScipyIVP_3D
from atomsmltr.environment import Limits
from atomsmltr.atoms import Ytterbium

# atom
atom = Ytterbium()

# zones
xlim = Limits(-10, 10, axis=0, target="position", action="stop")
ylim = Limits(-10, 10, axis=1, target="position", action="stop")

# config
config = Configuration()
config.atom = atom
config.add_objects([xlim, ylim])

# simulator

sim = ScipyIVP_3D(config)

# set initial conditions
t = np.linspace(0, 10, 1000)
theta_list = np.linspace(0, 2 * np.pi, 50)
u0 = [(0, 0, 0, 10 * np.cos(th), 10 * np.sin(th), 0) for th in theta_list]

# run
res_coll = sim.run(t, u0, npools=5)

# compare with case without zones
config.rm_all_zones()
sim.config = config
res_coll_nozones = sim.run(t, u0, npools=5)

# plot
plt.figure(figsize=(4, 4), tight_layout=True)
for res in res_coll_nozones:
    plt.plot(res.y[0], res.y[1], color="C1")
for res in res_coll:
    plt.plot(res.y[0], res.y[1], color="C0")
plt.xlabel("x (m)")
plt.ylabel("y (m/)")
plt.xlim(-20, 20)
plt.ylim(-20, 20)
plt.grid()
plt.show()
../../_images/3b7ec8c6ae2e8ba99074b577073dea96202210c20652dbf95b6f6444aeffe717.png