gpe.bec_basic ============= .. py:module:: gpe.bec_basic .. autoapi-nested-parse:: Code for simulating a single component BEC. This is the simplest code in the project, implementing all of the basic features without excess complications. It is somewhat excessively commented to provide a pedagogical introduction for users unfamiliar with the techniques and using python for scientific computing. A more complete implementation can be found in bec.py, but this relies on external packages to implement, for example, different bases, making it less suitable for teaching. Classes ------- .. autoapisummary:: gpe.bec_basic.Units gpe.bec_basic.State Module Contents --------------- .. py:class:: Units Units and physical constants. This class is simply a container for physical constants and units. We have chosen values here that are relevant for 87Rb which is studied at WSU in Peter Engels' lab. .. py:attribute:: hbar :value: 1.0 .. py:attribute:: micron :value: 1.0 .. py:attribute:: mm :value: 1000.0 .. py:attribute:: cm :value: 10000.0 .. py:attribute:: nm :value: 0.001 .. py:attribute:: meter :value: 1000000.0 .. py:attribute:: u :value: 1.0 .. py:attribute:: kg :value: 6.0221408585491615e+26 .. py:attribute:: G :value: 1.0 .. py:attribute:: mG :value: 0.001 .. py:attribute:: m :value: 86.909187 .. py:attribute:: a_B :value: 5.291772105498088e-05 .. py:attribute:: Hz :value: 1.5746097513521148e-05 .. py:attribute:: kHz :value: 0.015746097513521146 .. py:attribute:: MHz :value: 15.746097513521146 .. py:attribute:: s :value: 63507.79925891489 .. py:attribute:: ms :value: 63.50779925891489 .. py:attribute:: nK :value: 0.002061483743918704 .. py:attribute:: scattering_lengths .. py:attribute:: a .. py:attribute:: mu_B .. py:class:: State(Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), ws=np.array([np.sqrt(8) * 126.0, 126.0, 126.0]) * u.Hz, g=4 * np.pi * u.hbar**2 * u.a / u.m, m=u.m, mu=None, hbar=u.hbar, x_TF=None, cooling_phase=1.0, symmetric_grid=False, t=0) Bases: :py:obj:`gpe.mixins.StateMixin`, :py:obj:`pytimeode.mixins.ArrayStateMixin` State class implementing the GPE for a single-component BEC. To provide time-evolution functionality, the `pytimeode` package requires that you provide a `State` class that implements the required behaviour of at least one of the `IStateFor*Evolvers` interfaces. Most of the required behaviour can be provided by the `mixins.ArrayStateMixin` class which we inherit from as long as we provide an array called `data` which represents the wavefunction. On top of this, we need to provide the following three methods used by the `pytimeode.evolvers`: * `compute_dy_dt(dy)`: Required by interface IStateForABMEvolvers This computes the rhs of the time-dependent differential equation and is all that is needed by the `pytimeode.evolvers.EvolverABM` evolver. This is a high-order generic evolver. It requires a rather small step size to ensure that it behaves well, but once the step size is small enough, it is converges to high accuracy. For example: there is often a critical step size above which the evolution will blow up. If you reduce this by a factor of 2, then you will typically obtain accuracy to 5 digits or so, at least for modest evolution lengths. The ABM evolver requires about 10 copies of the state to start up, and 8 copies of the state to be maintained for evolution. This can be a memory issue if the states are very large. * `apply_exp_K(dt)`, `apply_exp_V(dt, state)`: Required by interface `IStateForSplitEvolvers` These two methods allow you to use `pytimeode.evolvers.EvolverSplit` based on the Trotter decomposition. This is rather specialized to systems like the GPE where the exponential of the kinetic and potential terms can be computed exactly. The ODE solver is not as high order as the ABM method, but is manifestly unitary. Thus, you can often get away with very large step sizes and the system will still behave "reasonably" meaning that the evolution will not be accurate quantitatively, but will be qualitatively reasonable so you can get an idea of what is happening. Unfortunately, to achieve quantitative accuracy, you must generally go to very small step sizes. Another advantage of the split operator methods is that they require use only 2 or 3 copies of the state, and so are better for large states that might be a memory issue * `normalize()`: Required by interface `IStateWithNormalize` This method allows the evolvers to continually normalize the state, which can improve numerical performance, or be of use when finding initial states. The reset of the methods in the code play a supporting role to these methods required by the interfaces and for our convenience (for example, plotting the state). :param x_TF: Thomas Fermi "radius" for setting the initial state. The initial state will be set so that the density will fall to zero at this point `x=x_TF` (with y and z in the middle of the trap). If `None`, then we default to 80% of the length. :type x_TF: float .. py:attribute:: _use_fftw :value: True .. py:attribute:: t :value: 0 .. py:attribute:: m :value: 86.909187 .. py:attribute:: mu :value: None .. py:attribute:: hbar :value: 1.0 .. py:attribute:: Nxyz .. py:attribute:: Lxyz .. py:attribute:: ws .. py:attribute:: x_TF :value: None .. py:attribute:: g .. py:attribute:: cooling_phase :value: 1.0 .. py:attribute:: symmetric_grid :value: False .. py:attribute:: data .. py:method:: init() Initialize the state. This method defines the basis positions, momenta, etc. for use later on. We define these here rather than in the constructor `__init__()` so that the user can change them later and the reinitialize the state. We also call this function from the `pre_evolve_hook()` so that it is called before any evolution takes place. For this reason, we should not modify the state here. .. py:property:: E_max Return the maximum kinetic energy in the basis. This is useful when using evolvers as convergence should be obtained when the time-step is roughly:: dt = 0.1 * state.hbar / state.E_max See `t_scale`. .. py:property:: t_scale Return the smallest time-scale for the problem. Evolvers - especially the ABM evolvers - should use a `dt=0.1*t_scale` or so. If much smaller `dt` values are required for convergence, then it usually indicates that your lattice spacing is too larger. Likewise, if you can get away with much larger `dt` values, then your lattice spacing might be unnecessarily small. Required by `Simulation`. .. py:method:: pre_evolve_hook() This method is called by the evolvers at the start of evolution to ensure that the state is properly initialized. .. py:method:: set_initial_state(mu=None, x_TF=None) Set the state using the Thomas Fermi (TF) approximation from either `mu` or `x_TF` (pick only one or the other). :param mu: Fixed chemical potential. :type mu: float :param x_TF: Position defining the Thomas Fermi "radius". (The external potential is evaluated at this position and this is used to set `mu`.) :type x_TF: float .. py:method:: get_n_TF(mu, V_ext=None) Return the Thomas Fermi density profile n from mu. .. py:method:: get_mu_TF(x_TF, V_ext=None) Return the Thomas Fermi chemical potential at x_TF. :param x_TF: Position defining the Thomas Fermi "radius". (The external potential is evaluated at this position and this is used to get `mu`.) :type x_TF: float .. py:method:: get_mu_HO(N) Return the chemical potential required to set the particle number in the Thomas Fermi (TF) approximation for a 3D harmonic oscillator. .. py:property:: dim Dimension of the state. .. py:property:: shape Shape of the state Nxyz. .. py:method:: get_Vext() Return the external potential. Overload this method if you want to change the external potential. If the potential should be time dependent, use `self.t` which will be updated by the evolvers. If a chemical potential `self.mu` is defined, then this is subtracted from `Vext`. This allows `gpe.minimize` to find states at a constant chemical potential. Note: for general evolution, it is better not to set the chemical potential as this is automatically set by `compute_dy_dt` and will then cause `get_energy` to return the thermodynamic potential instead. .. py:method:: get_density() Return the density of the state. .. py:method:: get_N() Return the total D-dimensional particle number in the state. .. py:method:: integrate(a) Return the integral of `a` over the box. .. py:method:: braket(a, b) Return the inner product of the functions `a` and `b`. Note: `a` is conjugated and the metric is applied. .. py:method:: get_V() Return the complete potential `V` - internal and external. This function defines the non-linear interaction of the equations. For the GPE, the energy-density has $gn^2/2$, so we have the derivative `gn` here. .. py:method:: fft(y) Return the FFT of y over the spatial indices. .. py:method:: ifft(y) Return the IFFT of y over the spatial indices. .. py:method:: _fftn(y) Return the FFT of y over the spatial indices. These versions with an underscore may be overloaded for performance. .. py:method:: _ifftn(y) Return the IFFT of y over the spatial indices. These versions with an underscore may be overloaded for performance. .. py:method:: compute_dy_dt(dy, subtract_mu=True) Return `dy_dt` storing the results in `dy`. :param subtract_mu: If `True`, then subtract the chemical potential such that `dy_dt` is orthogonal to the original state `y`. This will minimize the evolution of the overall phase during real-time evolution (which can reduce numerical errors) and will ensure that evolution under imaginary or complex time will preserve particle number. This should not be set if computing physical energy of the state, however, which is why it is a parameter. :type subtract_mu: bool .. py:attribute:: linear :value: False .. py:method:: apply_exp_K(dt) Apply $e^{-i K dt}$ in place. .. py:method:: apply_exp_V(dt, state) Apply $e^{-i V dt}$ in place using `state` for any nonlinear dependence in V. (Linear problems should ignore `state`.) .. py:method:: normalize(s=None) Normalize the state, return the scale factors and number `(s, N)`. .. py:method:: get_energy() Return the energy of the state. Useful for minimization. .. py:method:: get_energy_density() Return the energy density. .. py:method:: get_H(subtract_mu=False) Return the single-particle Hamiltonian for convenience only. .. note:: May be huge! This is a square matrix whose sides are of dimension `prod(Nxyz)` which can be huge. It is not needed for general computations, but might be useful for analysis such as exploring the BdG in 1D systems. .. py:method:: get_mu() Compute the chemical potential for convenience only. .. py:method:: get_Hy(subtract_mu=False) Return `H(y)` for convenience only. .. py:method:: plot(log=False, label=None)