Welcome to XMDS2!¶

Installers¶

If you’re using a Mac with OSX >= 10.6, there is a simple drag-and-drop installer available.

If you are using a recent version of Ubuntu or Debian, XMDS2 is available as a package in your package / software manager. Note that if using the Ubuntu Software Centre rather than Synaptic, you must explicitly search for “xmds2” rather than “xmds”, or you will only get XMDS1 as a result.

If you’re using a version of Linux that’s not Ubuntu / Debian, the easiest way to get started is to use the install shell script linked to below.

If we don’t have an installer for your system, follow the manual installation instructions.

If you have one of the supported operating systems listed above, but you find the installer doesn’t work for you, please let us know by emailing xmds-devel <at> lists.sourceforge.net. If you’d like to tweak the linux installer to work on a distribution we haven’t tested, we’d love you to do that and let us know!

Linux installer instructions¶

Note: If you’re using Debian / Ubuntu, unless you want a developer install, you’re probably better off installing XMDS2 via your package manager rather than this install script.

The linux installer has currently only been tested with Ubuntu, Debian, Fedora, and Red Hat. Download the installer here: http://svn.code.sf.net/p/xmds/code/trunk/xpdeint/admin/linux_installer.sh

Once you have downloaded it, make the installer executable and run it by typing the following into a terminal:

chmod u+x linux_installer.sh
./linux_installer.sh

Alternatively, if you wish to download and run the installer in a single step, you can use the following command:

/bin/bash -c "$(wget -qO - http://svn.code.sf.net/p/xmds/code/trunk/xpdeint/admin/linux_installer.sh)"

The linux installer installs all XMDS2 dependencies from your native package manager where possible (apt-get for Ubuntu/Debian, yum for Fedora/Red Hat) but will download and compile the source code for libraries not available through the package manager. This means you’ll need to be connected to the internet when running the installer. The installer should not be run with administrative privileges; it will ask you to enter your admin password at the appropriate point.

For instructions on how to install XMDS2 on systems where you lack administrative rights, see Manual installation from source.

By default, this installer will install a known stable version of XMDS, which can be updated at any time by navigating to the XMDS directory and typing ‘make update’. To install the latest developer version at the beginning, simply run the installer with the --develop option.

Once XMDS2 has been installed, you can run it from the terminal by typing xmds2. See the Quickstart Tutorial for next steps.

Mac OS X Installation¶

Manual installation from source¶

This installation guide will take you through a typical full install step by step. A large part of this procedure is obtaining and installing other libraries that XMDS2 requires, before installing XMDS2 itself.

While the instructions below detail these packages individually, if you have administrative privileges (or can request packages from your administrator) and if you are using an Ubuntu, Debian, Fedora or Red Hat linux distribution, you can install all required and optional dependencies (but not XMDS2 itself) via

Ubuntu / Debian:

sudo apt-get install build-essential subversion libopenmpi-dev openmpi-bin python-dev python-setuptools python-cheetah python-numpy python-pyparsing python-lxml python-mpmath libhdf5-serial-dev libgsl0-dev python-sphinx python-h5py libatlas-base-dev

Fedora / Red Hat:

sudo yum install gcc gcc-c++ make automake subversion openmpi-devel python-devel python-setuptools python-cheetah numpy gsl-devel python-sphinx libxml2-devel libxslt-devel atlas-devel hdf5-devel pyparsing pyparsing python-lxml python-mpmath h5py

You will still have to download and build FFTW 3.3 from source (see below) since prebuilt packages with MPI and AVX support are not currently available in the repositories.

Also note that this guide adds extra notes for users wishing to install XMDS2 using the SVN repository. This requires a few extra steps, but allows you to edit your copy, and/or update your copy very efficiently (with all the usual advantages and disadvantages of using unreleased material).

0. You will need a copy of XMDS2.

   The current release can be found at Sourceforge, and downloaded as a single file. Download this file, and expand it in a directory where you want to keep the program files.

   • Developer-only instructions: You can instead check out a working copy of the source using SVN. In a directory where you want to check out the repository, run: svn checkout https://svn.code.sf.net/p/xmds/code/trunk/xpdeint . (Only do this once. To update your copy, type svn up or make update in the same directory, and then repeat any developer-only instructions below). A checkout with read/write permissions requires a checkout with your login details. e.g.: svn checkout --username=myusername https://myusername@svn.code.sf.net/p/xmds/code/trunk/xpdeint .

1. You will need a working C++ compiler.

   For Mac OS X, this means that the developer tools (XCode) should be installed. One common free compiler is gcc. It can be downloaded using your favourite package manager. XMDS2 can also use Intel’s C++ compiler if you have it. Intel’s compiler typically generates faster code than gcc, but it isn’t free.

2. You will need a python distribution.

   • Mac OS X: It is pre-installed on Mac OS X 10.5 or later.
   • Linux: It should be pre-installed. If not, install using your favourite package manager.

   We require python 2.4 or greater. XMDS2 does not support Python 3.

3. Install setuptools.

   If you have root (sudo) access, the easy way to install this is by executing ez_setup.py from the repository. Simply type sudo python ez_setup.py

   If you want to install into your home directory without root access, this is more complex:

   1. First create the path ~/lib/python2.5/site-packages (assuming you installed python version 2.5) and ~/bin Add “export PYTHONPATH=~/lib/python2.5/site-packages:$PYTHONPATH” and “export PATH=~/bin:$PATH” (if necessary) to your .bashrc file (and run ”. ~/.bashrc”)
   2. If necessary install setuptools, by executing ez_setup.py from the repository. python ez_setup.py --prefix=~

   If you use Mac OS X 10.5 or later, or installed the Enthought Python Distribution on Windows, then setuptools is already installed. Though if the next step fails, you may need to upgrade setuptools. To do that, type sudo easy_install -U setuptools

4. Install HDF5 and FFTW3 (and optionally MPI).

   1. HDF5 is a library for reading and writing the Hierarchical Data Format.

      This is a standardised data format which it is suggested that people use in preference to the older ‘binary’ output (which is compatible with xmds-1). The advantage of HDF5 is that this data format is understood by a variety of other tools. xsil2graphics2 provides support for loading data created in this format into Mathematica and Matlab.

      XMDS2 only requires the single process version of HDF5, so there is no need to install the MPI version.

      * Sidebar: Installing HDF5 from source follows a common pattern, which you may find yourself repeating later:

      1. After extracting the source directory, type configure and then add possible options.

         (For HDF5, install with the --prefix=/usr/local/ option if you want XMDS2 to find the library automatically. This is rarely needed for other packages.)

      2. Once that is finished, type make. Then wait for that to finish, which will often be longer than you think.

      3. Finally, type sudo make install to install it into the appropriate directory.

   2. FFTW is the library XMDS2 uses for Fourier transforms.

      This is the transform most people will use in their simulations. If you need support for MPI distributed simulations, you must configure FFTW to use MPI.

      FFTW is available for free at the FFTW website. To configure and compile it, follow the steps described in the HDF5 sidebar above. You may wish to add the --enable-mpi --disable-fortran options to the configure command.

   3. MPI is an API for doing parallel processing.

      XMDS2 can use MPI to parallelise simulations on multi-processor/multi-core computers, or clusters of computers. Many supercomputing systems come with MPI libraries pre-installed. The Open MPI project has free distributions of this library available.

      If you intend to take advantage of XMDS2’s multi-processing features, you must install MPI, and configure FFTW3 to use it.

5. There are a range of optional installs. We recommend that you install them all if possible:

   1. A Matrix library like ATLAS, Intel’s MKL or the GNU Scientific library (GSL)

      These libraries allow efficient implementation of transform spaces other than Fourier space. Mac OS X comes with its own (fast) matrix library.

   2. numpy is a tool that XMDS2 uses for automated testing.

      It can be installed with sudo easy_install numpy.

      Mac OS X 10.5 and later come with numpy.

   3. lxml is used to validate the syntax of scripts passed to XMDS2.

      If you have root access, this can be installed with the command sudo easy_install lxml

      You will need to have ‘libxml2’ and ‘libxslt’ installed (via your choice of package manager) to install lxml. Sufficient versions are preinstalled on Mac OS X 10.6.

      If you don’t have root access or want to install into your home directory, use:

      easy_install --prefix=~ lxml

   4. h5py is needed for checking the results of XMDS2 tests that generate HDF5 output.

      h5py requires numpy version 1.0.3 or later.

      Upgrading h5py on Mac OS X is best done with the source of the package, as the easy_install option can get confused with multiple numpy versions. (Mac OS X Snow Leopard comes with version 1.2.1). After downloading the source, execute python ./setup.py build in the source directory, and then python ./setup.py install to install it.

6. Install XMDS2 into your python path by running (in the xmds-2.2.1/ directory):

   sudo ./setup.py develop

   If you want to install it into your home directory, type ./setup.py develop --prefix=~

   This step requires access to the net, as it downloads any dependent packages. If you are behind a firewall, you may need to set your HTTP_PROXY environment variable in order to do this.

   • Developer only instructions:

     The Cheetah templates (*.tmpl) must be compiled into python. To do this, run make in the xmds-2.2.1/ directory.

   • Developer-only instructions:

     If you have ‘numpy’ installed, test XMDS2 by typing ./run_tests.py in the xmds-2.2.1/ directory. The package ‘numpy’ is one of the optional packages, with installation instructions below.

   • Developer-only instructions:

     To build the user documentation, you first need to install sphinx, either via your package manager or: sudo easy_install Sphinx

     Then, to build the documentation, in the xmds-2.2.1/admin/userdoc-source/ directory run: make html

     If this results in an error, you may need to run sudo ./setup.py develop

     The generated html documentation will then be found at xmds-2.2.1/documentation/index.html

7. Configure XMDS2 by typing xmds2 --reconfigure. If XMDS2 is unable to find a library, you can tell XMDS2 where these libraries are located by adding include and lib search paths using the --include-path and --lib-path options. For example, if FFTW3 is installed in /apps/fftw3 with headers in /apps/fftw3/include/ and the libraries in /apps/fftw3/lib, (re)configure XMDS2 by typing:

   • xmds2 --reconfigure --include-path /apps/fftw3/include --lib-path /apps/fftw3/lib.

   If you need to use additional compiler or link flags for XMDS2 to use certain libraries, set the CXXFLAGS or LINKFLAGS environment variables before calling xmds2 --reconfigure. For example, to pass the compiler flag -pedantic and the link flag -lm, use:

   • CXXFLAGS="-pedantic" LINKFLAGS="-lm" xmds2 --reconfigure.

Congratulations! You should now have a fully operational copy of xmds2 and xsil2graphics2. You can test your copy using examples from the “xmds-2.2.1/examples” directory, and follow the worked examples in the Quickstart Tutorial and Worked Examples.

The nonlinear Schrödinger equation¶

This worked example will show a range of new features that can be used in an XMDS2 script, and we will also examine our first partial differential equation. We will take the one dimensional nonlinear Schrödinger equation, which is a common nonlinear wave equation. The equation describing this problem is:

where \(\phi\) is a complex-valued field, and \(\Gamma(\tau)\) is a \(\tau\)-dependent damping term. Let us look at an XMDS2 script that integrates this equation, and then examine it in detail.

<simulation xmds-version="2">
  <name>nlse</name>

  <author>Joe Hope</author>
  <description>
    The nonlinear Schrodinger equation in one dimension,
    which is a simple partial differential equation.
    We introduce several new features in this script.
  </description>

  <features>
      <benchmark />
      <bing />
      <fftw plan="patient" />
      <openmp />
      <auto_vectorise />
      <globals>
          <![CDATA[
          const double energy = 4;
          const double vel = 0.3;
          const double hwhm = 1.0;
          ]]>
       </globals>
     </features>

  <geometry>
      <propagation_dimension> xi </propagation_dimension>
      <transverse_dimensions>
        <dimension name="tau" lattice="128"  domain="(-6, 6)" />
      </transverse_dimensions>
   </geometry>

  <vector name="wavefunction" type="complex" dimensions="tau">
    <components> phi </components>
    <initialisation>
      <![CDATA[
      const double w0 = hwhm*sqrt(2/log(2));
      const double amp = sqrt(energy/w0/sqrt(M_PI/2));
      phi = amp*exp(-tau*tau/w0/w0)*exp(i*vel*tau);
      ]]>
    </initialisation>
  </vector>

  <vector name="dampingVector" type="real">
    <components> Gamma </components>
    <initialisation>
      <![CDATA[
      Gamma=1.0*(1-exp(-pow(tau*tau/4.0/4.0,10)));
      ]]>
    </initialisation>
  </vector>

  <sequence>
    <integrate algorithm="ARK45" interval="20.0" tolerance="1e-7">
      <samples>10 100 10</samples>
      <operators>
        <integration_vectors>wavefunction</integration_vectors>
        <operator kind="ex">
          <operator_names>Ltt</operator_names>
          <![CDATA[
            Ltt = -i*ktau*ktau*0.5;
          ]]>
        </operator>
        <![CDATA[
        dphi_dxi = Ltt[phi] - phi*Gamma + i*mod2(phi)*phi;
        ]]>
        <dependencies>dampingVector</dependencies>
      </operators>
    </integrate>
  </sequence>

  <output>
    <sampling_group basis="tau" initial_sample="yes">
      <moments>density</moments>
      <dependencies>wavefunction</dependencies>
      <![CDATA[
        density = mod2(phi);
      ]]>
    </sampling_group>

    <sampling_group basis="tau(0)" initial_sample="yes">
      <moments>normalisation</moments>
      <dependencies>wavefunction</dependencies>
      <![CDATA[
        normalisation = mod2(phi);
      ]]>
    </sampling_group>

    <sampling_group basis="ktau(32)" initial_sample="yes">
      <moments>densityK</moments>
      <dependencies>wavefunction</dependencies>
      <![CDATA[
        densityK = mod2(phi);
      ]]>
    </sampling_group>

  </output>
</simulation>

Let us examine the new items in the <features> element that we have demonstrated here. The existence of the <benchmark> element causes the simulation to be timed. The <bing> element causes the computer to make a sound upon the conclusion of the simulation. The <fftw> element is used to pass options to the FFTW libraries for fast Fourier transforms, which are needed to do spectral derivatives for the partial differential equation. Here we used the option plan=”patient”, which makes the simulation test carefully to find the fastest method for doing the FFTs. More information on possible choices can be found in the FFTW documentation.

Finally, we use two tags to make the simulation run faster. The <auto_vectorise> element switches on several loop optimisations that exist in later versions of the GCC compiler. The <openmp> element turns on threaded parallel processing using the OpenMP standard where possible. These options are not activated by default as they only exist on certain compilers. If your code compiles with them on, then they are recommended.

Let us examine the <geometry> element.

<geometry>
    <propagation_dimension> xi </propagation_dimension>
    <transverse_dimensions>
      <dimension name="tau" lattice="128"  domain="(-6, 6)" />
    </transverse_dimensions>
 </geometry>

This is the first example that includes a transverse dimension. We have only one dimension, and we have labelled it “tau”. It is a continuous dimension, but only defined on a grid containing 128 points (defined with the lattice variable), and on a domain from -6 to 6. The default is that transforms in continuous dimensions are fast Fourier transforms, which means that this dimension is effectively defined on a loop, and the “tau=-6” and “tau=6” positions are in fact the same. Other transforms are possible, as are discrete dimensions such as an integer-valued index, but we will leave these advanced possibilities to later examples.

Two vector elements have been defined in this simulation. One defines the complex-valued wavefunction “phi” that we wish to evolve. We define the transverse dimensions over which this vector is defined by the dimensions tag in the description. By default, it is defined over all of the transverse dimensions in the <geometry> element, so even though we have omitted this tag for the second vector, it also assumes that the vector is defined over all of tau.

The second vector element contains the component “Gamma” which is a function of the transverse variable tau, as specified in the equation of motion for the field. This second vector could have been avoided in two ways. First, the function could have been written explicitly in the integrate block where it is required, but calculating it once and then recalling it from memory is far more efficient. Second, it could have been included in the “wavefunction” vector as another component, but then it would have been unnecessarily complex-valued, it would have needed an explicit derivative in the equations of motion (presumably dGamma_dxi = 0;), and it would have been Fourier transformed whenever the phi component was transformed. So separating it as its own vector is far more efficient.

The <integrate> element for a partial differential equation has some new features:

<integrate algorithm="ARK45" interval="20.0" tolerance="1e-7">
  <samples>10 100 10</samples>
  <operators>
    <integration_vectors>wavefunction</integration_vectors>
    <operator kind="ex">
      <operator_names>Ltt</operator_names>
      <![CDATA[
        Ltt = -i*ktau*ktau*0.5;
      ]]>
    </operator>
    <![CDATA[
    dphi_dxi = Ltt[phi] - phi*Gamma + i*mod2(phi)*phi;
    ]]>
    <dependencies>dampingVector</dependencies>
  </operators>
</integrate>

There are some trivial changes from the tutorial script, such as the fact that we are using the ARK45 algorithm rather than ARK89. Higher order algorithms are often better, but not always. Also, since this script has multiple output groups, we have to specify how many times each of these output groups are sampled in the <samples> element, so there are three numbers there. Besides the vectors that are to be integrated, we also specify that we want to use the vector “dampingVector” during this integration. This is achieved by including the <dependencies> element inside the <operators> element.

The equation of motion as written in the CDATA block looks almost identical to our desired equation of motion, except for the term based on the second derivative, which introduces an important new concept. Inside the <operators> element, we can define any number of operators. Operators are used to define functions in the transformed space of each dimension, which in this case is Fourier space. The derivative of a function is equivalent to multiplying by \(i*k\) in Fourier space, so the \(\frac{i}{2}\frac{\partial^2 \phi}{\partial \tau^2}\) term in our equation of motion is equivalent to multiplying by \(-\frac{i}{2}k_\tau^2\) in Fourier space. In this example we define “Ltt” as an operator of exactly that form, and in the equation of motion it is applied to the field “phi”.

Operators can be explicit (kind="ex") or in the interaction picture (kind="ip"). The interaction picture can be more efficient, but it restricts the possible syntax of the equation of motion. Safe utilisation of interaction picture operators will be described later, but for now let us emphasise that explicit operators should be used unless the user is clear what they are doing. That said, XMDS2 will generate an error if the user tries to use interaction picture operators incorrectly. The constant="yes" option in the operator block means that the operator is not a function of the propagation dimension “xi”, and therefore only needs to be calculated once at the start of the simulation.

The output of a partial differential equation offers more possibilities than an ordinary differential equation, and we examine some in this example.

For vectors with transverse dimensions, we can sample functions of the vectors on the full lattice or a subset of the points. In the <sampling_group> element, we must add a string called “basis” that determines the space in which each transverse dimension is to be sampled, optionally followed by the number of points to be sampled in parentheses. If the number of points is not specified, it will default to a complete sampling of all points in that dimension. If a non-zero number of points is specified, it must be a factor of the lattice size for that dimension.

<sampling_group basis="tau" initial_sample="yes">
  <moments>density</moments>
  <dependencies>wavefunction</dependencies>
  <![CDATA[
    density = mod2(phi);
  ]]>
</sampling_group>

The first output group samples the mod square of the vector “phi” over the full lattice of 128 points.

If the lattice parameter is set to zero points, then the corresponding dimension is integrated.

<sampling_group basis="tau(0)" initial_sample="yes">
  <moments>normalisation</moments>
  <dependencies>wavefunction</dependencies>
  <![CDATA[
    normalisation = mod2(phi);
  ]]>
</sampling_group>

This second output group samples the normalisation of the wavefunction \(\int d\tau |\phi(\tau)|^2\) over the domain of \(\tau\). This output requires only a single real number per sample, so in the integrate element we have chosen to sample it many more times than the vectors themselves.

Finally, functions of the vectors can be sampled with their dimensions in Fourier space.

<sampling_group basis="ktau(32)" initial_sample="yes">
  <moments>densityK</moments>
  <dependencies>wavefunction</dependencies>
  <![CDATA[
    densityK = mod2(phi);
  ]]>
</sampling_group>

The final output group above samples the mod square of the Fourier-space wavefunction phi on a sample of 32 points.

Kubo Oscillator¶

This example demonstrates the integration of a stochastic differential equation. We examine the Kubo oscillator, which is a complex variable whose phase is evolving according to a Wiener noise. In a suitable rotating frame, the equation of motion for the variable is

where \(\eta(t)\) is the Wiener differential, and we interpret this as a Stratonovich equation. In other common notation, this is sometimes written:

Most algorithms employed by XMDS require the equations to be input in the Stratonovich form. Ito differential equations can always be transformed into Stratonovich equations, and in this case the difference is equivalent to the choice of rotating frame. This equation is solved by the following XMDS2 script:

<simulation xmds-version="2">
  <name>kubo</name>
  <author>Graham Dennis and Joe Hope</author>
  <description>
    Example Kubo oscillator simulation
  </description>

  <geometry>
    <propagation_dimension> t </propagation_dimension>
  </geometry>

  <driver name="multi-path" paths="10000" />

  <features>
    <error_check />
    <benchmark />
  </features>

  <noise_vector name="drivingNoise" dimensions="" kind="wiener" type="real" method="dsfmt" seed="314 159 276">
    <components>eta</components>
  </noise_vector>

  <vector name="main" type="complex">
    <components> z </components>
    <initialisation>
      <![CDATA[
        z = 1.0;
      ]]>
    </initialisation>
  </vector>

  <sequence>
    <integrate algorithm="SI" interval="10" steps="1000">
      <samples>100</samples>
      <operators>
        <integration_vectors>main</integration_vectors>
        <dependencies>drivingNoise</dependencies>
        <![CDATA[
          dz_dt = i*z*eta;
        ]]>
      </operators>
    </integrate>
  </sequence>

  <output>
    <sampling_group initial_sample="yes">
      <moments>zR zI</moments>
      <dependencies>main</dependencies>
      <![CDATA[
        zR = z.Re();
        zI = z.Im();
      ]]>
    </sampling_group>
  </output>
</simulation>

The first new item in this script is the <driver> element. This element enables us to change top level management of the simulation. Without this element, XMDS2 will integrate the stochastic equation as described. With this element and the option name="multi-path", it will integrate it multiple times, using different random numbers each time. The output will then contain the mean values and standard errors of your output variables. The number of integrations included in the averages is set with the paths variable.

In the <features> element we have included the <error_check> element. This performs the integration first with the specified number of steps (or with the specified tolerance), and then with twice the number of steps (or equivalently reduced tolerance). The output then includes the difference between the output variables on the coarse and the fine grids as the ‘error’ in the output variables. This error is particularly useful for stochastic integrations, where algorithms with adaptive step-sizes are less safe, so the number of integration steps must be user-specified.

We define the stochastic elements in a simulation with the <noise_vector> element.

<noise_vector name="drivingNoise" dimensions="" kind="wiener" type="real" method="dsfmt" seed="314 159 276">
 <components>eta</components>
</noise_vector>

This defines a vector that is used like any other, but it will be randomly generated with particular statistics and characteristics rather than initialised. The name, dimensions and type tags are defined just as for normal vectors. The names of the components are also defined in the same way. The noise is defined as a Wiener noise here (kind = "wiener"), which is a zero-mean Gaussian random noise with an average variance equal to the discretisation volume (here it is just the step size in the propagation dimension, as it is not defined over transverse dimensions). Other noise types are possible, including uniform and Poissonian noises, but we will not describe them in detail here.

We may also define a noise method to choose a non-default pseudo random number generator, and a seed for the random number generator. Using a seed can be very useful when debugging the behaviour of a simulation, and many compilers have pseudo-random number generators that are superior to the default option (posix).

The integrate block is using the semi-implicit algorithm (algorithm="SI"), which is a good default choice for stochastic problems, even though it is only second order convergent for deterministic equations. More will be said about algorithm choice later, but for now we should note that adaptive algorithms based on Runge-Kutta methods are not guaranteed to converge safely for stochastic equations. This can be particularly deceptive as they often succeed, particularly for almost any problem for which there is a known analytic solution.

We include elements from the noise vector in the equation of motion just as we do for any other vector. The default SI and Runge-Kutta algorithms converge to the Stratonovich integral. Ito stochastic equations can be converted to Stratonovich form and vice versa.

Executing the generated program ‘kubo’ gives slightly different output due to the “multi-path” driver.

$ ./kubo
Beginning full step integration ...
Starting path 1
Starting path 2

... many lines omitted ...

Starting path 9999
Starting path 10000
Beginning half step integration ...
Starting path 1
Starting path 2

... many lines omitted ...

Starting path 9999
Starting path 10000
Generating output for kubo
Maximum step error in moment group 1 was 4.942549e-04
Time elapsed for simulation is: 2.71 seconds

The maximum step error in each moment group is given in absolute terms. This is the largest difference between the full step integration and the half step integration. While a single path might be very stochastic:

The mean value of the real and imaginary components of the z variable for a single path of the simulation.

The average over multiple paths can be increasingly smooth.

The mean and standard error of the z variable averaged over 10000 paths, as given by this simulation. It agrees within the standard error with the expected result of \(\exp(-t/2)\).

Fibre Noise¶

This simulation is a stochastic partial differential equation, in which a one-dimensional damped field is subject to a complex noise. This script can be found in examples/fibre.xmds.

where the noise terms \(\eta_j(x,t)\) are Wiener differentials and the equation is interpreted as a Stratonovich differential equation. On a finite grid, these increments have variance \(\frac{1}{\Delta x \Delta t}\).

<simulation xmds-version="2">
  <name>fibre</name>
  <author>Joe Hope and Graham Dennis</author>
  <description>
    Example fibre noise simulation
  </description>

  <geometry>
    <propagation_dimension> t </propagation_dimension>
    <transverse_dimensions>
      <dimension name="x" lattice="64"  domain="(-5, 5)" />
    </transverse_dimensions>
  </geometry>

  <driver name="mpi-multi-path" paths="8" />

  <features>
    <auto_vectorise />
    <benchmark />
    <error_check />
    <globals>
      <![CDATA[
      const real ggamma = 1.0;
      const real beta = sqrt(M_PI*ggamma/10.0);
      ]]>
    </globals>
  </features>

  <noise_vector name="drivingNoise" dimensions="x" kind="wiener" type="complex" method="dsfmt" seed="314 159 276">
    <components>Eta</components>
  </noise_vector>

  <vector name="main" initial_basis="x" type="complex">
    <components>phi</components>
    <initialisation>
      <![CDATA[
        phi = 0.0;
      ]]>
    </initialisation>
  </vector>

  <sequence>
    <integrate algorithm="SI" iterations="3" interval="2.5" steps="200000">
      <samples>50</samples>
      <operators>
        <operator kind="ex">
          <operator_names>L</operator_names>
          <![CDATA[
            L = -i*kx*kx;
          ]]>
        </operator>
        <dependencies>drivingNoise</dependencies>
        <integration_vectors>main</integration_vectors>
        <![CDATA[
          dphi_dt = L[phi] - ggamma*phi + beta*Eta;
        ]]>
      </operators>
    </integrate>
  </sequence>

  <output>
    <sampling_group basis="kx" initial_sample="yes">
      <moments>pow_dens</moments>
      <dependencies>main</dependencies>
      <![CDATA[
        pow_dens = mod2(phi);
      ]]>
    </sampling_group>
  </output>
</simulation>

Note that the noise vector used in this example is complex-valued, and has the argument dimensions="x" to define it as a field of delta-correlated noises along the x-dimension.

This simulation demonstrates the ease with which XMDS2 can be used in a parallel processing environment. Instead of using the stochastic driver “multi-path”, we simply replace it with “mpi-multi-path”. This instructs XMDS2 to write a parallel version of the program based on the widespread MPI standard. This protocol allows multiple processors or clusters of computers to work simultaneously on the same problem. Free open source libraries implementing this standard can be installed on a linux machine, and come standard on Mac OS X. They are also common on many supercomputer architectures. Parallel processing can also be used with deterministic problems to great effect, as discussed in the later example Wigner Function.

Executing this program is slightly different with the MPI option. The details can change between MPI implementations, but as an example:

$xmds2 fibre.xmds
xmds2 version 2.1 "Happy Mollusc" (r2543)
Copyright 2000-2012 Graham Dennis, Joseph Hope, Mattias Johnsson
                    and the xmds team
Generating source code...
... done
Compiling simulation...
... done. Type './fibre' to run.

Note that different compile options (and potentially a different compiler) are used by XMDS2, but this is transparent to the user. MPI simulations will have to be run using syntax that will depend on the MPI implementation. Here we show the version based on the popular open source Open-MPI implementation.

$ mpirun -np 4 ./fibre
Found enlightenment... (Importing wisdom)
Planning for x <---> kx transform... done.
Beginning full step integration ...
Rank[0]: Starting path 1
Rank[1]: Starting path 2
Rank[2]: Starting path 3
Rank[3]: Starting path 4
Rank[3]: Starting path 8
Rank[0]: Starting path 5
Rank[1]: Starting path 6
Rank[2]: Starting path 7
Rank[3]: Starting path 4
Beginning half step integration ...
Rank[0]: Starting path 1
Rank[2]: Starting path 3
Rank[1]: Starting path 2
Rank[3]: Starting path 8
Rank[0]: Starting path 5
Rank[2]: Starting path 7
Rank[1]: Starting path 6
Generating output for fibre
Maximum step error in moment group 1 was 4.893437e-04
Time elapsed for simulation is: 20.99 seconds

In this example we used four processors. The different processors are labelled by their “Rank”, starting at zero. Because the processors are working independently, the output from the different processors can come in a randomised order. In the end, however, the .xsil and data files are constructed identically to the single processor outputs.

The analytic solution to the stochastic averages of this equation is given by

where \(L_x\) is the length of the x domain. We see that a single integration of these equations is quite chaotic:

The momentum space density of the field as a function of time for a single path realisation.

while an average of 1024 paths (change paths="8" to paths="1024" in the <driver> element) converges nicely to the analytic solution:

The momentum space density of the field as a function of time for an average of 1024 paths.

Integer Dimensions¶

This example shows how to handle systems with integer-valued transverse dimensions. We will integrate the following set of equations

where \(x_j\) are complex-valued variables defined on a ring, such that \(j\in \{0,j_{max}\}\) and the \(x_{j_{max}+1}\) variable is identified with the variable \(x_{0}\), and the variable \(x_{-1}\) is identified with the variable \(x_{j_{max}}\).

<simulation xmds-version="2">
  <name>integer_dimensions</name>
  <author>Graham Dennis</author>
  <description>
    XMDS2 script to test integer dimensions.
  </description>

  <features>
    <benchmark />
    <error_check />
    <bing />
    <diagnostics /> <!-- This will make sure that all nonlocal accesses of dimensions are safe -->
  </features>

  <geometry>
    <propagation_dimension> t </propagation_dimension>
    <transverse_dimensions>
      <dimension name="j" type="integer" lattice="5" domain="(0,4)" />
    </transverse_dimensions>
  </geometry>

  <vector name="main" type="complex">
    <components> x </components>
    <initialisation>
      <![CDATA[
      x = 1.0e-3;
      x(j => 0) = 1.0;
      ]]>
    </initialisation>
  </vector>

  <sequence>
    <integrate algorithm="ARK45" interval="60" steps="25000" tolerance="1.0e-9">
      <samples>1000</samples>
      <operators>
        <integration_vectors>main</integration_vectors>
        <![CDATA[
        long j_minus_one = (j-1) % _lattice_j;
        if (j_minus_one < 0)
          j_minus_one += _lattice_j;
        long j_plus_one  = (j+1) % _lattice_j;
        dx_dt(j => j) = x(j => j)*(x(j => j_minus_one) - x(j => j_plus_one));
        ]]>
      </operators>
    </integrate>
  </sequence>

  <output>
    <sampling_group basis="j" initial_sample="yes">
      <moments>xR</moments>
      <dependencies>main</dependencies>
      <![CDATA[
        xR = x.Re();
      ]]>
    </sampling_group>
  </output>
</simulation>

The first extra feature we have used in this script is the <diagnostics> element. It performs run-time checking that our generated code does not accidentally attempt to access a part of our vector that does not exist. Removing this tag will increase the speed of the simulation, but its presence helps catch coding errors.

The simulation defines a vector with a single transverse dimension labelled “j”, of type “integer” (“int” and “long” can also be used as synonyms for “integer”). In the absence of an explicit type, the dimension is assumed to be real-valued. The dimension has a “domain” argument as normal, defining the minimum and maximum values of the dimension’s range. The lattice element, if specified, is used as a check on the size of the domain, and will create an error if the two do not match.

Integer-valued dimensions can be called non-locally. Real-valued dimensions are typically coupled non-locally only through local operations in the transformed space of the dimension, but can be called non-locally in certain other situations as described in the reference. The syntax for calling integer dimensions non-locally can be seen in the initialisation CDATA block:

x = 1.0e-3;
x(j => 0) = 1.0;

where the syntax x(j => 0) is used to reference the variable \(x_0\) directly. We see a more elaborate example in the integrate CDATA block:

dx_dt(j => j) = x(j => j)*(x(j => j_minus_one) - x(j => j_plus_one));

where the vector “x” is called using locally defined variables. This syntax is chosen so that multiple dimensions can be addressed non-locally with minimal possibility for confusion.

Wigner Function¶

This example integrates the two-dimensional partial differential equation

with the added restriction that the derivative is forced to zero outside a certain radius. This extra condition helps maintain the long-term stability of the integration. The script can be found in examples/wigner_arguments_mpi.xmds under your XMDS2 installation directory.

<simulation xmds-version="2">
  <name>wigner</name>
  <author>Graham Dennis and Joe Hope</author>
  <description>
    Simulation of the Wigner function for an anharmonic oscillator with the initial state
    being a coherent state.
  </description>
  <features>
    <benchmark />
    <globals>
      <![CDATA[
        real Uint_hbar_on16;
      ]]>
    </globals>
    <arguments>
      <argument name="omega" type="real" default_value="0.0" />
      <argument name="alpha_0"     type="real" default_value="3.0" />
      <argument name="absorb"     type="real" default_value="8.0" />
      <argument name="width" type="real" default_value="0.3" />
      <argument name="Uint_hbar" type="real" default_value="1.0" />
      <![CDATA[
        /* derived constants */
        Uint_hbar_on16 = Uint_hbar/16.0;
      ]]>
    </arguments>
    <bing />
    <fftw plan="patient" />
    <openmp />
  </features>

  <driver name="distributed-mpi" />

  <geometry>
    <propagation_dimension> t </propagation_dimension>
    <transverse_dimensions>
      <dimension name="x" lattice="128"  domain="(-6, 6)" />
      <dimension name="y" lattice="128"  domain="(-6, 6)" />
    </transverse_dimensions>
  </geometry>

  <vector name="main" initial_basis="x y" type="complex">
    <components> W </components>
    <initialisation>
      <![CDATA[
        W = 2.0/M_PI * exp(-2.0*(y*y + (x-alpha_0)*(x-alpha_0)));
      ]]>
    </initialisation>
  </vector>

  <vector name="dampConstants" initial_basis="x y" type="real">
    <components>damping</components>
    <initialisation>
      <![CDATA[
      if (sqrt(x*x + y*y) > _max_x-width)
        damping = 0.0;
      else
        damping = 1.0;
      ]]>
    </initialisation>
  </vector>

  <sequence>
    <integrate algorithm="ARK89" tolerance="1e-7" interval="7.0e-4" steps="100000">
      <samples>50</samples>
      <operators>
        <operator kind="ex">
          <operator_names>Lx Ly Lxxx Lxxy Lxyy Lyyy</operator_names>
          <![CDATA[
            Lx = i*kx;
            Ly = i*ky;
            Lxxx = -i*kx*kx*kx;
            Lxxy = -i*kx*kx*ky;
            Lxyy = -i*kx*ky*ky;
            Lyyy = -i*ky*ky*ky;
          ]]>
        </operator>
        <integration_vectors>main</integration_vectors>
        <dependencies>dampConstants</dependencies>
        <![CDATA[
        real rotation = omega + Uint_hbar*(-1.0 + x*x + y*y);

        dW_dt = damping * ( rotation * (x*Ly[W] - y*Lx[W])
                    - Uint_hbar_on16*( x*(Lxxy[W] + Lyyy[W]) - y*(Lxyy[W] + Lxxx[W]) )
                );
        ]]>
      </operators>
    </integrate>
  </sequence>

  <output>
    <sampling_group basis="x y" initial_sample="yes">
      <moments>WR WI</moments>
      <dependencies>main</dependencies>
      <![CDATA[
        _SAMPLE_COMPLEX(W);
      ]]>
    </sampling_group>
  </output>
</simulation>

This example demonstrates two new features of XMDS2. The first is the use of parallel processing for a deterministic problem. The FFTW library only allows MPI processing of multidimensional vectors. For multidimensional simulations, the generated program can be parallelised simply by adding the name="distributed-mpi" argument to the <driver> element.

$ xmds2 wigner_argument_mpi.xmds
xmds2 version 2.1 "Happy Mollusc" (r2680)
Copyright 2000-2012 Graham Dennis, Joseph Hope, Mattias Johnsson
                    and the xmds team
Generating source code...
... done
Compiling simulation...
... done. Type './wigner' to run.

To use multiple processors, the final program is then called using the (implementation specific) MPI wrapper:

$ mpirun -np 2 ./wigner
Planning for (distributed x, y) <---> (distributed ky, kx) transform... done.
Planning for (distributed x, y) <---> (distributed ky, kx) transform... done.
Sampled field (for moment group #1) at t = 0.000000e+00
Current timestep: 5.908361e-06
Sampled field (for moment group #1) at t = 1.400000e-05
Current timestep: 4.543131e-06

...

The possible acceleration achievable when parallelising a given simulation depends on a great many things including available memory and cache. As a general rule, it will improve as the simulation size gets larger, but the easiest way to find out is to test. The optimum speed up is obviously proportional to the number of available processing cores.

The second new feature in this simulation is the <arguments> element in the <features> block. This is a way of specifying global variables with a given type that can then be input at run time. The variables are specified in a self explanatory way

<arguments>
  <argument name="omega" type="real" default_value="0.0" />
    ...
  <argument name="Uint_hbar" type="real" default_value="1.0" />
</arguments>

where the “default_value” is used as the valuable of the variable if no arguments are given. In the absence of the generating script, the program can document its options with the --help argument:

$ ./wigner --help
Usage: wigner --omega <real> --alpha_0 <real> --absorb <real> --width <real> --Uint_hbar <real>

Details:
Option              Type            Default value
-o,  --omega        real            0.0
-a,  --alpha_0      real            3.0
-b,  --absorb       real            8.0
-w,  --width        real            0.3
-U,  --Uint_hbar    real            1.0

We can change one or more of these variables’ values in the simulation by passing it at run time.

$ mpirun -np 2 ./wigner --omega 0.1 --alpha_0 2.5 --Uint_hbar 0
Found enlightenment... (Importing wisdom)
Planning for (distributed x, y) <---> (distributed ky, kx) transform... done.
Planning for (distributed x, y) <---> (distributed ky, kx) transform... done.
Sampled field (for moment group #1) at t = 0.000000e+00
Current timestep: 1.916945e-04

...

The values that were used for the variables, whether default or passed in, are stored in the output file (wigner.xsil).

<info>
Script compiled with XMDS2 version 2.1 "Happy Mollusc" (r2680)
See http://www.xmds.org for more information.

Variables that can be specified on the command line:
  Command line argument omega = 1.000000e-01
  Command line argument alpha_0 = 2.500000e+00
  Command line argument absorb = 8.000000e+00
  Command line argument width = 3.000000e-01
  Command line argument Uint_hbar = 0.000000e+00
</info>

Finally, note the shorthand used in the output group

<![CDATA[
  _SAMPLE_COMPLEX(W);
]]>

which is short for

<![CDATA[
  WR = W.Re();
  WI = W.Im();
]]>

Finding the Ground State of a BEC (continuous renormalisation)¶

This simulation solves another partial differential equation, but introduces several powerful new features in XMDS2. The nominal problem is the calculation of the lowest energy eigenstate of a non-linear Schrödinger equation:

which can be found by evolving the above equation in imaginary time while keeping the normalisation constant. This causes eigenstates to exponentially decay at the rate of their eigenvalue, so after a short time only the state with the lowest eigenvalue remains. The evolution equation is straightforward:

but we will need to use new XMDS2 features to manage the normalisation of the function \(\phi(y,t)\). The normalisation for a non-linear Schrödinger equation is given by \(\int dy |\phi(y,t)|^2 = N_{particles}\), where \(N_{particles}\) is the number of particles described by the wavefunction.

The code for this simulation can be found in examples/groundstate_workedexamples.xmds:

<simulation xmds-version="2">
  <name>groundstate</name>
  <author>Joe Hope</author>
  <description>
    Calculate the ground state of the non-linear Schrodinger equation in a harmonic magnetic trap.
    This is done by evolving it in imaginary time while re-normalising each timestep.
  </description>

  <features>
    <auto_vectorise />
    <benchmark />
    <bing />
    <fftw plan="exhaustive" />
    <globals>
      <![CDATA[
        const real Uint = 2.0;
        const real Nparticles = 5.0;
      ]]>
    </globals>
  </features>

  <geometry>
    <propagation_dimension> t </propagation_dimension>
    <transverse_dimensions>
      <dimension name="y" lattice="256"  domain="(-15.0, 15.0)" />
    </transverse_dimensions>
  </geometry>

  <vector name="potential" initial_basis="y" type="real">
    <components> V1 </components>
    <initialisation>
      <![CDATA[
        V1  = 0.5*y*y;
      ]]>
    </initialisation>
  </vector>

  <vector name="wavefunction" initial_basis="y" type="complex">
    <components> phi </components>
    <initialisation>
      <![CDATA[
        if (fabs(y) < 3.0) {
          phi = 1.0;
          // This will be automatically normalised later
        } else {
          phi = 0.0;
        }
            ]]>
    </initialisation>
  </vector>

  <computed_vector name="normalisation" dimensions="" type="real">
    <components> Ncalc </components>
    <evaluation>
      <dependencies basis="y">wavefunction</dependencies>
      <![CDATA[
        // Calculate the current normalisation of the wave function.
        Ncalc = mod2(phi);
      ]]>
    </evaluation>
  </computed_vector>

  <sequence>
      <filter>
        <![CDATA[
          printf("Hello world from a filter segment!\n");
        ]]>
      </filter>

    <filter>
        <dependencies>normalisation wavefunction</dependencies>
      <![CDATA[
        phi *= sqrt(Nparticles/Ncalc);
      ]]>
    </filter>

    <integrate algorithm="ARK45" interval="1.0" steps="4000" tolerance="1e-10">
      <samples>25 4000</samples>
      <filters where="step end">
        <filter>
          <dependencies>wavefunction normalisation</dependencies>
          <![CDATA[
            // Correct normalisation of the wavefunction
            phi *= sqrt(Nparticles/Ncalc);
          ]]>
        </filter>
      </filters>
      <operators>
        <operator kind="ip">
          <operator_names>T</operator_names>
          <![CDATA[
            T = -0.5*ky*ky;
          ]]>
        </operator>
        <integration_vectors>wavefunction</integration_vectors>
        <dependencies>potential</dependencies>
        <![CDATA[
          dphi_dt = T[phi] - (V1 + Uint*mod2(phi))*phi;
        ]]>
      </operators>
    </integrate>

    <breakpoint filename="groundstate_break.xsil">
      <dependencies basis="ky">wavefunction </dependencies>
    </breakpoint>

  </sequence>

  <output>
    <sampling_group basis="y" initial_sample="yes">
      <moments>norm_dens</moments>
      <dependencies>wavefunction normalisation</dependencies>
      <![CDATA[
        norm_dens = mod2(phi);
      ]]>
    </sampling_group>

    <sampling_group initial_sample="yes">
      <moments>norm</moments>
      <dependencies>normalisation</dependencies>
      <![CDATA[
        norm = Ncalc;
      ]]>
    </sampling_group>
  </output>
</simulation>

We have used the plan="exhaustive" option in the <fftw> element to ensure that the absolute fastest transform method is found. Because the FFTW package stores the results of its tests (by default in the ~/.xmds/wisdom directory), this option does not cause significant computational overhead, except perhaps on the very first run of a new program.

This simulation introduces the first example of a very powerful feature in XMDS2: the <computed_vector> element. This has syntax like any other vector, including possible dependencies on other vectors, and an ability to be used in any element that can use vectors. The difference is that, much like noise vectors, computed vectors are recalculated each time they are required. This means that a computed vector can never be used as an integration vector, as its values are not stored. However, computed vectors allow a simple and efficient method of describing complicated functions of other vectors. Computed vectors may depend on other computed vectors, allowing for spectral filtering and other advanced options. See for example, the Advanced Topics section on Convolutions and Fourier transforms.

The difference between a computed vector and a stored vector is emphasised by the replacement of the <initialisation> element with an <evaluation> element. Apart from the name, they have virtually identical purpose and syntax.

<computed_vector name="normalisation" dimensions="" type="real">
  <components> Ncalc </components>
  <evaluation>
    <dependencies basis="y">wavefunction</dependencies>
    <![CDATA[
      // Calculate the current normalisation of the wave function.
      Ncalc = mod2(phi);
    ]]>
  </evaluation>
</computed_vector>

Here, our computed vector has no transverse dimensions and depends on the components of “wavefunction”, so the extra transverse dimensions are integrated out. This code therefore integrates the square modulus of the field, and returns it in the variable “Ncalc”. This will be used below to renormalise the “phi” field. Before we examine that process, we have to introduce the <filter> element.

The <filter> element can be placed in the <sequence> element, or inside <integrate> elements as we will see next. Elements placed in the <sequence> element are executed in the order they are found in the .xmds file. Filter elements place the included CDATA block directly into the generated program at the designated position. If the element does not contain any dependencies, like in our first example, then the code is placed alone:

<filter>
  <![CDATA[
    printf("Hello world from a filter segment!\n");
  ]]>
</filter>

This filter block merely prints a string into the output when the generated program is run. If the <filter> element contains dependencies, then the variables defined in those vectors (or computed vectors, or noise vectors) will be available, and the CDATA block will be placed inside loops that run over all the transverse dimensions used by the included vectors. The second filter block in this example depends on both the “wavefunction” and “normalisation” vectors:

<filter>
    <dependencies>normalisation wavefunction</dependencies>
  <![CDATA[
    phi *= sqrt(Nparticles/Ncalc);
  ]]>
</filter>

Since this filter depends on a vector with the transverse dimension “y”, this filter will execute for each point in “y”. This code multiplies the value of the field “phi” by the factor required to produce a normalised function in the sense that \(\int dy |\phi(y,t)|^2 = N_{particles}\).

The next usage of a <filter> element in this program is inside the <integrate> element, where all filters are placed inside a <filters> element.

<filters where="step end">
  <filter>
    <dependencies>wavefunction normalisation</dependencies>
    <![CDATA[
      // Correct normalisation of the wavefunction
      phi *= sqrt(Nparticles/Ncalc);
    ]]>
  </filter>
</filters>

Filters placed in an integration block are applied each integration step. The “where” flag is used to determine whether the filter should be applied directly before or directly after each integration step. The default value for the where flag is where="step start", but in this case we chose “step end” to make sure that the final output was normalised after the last integration step.

At the end of the sequence element we introduce the <breakpoint> element. This serves two purposes. The first is a simple matter of convenience. Often when we manage our input and output from a simulation, we are interested solely in storing the exact state of our integration vectors. A breakpoint element does exactly that, storing the components of any vectors contained within, taking all the normal options of the <output> element but not requiring any <sampling_group> elements as that information is assumed.

<breakpoint filename="groundstate_break.xsil">
  <dependencies basis="ky">wavefunction</dependencies>
</breakpoint>

If the filename argument is omitted, the output filenames are numbered sequentially. Any given <breakpoint> element must only depend on vectors with identical dimensions.

This program begins with a very crude guess to the ground state, but it rapidly converges to the lowest eigenstate.

The shape of the ground state rapidly approaches the lowest eigenstate. For weak nonlinearities, it is nearly Gaussian.

When the nonlinear term is larger (\(U=20\)), the ground state is wider and more parabolic.

Finding the Ground State of a BEC again¶

Here we repeat the same simulation as in the Finding the Ground State of a BEC (continuous renormalisation) example, using a different transform basis. While spectral methods are very effective, and Fourier transforms are typically very efficient due to the Fast Fourier transform algorithm, it is often desirable to describe nonlocal evolution in bases other than the Fourier basis. The previous calculation was the Schrödinger equation with a harmonic potential and a nonlinear term. The eigenstates of such a system are known analytically to be Gaussians multiplied by the Hermite polynomials.

where

where \(H_n(u)\) are the physicist’s version of the Hermite polynomials. Rather than describing the derivatives as diagonal terms in Fourier space, we therefore have the option of describing the entire \(-\frac{\hbar}{2 m}\frac{\partial^2}{\partial x^2} + \frac{1}{2}\omega^2 x^2\) term as a diagonal term in the Hermite-Gaussian basis. Here is an XMDS2 simulation that performs the integration in this basis. The following is a simplified version of the examples/hermitegauss_groundstate.xmds script.

<simulation xmds-version="2">
  <name>hermitegauss_groundstate</name>
  <author>Graham Dennis</author>
  <description>
    Solve for the groundstate of the Gross-Pitaevskii equation using the Hermite-Gauss basis.
  </description>

  <features>
    <benchmark />
    <bing />
    <validation kind="run-time" />
    <globals>
      <![CDATA[
        const real omegaz = 2*M_PI*20;
        const real omegarho = 2*M_PI*200;
        const real hbar = 1.05457148e-34;
        const real M = 1.409539200000000e-25;
        const real g = 9.8;
        const real scatteringLength = 5.57e-9;
        const real transverseLength = 1e-5;
        const real Uint = 4.0*M_PI*hbar*hbar*scatteringLength/M/transverseLength/transverseLength;
        const real Nparticles = 5.0e5;

        /* offset constants */
        const real EnergyOffset = 0.3*pow(pow(3.0*Nparticles/4*omegarho*Uint,2.0)*M/2.0,1/3.0); // 1D
      ]]>
    </globals>
  </features>

  <geometry>
    <propagation_dimension> t </propagation_dimension>
    <transverse_dimensions>
      <dimension name="x" lattice="100" length_scale="sqrt(hbar/(M*omegarho))" transform="hermite-gauss" />
    </transverse_dimensions>
  </geometry>

  <vector name="wavefunction" initial_basis="x" type="complex">
    <components> phi </components>
    <initialisation>
      <![CDATA[
      phi = sqrt(Nparticles) * pow(M*omegarho/(hbar*M_PI), 0.25) * exp(-0.5*(M*omegarho/hbar)*x*x);
      ]]>
    </initialisation>
  </vector>

  <computed_vector name="normalisation" dimensions="" type="real">
    <components> Ncalc </components>
    <evaluation>
      <dependencies basis="x">wavefunction</dependencies>
      <![CDATA[
        // Calculate the current normalisation of the wave function.
        Ncalc = mod2(phi);
      ]]>
    </evaluation>
  </computed_vector>

  <sequence>
    <integrate algorithm="ARK45" interval="1.0e-2" steps="4000"  tolerance="1e-10">
      <samples>100 100</samples>
      <filters>
        <filter>
          <dependencies>wavefunction normalisation</dependencies>
          <![CDATA[
            // Correct normalisation of the wavefunction
            phi *= sqrt(Nparticles/Ncalc);
          ]]>
        </filter>
      </filters>
      <operators>
        <operator kind="ip" type="real">
          <operator_names>L</operator_names>
          <![CDATA[
            L = EnergyOffset/hbar - (nx + 0.5)*omegarho;
          ]]>
        </operator>
        <integration_vectors>wavefunction</integration_vectors>
        <![CDATA[
          dphi_dt = L[phi] - Uint/hbar*mod2(phi)*phi;
        ]]>
      </operators>
    </integrate>

    <filter>
        <dependencies>normalisation wavefunction</dependencies>
      <![CDATA[
        phi *= sqrt(Nparticles/Ncalc);
      ]]>
    </filter>

    <breakpoint filename="hermitegauss_groundstate_break.xsil" format="ascii">
      <dependencies basis="nx">wavefunction</dependencies>
    </breakpoint>
  </sequence>

  <output>
    <sampling_group basis="x" initial_sample="yes">
      <moments>dens</moments>
      <dependencies>wavefunction</dependencies>
      <![CDATA[
        dens = mod2(phi);
      ]]>
    </sampling_group>
    <sampling_group basis="kx" initial_sample="yes">
      <moments>dens</moments>
      <dependencies>wavefunction</dependencies>
      <![CDATA[
        dens = mod2(phi);
      ]]>
    </sampling_group>
  </output>
</simulation>

The major difference in this simulation code, aside from the switch back from dimensionless units, is the new transverse dimension type in the <geometry> element.

<dimension name="x" lattice="100" length_scale="sqrt(hbar/(M*omegarho))" transform="hermite-gauss" />

We have explicitly defined the “transform” option, which by default expects the Fourier transform. The transform="hermite-gauss" option requires the ‘mpmath’ package installed, just as Fourier transforms require the FFTW package to be installed. The “lattice” option details the number of hermite-Gaussian eigenstates to include, and automatically starts from the zeroth order polynomial and increases. The number of hermite-Gaussian modes fully determines the irregular spatial grid up to an overall scale given by the length_scale parameter.

The length_scale="sqrt(hbar/(M*omegarho))" option requires a real number, but since this script defines it in terms of variables, XMDS2 is unable to verify that the resulting function is real-valued at the time of generating the code. XMDS2 will therefore fail to compile this program without the feature:

<validation kind="run-time" />

which disables many of these checks at the time of writing the C-code.

Multi-component Schrödinger equation¶

This example demonstrates a simple method for doing matrix calculations in XMDS2. We are solving the multi-component PDE

where the last term is more commonly written as a matrix multiplication. Writing this term out explicitly is feasible for a small number of components, but when the number of components becomes large, or perhaps \(V_{j k}\) should be precomputed for efficiency reasons, it is useful to be able to perform this sum over the integer dimensions automatically. This example show how this can be done naturally using a computed vector. The XMDS2 script is as follows:

<simulation xmds-version="2">
  <name>2DMSse</name>

  <author>Joe Hope</author>
  <description>
    Schroedinger equation for multiple internal states in two spatial dimensions.
  </description>

  <features>
      <benchmark />
      <bing />
      <fftw plan="patient" />
      <openmp />
      <auto_vectorise />
     </features>

  <geometry>
      <propagation_dimension> t </propagation_dimension>
      <transverse_dimensions>
          <dimension name="x" lattice="32"  domain="(-6, 6)" />
          <dimension name="y" lattice="32"  domain="(-6, 6)" />
          <dimension name="j" type="integer" lattice="2" domain="(0,1)" aliases="k"/>
      </transverse_dimensions>
   </geometry>

  <vector name="wavefunction" type="complex" dimensions="x y j">
    <components> phi </components>
    <initialisation>
      <![CDATA[
      phi = j*sqrt(2/sqrt(M_PI/2))*exp(-(x*x+y*y)/4)*exp(i*0.1*x);
      ]]>
    </initialisation>
  </vector>

  <vector name="spatialInteraction" type="real" dimensions="x y">
    <components> U </components>
    <initialisation>
      <![CDATA[
      U=exp(-(x*x+y*y)/4);
      ]]>
    </initialisation>
  </vector>

  <vector name="internalInteraction" type="real" dimensions="j k">
    <components> V </components>
    <initialisation>
      <![CDATA[
      V=3*(j*(1-k)+(1-j)*k);
      ]]>
    </initialisation>
  </vector>

  <computed_vector name="coupling" dimensions="x y j" type="complex">
    <components>
      VPhi
    </components>
    <evaluation>
      <dependencies basis="x y j k">internalInteraction wavefunction</dependencies>
      <![CDATA[
        // Calculate the current normalisation of the wave function.
        VPhi = V*phi(j => k);
      ]]>
    </evaluation>
  </computed_vector>

  <sequence>
    <integrate algorithm="ARK45" interval="2.0" tolerance="1e-7">
      <samples>20 100</samples>
      <operators>
        <integration_vectors>wavefunction</integration_vectors>
        <operator kind="ip" dimensions="x">
          <operator_names>Lx</operator_names>
          <![CDATA[
            Lx = -i*kx*kx*0.5;
          ]]>
        </operator>
        <operator kind="ip" dimensions="y">
          <operator_names>Ly</operator_names>
          <![CDATA[
            Ly = -i*ky*ky*0.5;
          ]]>
        </operator>
        <![CDATA[
        dphi_dt = Lx[phi] + Ly[phi] -i*U*VPhi;
        ]]>
        <dependencies>spatialInteraction coupling</dependencies>
      </operators>
    </integrate>
  </sequence>

  <output>
    <sampling_group basis="x y j" initial_sample="yes">
      <moments>density</moments>
      <dependencies>wavefunction</dependencies>
      <![CDATA[
        density = mod2(phi);
      ]]>
    </sampling_group>
    <sampling_group basis="x(0) y(0) j" initial_sample="yes">
      <moments>normalisation</moments>
      <dependencies>wavefunction</dependencies>
      <![CDATA[
        normalisation = mod2(phi);
      ]]>
    </sampling_group>
  </output>
</simulation>

The only truly new feature in this script is the “aliases” option on a dimension. The integer-valued dimension in this script indexes the components of the PDE (in this case only two). The \(V_{j k}\) term is required to be a square array of dimension of this number of components. If we wrote the k-index of \(V_{j k}\) using a separate <dimension> element, then we would not be enforcing the requirement that the matrix be square. Instead, we note that we will be using multiple ‘copies’ of the j-dimension by using the “aliases” tag.

<dimension name="j" type="integer" lattice="2" domain="(0,1)" aliases="k"/>

This means that we can use the index “k”, which will have exactly the same properties as the “j” index. This is used to define the “V” function in the “internalInteraction” vector. Now, just as we use a computed vector to perform an integration over our fields, we use a computed vector to calculate the sum.

<computed_vector name="coupling" dimensions="x y j" type="complex">
  <components>
    VPhi
  </components>
  <evaluation>
    <dependencies basis="x y j k">internalInteraction wavefunction</dependencies>
    <![CDATA[
      // Calculate the current normalisation of the wave function.
      VPhi = V*phi(j => k);
    ]]>
  </evaluation>
</computed_vector>

Since the output dimensions of the computed vector do not include a “k” index, this index is integrated. The volume element for this summation is the spacing between neighbouring values of “j”, and since this spacing is one, this integration is just a sum over k, as required.

This example also demonstrates an optimisation for the IP operators by separating the \(x\) and \(y\) parts of the operator (see Optimising with the Interaction Picture (IP) operator). This gives an approximately 30% speed improvement over the more straightforward implementation:

<operator kind="ip">
  <operator_names>L</operator_names>
  <![CDATA[
    L = -i*(kx*kx + ky*ky)*0.5;
  ]]>
</operator>
<![CDATA[
  dphi_dt = L[phi] - i*U*VPhi;
]]>

By this point, we have introduced most of the important features in XMDS2. More details on other transform options and rarely used features can be found in the Advanced Topics section.

Importing data¶

There are many cases where it is advantageous to import previously acquired data into XMDS2. For example, the differential equation you wish to solve may depend on a complicated functional form, which is more easily obtained via an analytical package such as Mathematica or Maple. Furthermore, importing data from another source can be quicker than needlessly performing calculations in XMDS2. In this tutorial, we shall consider an example of importing into XMDS2 a function generated in Mathematica, version 6.0. Note, however, that in order to do this it is required that hdf5 is installed (see http://www.hdfgroup.org/HDF5/).

Suppose we want to import the following function into XMDS2:

The first step is to create an hdf5 file, from XMDS2, which specifies the dimensions of the grid for the x dimension. Create and save a new XMDS2 file. For the purposes of this tutorial we shall call it “grid_specifier.xmds” with name “grid_specifier”. Within this file, enter the following “dummy” vector - which we shall call “gen_dummy” - which depends on the x dimension:

<vector type="real" name="gen_dummy" dimensions="x">
  <components>dummy</components>
  <initialisation>
  <![CDATA[
    dummy = x;
      ]]>
  </initialisation>
</vector>

What “dummy” is is not actually important. It is only necessary that it is a function of \(x\). However, it is important that the domain and lattice for the \(x\) dimension are identical to those in the XMDS2 you plan to import the function into. We output the following xsil file (in hdf5 format) by placing a breakpoint in the sequence block as follows:

<sequence>
  <breakpoint filename="grid.xsil" format="hdf5">
      <dependencies>
        gen_dummy
      </dependencies>
  </breakpoint>

In terminal, compile the file “grid_specifier.xmds” in XMDS2 and run the C code as usual. This creates two files - “grid.xsil” and “grid.h5”. The file “grid.h5” contains the list of points which make up the grids for the x dimensions. This data can now be used to ensure that the function \(f(x)\) which we will import into XMDS2 is compatible with the the specified grid in your primary XMDS2 file.

In order to read the “grid.h5” data into Mathematica version 6.0, type the following command into terminal:.. code-block:

xsil2graphics2 -e grid.xsil

This creates the Mathematica notebook “grid.nb”. Open this notebook in Mathematica and evaluate the first set of cells. This has loaded the grid information into Mathematica. For example, suppose you have specified that the \(x\) dimension has a lattice of 128 points and a domain of (-32, 32). Then calling “x1” in Mathematica should return the following list:

{-32., -31.5, -31., -30.5, -30., -29.5, -29., -28.5, -28., -27.5,
-27., -26.5, -26., -25.5, -25., -24.5, -24., -23.5, -23., -22.5,
-22., -21.5, -21., -20.5, -20., -19.5, -19., -18.5, -18., -17.5,
-17., -16.5, -16., -15.5, -15., -14.5, -14., -13.5, -13., -12.5,
-12., -11.5, -11., -10.5, -10., -9.5, -9., -8.5, -8., -7.5, -7.,
-6.5, -6., -5.5, -5., -4.5, -4., -3.5, -3., -2.5, -2., -1.5, -1.,
-0.5, 0., 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5., 5.5, 6., 6.5,
7., 7.5, 8., 8.5, 9., 9.5, 10., 10.5, 11., 11.5, 12., 12.5, 13.,
13.5, 14., 14.5, 15., 15.5, 16., 16.5, 17., 17.5, 18., 18.5, 19.,
19.5, 20., 20.5, 21., 21.5, 22., 22.5, 23., 23.5, 24., 24.5, 25.,
25.5, 26., 26.5, 27., 27.5, 28., 28.5, 29., 29.5, 30., 30.5, 31.,
31.5}

This is, of course, the list of points which define our grid.

We are now in a position to define the function \(f(x)\) in Mathematica. Type the following command into a cell in the Mathematica notebook “grid.nb”:

f[x_]:= x^2

At this stage this is an abstract mathematical function as defined in Mathematica. What we need is a list of values for \(f(x)\) corresponding to the specified grid points. We will call this list “func”. This achieved by simply acting the function on the list of grid points “x1”:

func := f[x1]

For the example grid mentioned above, calling “func” gives the following list:

{1024., 992.25, 961., 930.25, 900., 870.25, 841., 812.25, 784.,
756.25, 729., 702.25, 676., 650.25, 625., 600.25, 576., 552.25, 529.,
506.25, 484., 462.25, 441., 420.25, 400., 380.25, 361., 342.25, 324.,
306.25, 289., 272.25, 256., 240.25, 225., 210.25, 196., 182.25, 169.,
156.25, 144., 132.25, 121., 110.25, 100., 90.25, 81., 72.25, 64.,
56.25, 49., 42.25, 36., 30.25, 25., 20.25, 16., 12.25, 9., 6.25, 4.,
2.25, 1., 0.25, 0., 0.25, 1., 2.25, 4., 6.25, 9., 12.25, 16., 20.25,
25., 30.25, 36., 42.25, 49., 56.25, 64., 72.25, 81., 90.25, 100.,
110.25, 121., 132.25, 144., 156.25, 169., 182.25, 196., 210.25, 225.,
240.25, 256., 272.25, 289., 306.25, 324., 342.25, 361., 380.25, 400.,
420.25, 441., 462.25, 484., 506.25, 529., 552.25, 576., 600.25, 625.,
650.25, 676., 702.25, 729., 756.25, 784., 812.25, 841., 870.25, 900.,
930.25, 961., 992.25}

The next step is to export the list “func” as an h5 file that XMDS2 can read. This is done by typing the following command into a Mathematica cell:

SetDirectory[NotebookDirectory[]];
Export["func.h5", {func, x1}, {"Datasets", { "function_x", "x"}}]

In the directory containing the notebook “grid.nb” you should now see the file “func.h5”. This file essentially contains the list {func, x1}. However, the hdf5 format stores func and x1 as separate entities called “Datasets”. For importation into XMDS2 it is necessary that these datasets are named. This is precisely what the segment {"Datasets", { "function_x", "x"}} in the above Mathematica command does. The dataset corresponding to the grid x1 needs to be given the name of the dimension that will be used in XMDS2 - in our case this is “x”. It does not matter what the name of the dataset corresponding to the list “func” is; in our case it is “function_x”.

The final step is to import the file “func.h5” into your primary XMDS2 file. This data will be stored as a vector called “gen_function_x”, in component “function_x”.

<vector type="real" name="gen_function_x" dimensions="x">
  <components>function_x</components>
  <initialisation kind="hdf5">
    <filename> function_x.h5 </filename>
  </initialisation>
</vector>

You’re now done. Anytime you want to use \(f(x)\) you can simply refer to “function_x” in the vector “gen_function_x”.

The situation is slightly more complicated if the function you wish to import depends on more than one dimension. For example, consider

As for the single dimensional case, we need to export an hdf5 file from XMDS2 which specifies the dimensions of the grid. As in the one dimensional case, this is done by creating a dummy vector which depends on all the relevant dimensions:

<vector type="real" name="gen_dummy" dimensions="x y">
  <components>dummy</components>
  <initialisation>
  <![CDATA[
    dummy = x;
      ]]>
  </initialisation>
</vector>

and exporting it as shown above.

After importing the grid data into Mathematica, define the multi-dimensional function which you wish to import into XMDS2:

g[x_,y_]:= x*Sin[y]

We need to create a 2x2 array of data which depends upon the imported lists x1 and y1. This can be done by using the Table function:

func := Table[g[x, p], {x, x1}, {p, p1}]

This function can be exported as an h5 file,

SetDirectory[NotebookDirectory[]];
Export["func.h5", {func, x1, y1}, {"Datasets", { "function_x", "x", "y"}}]

and imported into XMDS2 as outlined above.

Convolutions and Fourier transforms¶

When evaluating a numerical Fourier transform, XMDS2 doesn’t behave as expected. While many simulations have ranges in their spatial coordinate (here assumed to be x) that range from some negative value \(x_\text{min}\) to some positive value \(x_\text{max}\), the Fourier transform used in XMDS2 treats all spatial coordinates as starting at zero. The result of this is that a phase factor of the form \(e^{-i x_\text{min} k}\) is applied to the Fourier space functions after all forward (from real space to Fourier space) Fourier transforms, and its conjugate is applied to the Fourier space functions before all backward (from Fourier space to real space) Fourier transforms.

The standard Fourier transform is

The XMDS2 Fourier transform is

When the number of forward Fourier transforms and backwards Fourier transforms are unequal a phase factor is required. Some examples of using Fourier transforms in XMDS2 are shown below.

Example 1¶

When data is input in Fourier space and output in real space there is one backwards Fourier transform is required. Therefore the Fourier space data must be multiplied by a phase factor before the backwards Fourier transform is applied.

\[\mathcal{F}^{-1}[F(k)](x) = \tilde{\mathcal{F}}[e^{i x_\text{min} k} F(k)](x)\]

Example 2¶

Functions of the form \(h(x) = \int f(x') g(x-x') dx'\) can be evaluated using the convolution theorem:

\[\mathcal{F}[h(x)](k) = \mathcal{F}[f(x)](k) \times \mathcal{F}[g(x)](k)\]

This requires two forward Fourier transforms to get the two functions f and g into Fourier space, and one backwards Fourier transform to get the resulting product back into real space. Thus in Fourier space the product needs to be multiplied by a phase factor \(e^{-i x_\text{min} k}\)

Example 3¶

Sometimes when the convolution theorem is used one of the forward Fourier transforms is calculated analytically and input in Fourier space. In this case only one forward numerical Fourier transform and one backward numerical Fourier transform is used. The number of forward and backward transforms are equal, so no phase factor is required.

‘Loose’ geometry_matching_mode¶

Dimension aliases¶

Dimension aliases specify that two or more dimensions have exactly the same lattice, domain and transform. This can be useful in situations where the problem enforces this, for example when computing correlation functions or representing square matrices.

Dimension aliases are not just a short-hand for defining an additional dimension, they also permit dimensions to be accessed non-locally, which is essential when computing spatial correlation functions.

Here is how to compute a spatial correlation function \(g^{(1)}(x, x') = \psi^*(x) \psi(x')\) of the quantity psi:

<simulation xmds-version="2">

  <!-- name, features block -->

  <geometry>
    <propagation_dimension> t </propagation_dimension>
    <transverse_dimensions>
      <dimension name="x" lattice="1024" domain="(-1.0, 1.0)" aliases="xp" />
    </transverse_dimensions>
  </geometry>

  <vector name="wavefunction" type="complex" >
    <components> psi </components>
    <initialisation>
      <!-- initialisation code -->
    </initialisation>
  </vector>

  <computed_vector name="correlation" dimensions="x xp" type="complex" >
    <components> g1 </components>
    <evaluation>
      <dependencies> wavefunction </dependencies>
      <![CDATA[
        g1 = conj(psi(x => x)) * psi(x => xp);
      ]]>
    </evaluation>
  </computed_vector>

  <!-- integration and sampling code -->

</simulation>

In this simulation note that the vector wavefunction defaults to only having the dimension “x” even though “xp” is also a dimension (implicitly declared through the aliases attribute). vector‘s without an explicit dimensions attribute will only have the dimensions that are explicitly listed in the transverse_dimensions block, i.e. this will not include aliases.

See the example groundstate_gaussian.xmds for a complete example.

XMDS scripts look complicated! How do I start?¶

If you’re unfamiliar with XMDS2, writing a script from scratch might seem difficult. In most cases, however, the best approach is to take an existing script and modify it for your needs. At the most basic level, you can simply take a script from the /examples directory that is similar to what you want to do, change the name of the integration variable(s) and replace the line describing the differential equation to use your DE instead. That’s all you need to do, and will ensure all the syntax is correct and all the required XML blocks are present.

You can then incrementally change things such as the number of output points, what quantities get output, number of grid points, and so on. Many XMDS2 users have never written a script from scratch, and just use their previous scripts and example scripts as a scaffold when they create a script for a new problem.

Where can I get help?¶

The documentation on this website is currently incomplete, but it still covers a fair bit and is worth reading. Similarly, the example scripts in the /examples directory cover most of the functionality of XMDS2, so it’s worth looking looking through them to see if any of them do something similar to what you’re trying to do.

You should also feel free to email questions to the XMDS users’ mailing list at xmds-users@lists.sourceforge.net, where the developers and other users can assist you. You can join the mailing list by going to http://sourceforge.net/projects/xmds/ and clicking on “mailing lists.” Also, if you look through the mailing list archives, your particular problem may already have been discussed.

How should I cite XMDS2?¶

If you publish work that has involved XMDS2, please cite it as: Comput. Phys. Commun. 184, 201-208 (2013).

I think I found a bug! Where should I report it?¶

Please report bugs to the developer mailing list at xmds-devel@lists.sourceforge.net. In your email, please include a description of the problem and attach the XMDS2 script that triggers the bug.

How do I put time dependence into my vectors?¶

Standard vectors can’t have time dependence (or, more accurately, depend on the propagation_dimension variable), but computed vectors can. So, for example, if you have set your propagation_dimension as “t”, you can simply use the variable “t” in your computed vector and it will work.

Alternatively, you can explicitly use the propagation_dimension variable in your differential equation inside the <operators> block.

Can I specify the range of my domain and number of grid points at run-time?¶

Yes, you can. In your script, specify the domain and number of grid points as arguments to be passed in at run-time, use those variables in your <geometry> block rather than explicitly specifying them, and use the <validation kind="run-time" /> feature. See the Validation entry in the Reference section for an example.

While the domain can always be specified in this way, specifying the lattice size at run-time is currently only allowed with the following transforms: ‘dct’, ‘dst’, ‘dft’ and ‘none’ (see Transforms in the Reference section).

Also note that for some multi-dimensional spaces using different transforms, XMDS2 will sometimes optimise the code it generates based on the relative sizes of the dimensions. If one or more of the lattices are specified at run-time it is unable to do this and will have to make guesses. In some situations this may result in slightly slower code.

When can I use IP operators (and why should I) and when must I use EX operators?¶

An <operator> that specifies named operators to be used in integration equations can have the kind="IP" or kind="EX" attribute, standing for ‘interaction picture’ and ‘explicit’ operators respectively. Explicit operators can be used in all situations, and simply construct and calculate a new vector of the form in the square brackets. IP operators use less memory and can improve speed by allowing larger timesteps, but have two important restrictions. Use of IP operators without understanding these restrictions can lead to incorrect code.

Some explanation is in order. The IP algorithm applies the operator separately to the rest of the evolution. The reason this can be so effective is that the separate evolution can be performed exactly. The solution of the equation \(\frac{d \psi}{dt} = L \psi\) is \(\psi(t+\Delta t) = exp(L \Delta t) \psi(t)\) for arbitrarily large timestep \(\Delta t\). For a diagonal linear L, the matrix exponential is straightforward. Also, when it is constant, then the exponential can be computed and stored prior to the integration, which makes the implementation of this operator very cheap. Thus, when IP operators are defined, XMDS2 reads the equations as written by the user, and determines which operators to apply to which fields. It then implements these operators separately, and the text describing the operator inside the equations (in this example, the L[psi] term) is replaced by the numeral zero.

Therefore, the limitations of IP operators themselves means that they can only be applied to to named components of one of the integration vectors, and not functions of those components. Furthermore, an IP operator acting on a component must only be used in the derivative for that particular component. Secondly, due to the implementation of IP operators in XMDS2, it is not safe to use them in comments, or in conjunction with declared variables. It is also not safe to multiply or divide them by any factors, functions or vectors. They must turn up in a purely additive way when defining the derivative of a component of an integration vector. The XMDS2 parser attempts to catch possible violations of these rules, and will produce warnings in some cases.

Visual Editors¶

In this section goes stuff about how to set up TextMate (or other editors to highlight xpdeint scripts).

Geometry and transform-based tricks¶

Reduce code complexity¶

Avoid transcendental functions like \(\sin(x)\) or \(\exp(x)\) in inner loops. Not all operations are made equal, use multiplication over division.

<operator kind="ip">
  <operator_names>T</operator_names>
  <![CDATA[
    T = -i*0.5*hbar/M*(kx*kx + ky*ky);
  ]]>
</operator>
<![CDATA[
  dpsi_dt = T[psi] + /* other terms */
]]>

<operator kind="ip" dimensions="x">
  <operator_names>Tx</operator_names>
  <![CDATA[
    Tx = -i*0.5*hbar/M*kx*kx;
  ]]>
</operator>
<operator kind="ip" dimensions="y">
  <operator_names>Ty</operator_names>
  <![CDATA[
    Ty = -i*0.5*hbar/M*ky*ky;
  ]]>
</operator>
<![CDATA[
  dpsi_dt = Tx[psi] + Ty[psi] + /* other terms */
]]>

<operators>
  <integration_vectors basis="kx">wavefunction</integration_vectors>
  <operator kind="ip" type="imaginary" >
    <operator_names>Lxx</operator_names>
    <![CDATA[
      Lxx = -i*0.5*hbar_M*(kx*kx);
    ]]>
  </operator>
  <![CDATA[

    dpsi0_dt = Lxx[psi0] - i*Omega*psi1;
    dpsi1_dt = Lxx[psi1] - i*Omega*psi0;

  ]]>
</operators>