% $Id: moreExamples.tex 1641 2008-01-10 02:55:05Z joehope $

% Copyright (C) 2000-2007
%
% Code contributed by Greg Collecutt, Joseph Hope and the xmds-devel team
%
% This file is part of xmds.
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software 
% Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.

\chapter{More Examples} 
\label{chap:moreExamples}

\section{ndparamp.xmds} 
\label{sec:ndparamp}

\begin{xmdsCode}
<?xml version="1.0"?>
<!--Non-Degenerate Parametric Amplifier-->
<!--Simulton formation for logical switching-->

<simulation>

  <name>ndparamp</name>
  <prop_dim>z</prop_dim>
  <error_check>yes</error_check>

  <globals>
    <![CDATA[
      const double e1 =350;
      const double e2 =350;
      const double r1 = 1;
      const double r2 = 1;
      const double vy1 =  0.5;
      const double vy2 = -0.5;
      const double yc1 = -0.2;
      const double yc2 =  0.2;
      const double tc1 = 0;
      const double tc2 = 0;

      double amp1=sqrt(e1/2/M_PI/r1/r1);
      double amp2=sqrt(e2/2/M_PI/r2/r2);
    ]]>
  </globals>

  <field>
    <name>main</name>
    <dimensions>   y       t     </dimensions>
    <lattice>     100     100    </lattice>
    <domains>  (-10,10) (-10,10) </domains>
    <samples>1 1</samples>

    <vector>
      <name>main</name>
      <type>complex</type>
      <components>ff1 ff2 sh</components>
      <fourier_space>no no</fourier_space>
      <![CDATA[
        ff1 = pcomplex(amp1*exp(-pow((y - yc1)/r1/2,2)
                        -pow((t - tc1)/r1/2,2)),+vy1*y);
        ff2 = pcomplex(amp2*exp(-pow((y - yc2)/r2/2,2)
                        -pow((t - tc2)/r2/2,2)),+vy2*y);
        sh  = rcomplex(0,0);
      ]]>
    </vector>

    <vector>
      <name>vc1</name>
      <type>double</type>
      <components>damping</components>
      <fourier_space>no no</fourier_space>
      <![CDATA[
        damping=1.0*(1-exp(-pow((y*y + t*t)/8/8,10)));
      ]]>
    </vector>
  </field>

  <sequence>
    <integrate>
      <algorithm>RK4IP</algorithm>
      <interval>10</interval>
      <lattice>500</lattice>
      <samples>50 50</samples>
      <k_operators>
        <constant>yes</constant>
        <operator_names>Lap1 Lap2</operator_names>
        <![CDATA[
          Lap1 = i*(-(ky*ky + kt*kt) - 1);
          Lap2 = i*(-(ky*ky + kt*kt)/2 - 1);
        ]]>
      </k_operators>
      <vectors>main vc1</vectors>
        <![CDATA[
          dff1_dz = Lap1[ff1] + i*~ff2*sh  - damping*ff1;
          dff2_dz = Lap1[ff2] + i*~ff1*sh  - damping*ff2;
          dsh_dz  = Lap2[sh]  + i* ff1*ff2 - damping* sh;
        ]]>
    </integrate>
  </sequence>

  <output>

    <group>
      <sampling>
        <fourier_space> no  no </fourier_space>
        <lattice>       50  0  </lattice>
        <moments>pow_dens</moments>
        <![CDATA[
          pow_dens = ~ff1*ff1 + ~ff2*ff2 + 2*~sh*sh;
        ]]>
      </sampling>
    </group>

    <group>
      <sampling>
        <fourier_space> no  no </fourier_space>
        <lattice>       0  0  </lattice>
        <moments>etot</moments>
        <![CDATA[
          etot = ~ff1*ff1 + ~ff2*ff2 + 2*~sh*sh;
        ]]>
      </sampling>
    </group>
  </output>
</simulation>
\end{xmdsCode}

This simulation describes how to solve for the evolution of three
bosonic fields governed by a non-degenerate parametric interaction, as
described in Equations~(\ref{eq:TIIa}) and~(\ref{eq:TIIb}). Perhaps
not surprisingly, the majority of the above script is to generate the
particular initialisation conditions and to define a number of output
moments---the portion concerned with actually implementing these
equations is not that great.
\begin{equation}
  \frac{\partial \phi _{j}}{\partial \xi } =i\left[ \left(
    \frac{\partial ^{2}}{\partial \tau ^{2}} + \frac{\partial
      ^{2}}{\partial \zeta ^{2}} + i\Gamma (\tau ,\zeta )-1\right) \phi
    _{j}+\phi ^{*}_{3-j} \phi _{3}\right], \quad j=1,2;
  \label{eq:TIIa}
\end{equation}
\begin{equation}
  \frac{\partial \phi _{3}}{\partial \xi } = i\left[ \left(
    \frac{1}{\sigma}\frac{\partial ^{2}}{\partial \tau ^{2}} +
    \frac{1}{2}\frac{\partial ^{2}}{\partial \zeta ^{2}} +i\Gamma (\tau
    ,\zeta )- \gamma \right) \phi _{3}+\phi _{1}\phi _{2}\right].
  \label{eq:TIIb}
\end{equation}

The main difference between this simulation and the \ttt{nlse.xmds}
simulation is that it now has a three component field which has two
transverse dimensions. Also two moment groups are being evaluated, one
being integrated over the ``\ttt{t}'' dimension, and the other
integrated over both transverse dimensions.
\begin{figure}[ht]
\centerline{%
\begin{tabular}{cc}
\includegraphics[width=\figwidth]{figures/ndparamp1} &
\includegraphics[width=\figwidth]{figures/ndparamp2}\\ a) Time
integrated intensity & b) Total energy (vs z)
\end{tabular}%
}
\caption{Results for ndparamp.xmds}
\label{fig:ndparamp}
\end{figure}

\section{kubo.xmds} 
\label{sec:kubo.xmds}

\begin{xmdsCode}
<?xml version="1.0"?>
<!--Example Kubo oscillator simulation-->

<simulation>

  <name>kubo</name>
  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <stochastic>yes</stochastic>
  <paths>1</paths>
  <use_mpi>no</use_mpi>
  <seed>1 2</seed>
  <noises>1</noises>
   
  <field>
    <samples>1</samples>
    <vector>
      <name>main</name>
      <type>complex</type>
      <components>z</components>
      <![CDATA[
        z = 1;
      ]]>
    </vector>
  </field>

  <sequence>
    <integrate>
      <algorithm>SIEX</algorithm>
      <interval>10</interval>
      <lattice>1000</lattice>
      <samples>100</samples>
      <iterations>3</iterations>
      <![CDATA[
        dz_dt = i*z*n_1;
      ]]>
    </integrate>

  </sequence>

  <output>
    <group>
      <sampling>
      <moments>realz</moments>
      <![CDATA[
	realz = z;
      ]]>
      </sampling>
    </group>
  </output>
</simulation>
\end{xmdsCode}

The kubo oscillator is described in \eqn{eq:kubo}, in which
the argument of the complex vector {\em z} is ``blown'' about by a
(real) Gaussian noise term, $\xi(t)$. This is a simple stochastic ODE.
\begin{equation}
\frac{\partial z}{\partial t } = i \xi(t) z.
\label{eq:kubo}
\end{equation}
Such Gaussian noise terms, in analytic form, are correlated in time
and space through Dirac delta functions, as shown in \eqn{eq:noiseDirac}.
\begin{equation}
\left< \xi_{i}(\vect{x})\xi_{j}(\vect{x'}) \right> = \delta_{i,j}
\Pi_{i=0}^{N} \delta(x^i-x'^i).
\label{eq:noiseDirac}
\end{equation}
However, when solving stochastic DEs numerically, algorithms work with
discrete time intervals and lattice spacings. Therefore these Dirac
delta correlations must be transformed to Kronecker delta correlations
using the integration time step and the spatial volume of the lattice,
as shown in \eqn{eq:noiseKronecker}:
\begin{equation}
\left< \xi_{i}(\vect{x})\xi_{j}(\vect{x'}) \right>
  = \frac{\delta_{i,j} \Pi_{i=0}^{N} \delta_{x^i,x'^i}}{\Pi_{i=0}^{N}
  \Delta x^i}.
\label{eq:noiseKronecker}
\end{equation}

The good news is that \xmds calculates this for the user---all that
has to be done is to specify, as a simple \xmdsTag{noises} assignment
within the \xmdsTag{simulation} element, the maximum number of noise
terms required in any one segment, and then reference them as
\ttt{n\_1}, \ttt{n\_2}, etc. as can be seen in this example. These noises are
available within the initialisation code for each field
\xmdsTag{vector}, within the main integration equations code (not in
the \xmdsTag{k\_operators} code), and in the code for any
\xmdsTag{filter} segments. Within the initialisation code and
\xmdsTag{filter} code the variances of the noises are determined by the
lattice cell volume product for the particular \xmdsTag{fourier\_space}
specification as shown in \eqn{eq:var}:
\begin{equation}
\left< n_{i}n_{j} \right> = \frac{\delta_{i,j}}{\Pi \Delta k^{m} \Pi
\Delta x^{n} },
\label{eq:var}
\end{equation}
where $m$ are the transverse dimensions in Fourier space and $n$ the
transverse dimensions in normal space. Also note that
\begin{equation}
\Delta k^i = \frac{2\pi}{x^i_{max} - x^i_{min}}.
\end{equation}

Within the main integration equations the variances must also reflect
the integration step size, as given by \eqn{eq:var_x-space}:
\begin{equation}
\left< n_{i}n_{j} \right> = \frac{\delta_{i,j}}{\Pi_{i=0}^{N} \Delta
x^i}.
\label{eq:var_x-space}
\end{equation}
\xmds uses the Box-Mueller technique, as shown in
\eqn{eq:box-mueller}, which generates a pair of Gaussian noises,
$\xi_1$ and $\xi_2$, from a pair of random numbers, $x_1$ and $x_2$,
that have a uniform distribution between zero and one.
\begin{equation}
  \xi_1 + i\xi_2 = \left[-2\Delta ln(x_1)\right]^\frac{1}{2} e^{i2\pi
    x_2} \; ; \; \mathcal P(x_i<y) = y \; ; \; y: 0 < y \leq 1;y\epsilon
  \mathcal R.
\label{eq:box-mueller}
\end{equation}

Estimating the error between full and half step sizes now poses an
interesting problem, since both evolutions must use the same
underlying noise (which is a function of both space and time) if the
difference between these paths is to be meaningful. The random number
generator must be reset before each of these integrations, but how is
the noise in the half-step case appropriated?  There are two
methods. The first is to use the same noise for both half-steps as is
used for the one whole step. This is undesirable since it makes sense
to use the half-step integration results (being the more accurate) for
the final output, and therefore independent noises must be used for
each step. So, the second solution is to do just this, and use the
average of the two noises when calculating the full-step
integration. Now, suppose the problem uses $N$ noises and has a
transverse lattice of $M$ points.  Within the SIEX, RK4IP, and the RK4EX
integration algorithms the main field vector is swept through severals
times in the course of each time step, thus a $M \times N$ vector of
noise must be calculated at the beginning of the time step and
referenced during the calculation of the main field vector's
derivatives. In the full-step case two such vectors are calculated,
and then averaged to provide the equivalent full-step noise. However,
the SIIP algorithm only sweeps though the main field vector once for
each time step, and so only $N$ noises need be calculated at a time
(thus saving on memory and RAM access), provided {\em two} independent
random number generators are used: the first generator is used for the
first half step, the second for the second half step, and the average
of both is used for the full step. In C this is done with the {\em
erand48(n)} function which uses the 48-bit integer $n$ to generate the
next random number, advancing $n$ to the next in sequence in the
process. The states of the two independent generators are simply two
independent integers, $n_1$ and $n_2$. The user supplies the initial
values for these integers in the \xmdsTag{seed} assignment.

Stochastic problems are very well suited to parallel computer
architectures, as different paths can run on different processors, and
do not have to transfer information until the integration is
complete. Multiple path stochastic problems such as this may be
parallised using MPI routines. If an MPI compiler was specified in the
configuration step of installation (this will often be \ttt{mpicc}),
then all that needs to be done is to toggle the optional
\xmdsTag{use\_mpi} assignment to \ttt{yes}. \xmds will then place the
appropriate MPI calls in the output code and compile it with the MPI
compiler. The executable should then be run through the MPI execution
handler (probably \ttt{mpirun}), with the number of processors option
supplied. For example if 16 processors are available then the final
command for execution would be
\begin{shellCode}
% mpirun -np 16 kubo
\end{shellCode}
Note that whole paths are assigned independently to the processors, so
there is no benefit in specifying more processors than there are paths
in the simulation. It is not necessary for the number of processors to
be a factor of the number of paths, some processors will simply do one
more paths than others.

Also note that for both of the semi-implicit algorithms, SIEX and
SIIP, an \xmdsTag{iterations} assignment may be used within the
\xmdsTag{integrate} element to specify the number of iterations to use
in the method (refer Sections~\ref{sec:siMethod},
\ref{sec:siexMethod}, and~\ref{sec:siipMethod}). This assignment is
optional, and will default to three when absent.

The reason that the kubo oscillator is used as an example is that it
has an analytic solution, as shown in \eqn{eq:kuboSolution}.
\fig{fig:kubo} shows the results for a single trajectory and an
averaged trajectory, illustrating the expected behaviour.
\begin{equation} 
\left< z(t) \right> = z_0 e^{-\frac{t}{2}}.
\label{eq:kuboSolution}
\end{equation}

\begin{figure}[ht]
\centerline{%
\begin{tabular}{cc}
\includegraphics[width=\figwidth]{figures/kubo1} &
\includegraphics[width=\figwidth]{figures/kubo2}\\ a) Single path & b)
1024 path mean
\end{tabular}%
}
\caption{Results for kubo.xmds}
\label{fig:kubo}
\end{figure}

Since the variance of the noise terms scale with the inverse of the
integration step size, the integration method suffers a loss of order
with regard to error vs step size. While this is normally a second
order method for non-stochastic problems, it becomes a first order
method for problems with noise, as was explained in
\Sec{sec:eulerMethods}.

\section{fibre.xmds} 
\label{sec:fibre.xmds}

\begin{xmdsCode}
<?xml version="1.0"?>
<!--Example fibre noise simulation-->

<simulation>

  <name>fibre</name>
  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <stochastic>yes</stochastic>
  <use_mpi>no</use_mpi>
  <paths>1</paths>
  <seed>1 2</seed>
  <noises>2</noises>

  <globals>
<![CDATA[
const double ggamma = 1;
const double beta = sqrt(2*2*M_PI*ggamma/10);
]]>
  </globals>

  <field>
    <name>main</name>
    <dimensions>   x   </dimensions>
    <lattice>     50  </lattice>
    <domains>  (-5,5)  </domains>
    <samples>1</samples>

    <vector>
      <name>main</name>
      <type>complex</type>
      <components>phi</components>
      <fourier_space>no</fourier_space>
<![CDATA[
phi=0;
]]>
    </vector>
  </field>

  <sequence>
    <integrate>
      <algorithm>SIIP</algorithm>
      <interval>2.5</interval>
      <lattice>5000</lattice>
      <samples>50</samples>
      <k_operators>
        <constant>yes</constant>
        <operator_names>L</operator_names>
<![CDATA[
L = i*(-kx*kx);
]]>
      </k_operators>
      <iterations>3</iterations>
<![CDATA[
dphi_dt = L[phi] - ggamma*phi + beta/sqrt(2)*complex(n_1,n_2);
]]>
    </integrate>
  </sequence>

  <output>
    <group>
      <sampling>
      <fourier_space> yes </fourier_space>
      <lattice>       50  </lattice>
      <moments>pow_dens</moments>
<![CDATA[
pow_dens = conj(phi)*phi;
]]>
      </sampling>
    </group>
  </output>
</simulation>
\end{xmdsCode}

This simulation solves \eqn{eq:fibre}, in which a one dimensional
damped field is subject to a complex noise. This is a stochastic PDE.
\begin{equation}
\frac{\partial \psi}{\partial t } = -i \frac{\partial^{2}
 \psi}{\partial x^{2}} - \gamma \psi + \frac{\beta}{\sqrt{2}}
 (\xi_1(x,t) + i\xi_2(x,t)).
\label{eq:fibre}
\end{equation}

Again the reason for using this as an example of a stochastic PDE is
that it has an analytic solution, as shown in \eqn{eq:fibreSolution}.
\fig{fig:fibre} displays the results of this simulation in Fourier
space for a single trajectory and an averaged trajectory, which appear
as expected.
\begin{equation} 
\left< | \psi (k,t) |^{2} \right> = e^{-2 \gamma t} |\psi_0 (k)|^{2} +
\frac{\beta^{2}L_{x}}{4 \pi \gamma} (1-e^{-2 \gamma t}),
\label{eq:fibreSolution}
\end{equation}
where $L_{x}$ is the length of the $x$ domain.
\begin{figure}[ht]
\centerline{%
\begin{tabular}{cc}
\includegraphics[width=\figwidth]{figures/fibre1} &
\includegraphics[width=\figwidth]{figures/fibre2}\\ a) Single path &
b) 1024 path mean
\end{tabular}
}%
\caption{Results for fibre.xmds}
\label{fig:fibre}
\end{figure}

One important issue here is that the variance of the noise terms now
scales with the product of the number of lattice points (for a given
domain). Hence changing to a finer lattice actually increases the
single trajectory error, and the relationship between the error and
the lattice product will depend on the order of the spatial
derivatives. The only way to overcome this is to reduce the
integration step size, which added to the fact that there are more
lattice points in the first place, means that fine lattice resolution
in a stochastic PDE is computationally {\em very} expensive.

\section{tla.xmds} 
\label{sec:tla.xmds}

\begin{xmdsCode}
<?xml version="1.0"?>
<!--Two Level Atom Example simulation to illustrate a
cross propagating field-->

<simulation>

  <prop_dim> z </prop_dim>

  <globals>
    <![CDATA[
      const double g  = 1;
      const double t0 = 1;
    ]]>
  </globals>

  <field>
    <dimensions>   t    </dimensions>
    <lattice>     100   </lattice>
    <domains> (-10, 15) </domains>
    <samples> 1 0 </samples>

    <vector>
      <name> main </name>
      <type> double </type>
      <components> E </components>
      <![CDATA[
        E = 2/t0/cosh(t/t0);
      ]]>
    </vector>

    <vector>
      <name> cross </name>
      <type> double </type>
      <components> P N </components>
      <![CDATA[
          P =  0;
          N = -1;
      ]]>
    </vector>
  </field>

  <sequence>
    <integrate>
      <algorithm> RK4EX </algorithm>
      <interval> 4 </interval>
      <lattice> 50 </lattice>
      <samples> 50 50 </samples>
      <vectors> main cross </vectors>
        <![CDATA[
          dE_dz =  g*P;
        ]]>
      <cross_propagation>
        <vectors> cross </vectors>
        <prop_dim> t </prop_dim>
        <![CDATA[
          dP_dt =  E*N;
          dN_dt = -E*P;
        ]]>
      </cross_propagation>
    </integrate>
  </sequence>

  <output>

    <group>
      <sampling>
        <lattice> 50 </lattice>
        <moments> pow_dens </moments>
        <![CDATA[
          pow_dens = E*E;
        ]]>
      </sampling>
    </group>

    <group>
      <sampling>
        <vectors> main cross </vectors>
        <lattice> 50 </lattice>
        <moments> P_out N_out </moments>
        <![CDATA[
          P_out = P;
          N_out = N;
        ]]>
      </sampling>
      </group>
  </output>
</simulation>
\end{xmdsCode}

This simulation solves for the propagation of an optical pulse through
a field of atoms having a transition frequency tuned to that of the
optical pulse centre frequency. The atoms are modelled as ``two
level'' atoms. The propagation equations, shown in \eqn{eq:tla}, are
deceptively simple.
\begin{align}
\frac{\partial E(t,z)}{\partial z} & = gP,\nonumber \\
\frac{\partial P(t,z)}{\partial t} & = EN,\nonumber \\
\frac{\partial N(t,z)}{\partial t} & = -EP.
\label{eq:tla}
\end{align}

The reality of this problem is that there are three components, two
propagating in the main propagation dimension, $t$, and the other in
the transverse dimension $z$. The component $E$ is the electric field
amplitude, $P$ the polarisation state of the atoms, and $N$ the
excitation state of the atoms (-1 being all in the ground state and +1
being all in the excited state). Lastly $g$ is a coupling constant
between the electric field and the atoms.

The curious feature of this set of PDEs, in fact the very reason why
it is chosen it as an example, is that there exists a soliton solution
for the electric field, as shown in \eqn{eq:tlaSolution}:
\begin{equation} 
E(t,z) = \frac{2}{gt_0} \mathrm{sech} \left(\frac{t-az}{t_0}\right),
\label{eq:tlaSolution}
\end{equation}
which is time lagged with propagation at the rate
\begin{equation}
a = \frac{g}{(t_0)^{2}}.
\end{equation}

The result for the electric field in the above simulation is shown in
\fig{fig:tla}, in which the soliton solution is evident.
\begin{figure}[ht]
\centerline{\includegraphics[width=\figwidth]{figures/tla}}
\caption{Results for tla.xmds}
\label{fig:tla}
\end{figure}

The cross propagating components, $N$ and $P$, may still be thought of
as existing in the same space as the main vector component, $E$, and
so they are declared as an extra vector in the \xmdsTag{field}
element. However, the equations governing the evolution of such cross
vectors are not allowed to include any transverse derivatives --
i.e. they are not allowed to be PDEs. Therefore, through the main
vector equations may be PDEs, the cross propagating vector need not be
transformed to Fourier space when the main vector is. So although the
cross vector components could be included as part of the main vector,
they are better defined as a separate vector for efficiency reasons.

In the SIIP integration algorithm the transverse evolution of the
cross vector is calculated simultaneously with the forward evolution
of the main vector, but in all other algorithms the cross vector is
calculated prior to calculating the main vector derivatives. Thus the
governing equations for the cross vector must be separated from those
for the main vector. This is done by including a
\xmdsTag{cross\_propagation} element within the main \xmdsTag{integrate}
element, and placing the cross vector equations within. Also required
within this element is a list of the \xmdsTag{vectors} that are to be
cross propagated, and a \xmdsTag{prop\_dim} assignment specifying the
dimension of cross propagation. The \xmdsTag{vectors} that were made
accessible for the main equations will also be accessible here.



\section{highdim.xmds} 
\label{sec:highdim}

\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>
  <name>highdim</name>

  <!-- Global system parameters and functionality -->
  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <use_mpi>yes</use_mpi> 
  <use_wisdom>yes</use_wisdom>
  <benchmark>yes</benchmark>
  
  <!-- Global variables for the simulation -->
  <globals>
  <![CDATA[
    const double noise = 0.0;
    const double hbar = 1.05500000000e-34;
    const double M = 1.409539200000000e-25;
    const double omegax = 0.58976353090742;
    const double omegay = 0.58976353090742;
    const double omegaz = 0.58976353090742/30;
    const double U11 = 2.974797272874263e-51;
    const double U13 = -1.417820412490823e-50;
    const double U33 = 2.974797272874263e-51;
    const double inum = 1.0e6;
    const double Uoh11 = U11/hbar;
    const double Uoh13 = U13/hbar;
    const double Uoh33 = U33/hbar;
    const double mu = pow(15*inum*U11*omegax*omegay
              *omegaz/M_PI/4,0.4)*pow(M,0.6)/2;
    const double delta = 1.0e9;
    const double F = 2.0e-2;
    const double g = sqrt(Uoh11*2.0/delta);
    const double loss11=1.0e-2;
    const double loss12=1.6e-22;
    const double loss31=1.0e-2;
    const double loss32=1.6e-22;
    const double loss132=8.0e-17;
    const double chi = F*g*delta;
    const double biggamma = g*g*delta/2;
    const double gam13 = Uoh13/chi;
    const double gam33 = Uoh33/chi;
    const double gameff = (Uoh11-biggamma)/chi;
    const double gamloss11=loss11/2/chi;
    const double gamloss12=loss12/chi;
    const double gamloss31=loss31/2/chi;
    const double gamloss32=loss32/chi;
    const double gamloss132=loss132/chi;
    const double cnoise =  noise/sqrt(2.0); 
  ]]>
  </globals>  

    <argv>
	   <arg>
	       <name>kjoek</name>
	       <type>double</type>
	       <default_value>-1.0e6</default_value>
	   </arg>
	   <arg>
	       <name>joekappamax</name>
	       <type>double</type>
	       <default_value>1.0e2</default_value>
	   </arg>
    </argv>

  <!-- Field to be integrated over -->
  <field>
    <dimensions>x y z</dimensions>
    <lattice>16 16 16</lattice>
    <domains>(-1.2e-4,1.2e-4) (-1.2e-4,1.2e-4) (-8.0e-3,8.0e-3)</domains>
    <samples>1 1 1</samples>
    
    <vector>
      <name> vc1 </name>
      <type>double</type>
      <components>vcore V1r V3r gV1r gV3r</components>
      <fourier_space>no no no</fourier_space>
      <![CDATA[
   vcore = (omegax*omegax*x*x+omegay*omegay*y*y+omegaz*omegaz*z*z);
   V1r = 0.5*M*vcore/hbar/chi -(gameff+gam13/2)/2/(dx*dy*dz);
   V3r = M*vcore/hbar/chi -(gam13/2+gam33)/2/(dx*dy*dz);
   gV1r = 0.5*M*vcore/hbar/chi;
   gV3r = M*vcore/hbar/chi;
      ]]>
    </vector>
    
    <vector>
      <name> main </name>
      <type>complex</type>
      <components>phi1a phi1b phi3a phi3b gphi1a gphi3a</components>
      <fourier_space>no no no</fourier_space>
      <vectors> vc1 </vectors>
      <![CDATA[
        const double realfn = (mu-0.5*M*vcore)/Uoh11/hbar;
        
        phi1a = realfn>0. ? complex(sqrt(realfn),0) : complex(0,0);
        phi1b = realfn>0. ? complex(sqrt(realfn),0) : complex(0,0);
        phi3a = complex(0,0);
        phi3b = complex(0,0);
        gphi1a = realfn>0. ? complex(sqrt(realfn),0) : complex(0,0);
        gphi3a = complex(0,0);
      ]]>
    </vector>
  </field>

  <!-- The sequence of integrations to perform -->
  <sequence>
    <integrate>
      <algorithm>ARK89IP</algorithm>
      <interval>1e-7</interval>
       <tolerance>1.0e-7</tolerance>
      <lattice>1000</lattice>
      <samples>10 10 1</samples>
      <k_operators>
        <constant>yes</constant>
        <operator_names> L2p L2n L4p L4n </operator_names>
        <![CDATA[
          L2p = complex(0,-hbar/M/2/chi*(kx*kx+ky*ky+kz*kz));
          L2n = complex(0, hbar/M/2/chi*(kx*kx+ky*ky+kz*kz));
          L4p = complex(0,-hbar/M/4/chi*(kx*kx+ky*ky+kz*kz));
          L4n = complex(0, hbar/M/4/chi*(kx*kx+ky*ky+kz*kz));
        ]]>
      </k_operators>
      
        <moment_group>
        <moments>chippy</moments>
        <integrate_dimension>yes yes yes</integrate_dimension>
	    <![CDATA[
		chippy += ~gphi1a*gphi1a;   
		]]>
        </moment_group>    
          
        <moment_group>
        <moments>ippy ichippy</moments>
        <integrate_dimension>no no no</integrate_dimension>
	    <![CDATA[
		ippy += phi1a;
		ichippy += gphi1a;   
		]]>
        </moment_group>      
        <moment_group>
        <moments>py ic</moments>
        <integrate_dimension>yes yes no</integrate_dimension>
	    <![CDATA[
		py += phi1a;
		ic += gphi1a;   
		]]>
        </moment_group>      
        <moment_group>
        <moments>ppy ichi</moments>
        <integrate_dimension>no yes yes</integrate_dimension>
	    <![CDATA[
		ppy += phi1a;
		ichi += gphi1a;   
		]]>
        </moment_group>      
           
      <vectors> main vc1 </vectors>
      <![CDATA[
        const complex dens1 = phi1b*phi1a;
        const complex dens3 = phi3b*phi3a;
        
    const double gdens1 = (gphi1a.re*gphi1a.re+gphi1a.im*gphi1a.im);
    const double gdens3 = (gphi3a.re*gphi3a.re+gphi3a.im*gphi3a.im);

        dphi1a_dt = L2p[phi1a] + (-i*V1r-gamloss11
          +(gamloss132/2+gamloss12)/2/(dx*dy*dz))*phi1a 
          + (-i*gameff-gamloss12)*dens1*phi1a 
          - (i*gam13+gamloss132)*dens3*phi1a -i*phi1b*phi3a 
          + i*chippy*ippy;
                               
        dphi1b_dt = L2n[phi1b] + (i*V1r-gamloss11
          +(gamloss132/2+gamloss12)/2/(dx*dy*dz))*phi1b 
          + (i*gameff-gamloss12)*dens1*phi1b 
          + (i*gam13-gamloss132)*dens3*phi1b +i*phi1a*phi3b;
                             
        dphi3a_dt = L4p[phi3a] + (-i*V3r-gamloss31
          +(gamloss132/2+gamloss32)/2/(dx*dy*dz))*phi3a 
          + (-i*gam33-gamloss32)*dens3*phi3a 
          - i*0.5*phi1a*phi1a -(i*gam13+gamloss132)*dens1*phi3a;

        dphi3b_dt = L4n[phi3b] + (i*V3r-gamloss31
          +(gamloss132/2+gamloss32)/2/(dx*dy*dz))*phi3b 
          + (i*gam33-gamloss32)*dens3*phi3b 
          + i*0.5*phi1b*phi1b +(i*gam13-gamloss132)*dens1*phi3b;
                               
        dgphi1a_dt = L2p[gphi1a] + (-i*gV1r-gamloss11)*gphi1a 
           +(-i*gameff-gamloss12)*gdens1*ichippy 
          - (i*gam13+gamloss132)*gdens3*gphi1a-i*conj(gphi1a)*gphi3a;
                                 
        dgphi3a_dt = L4p[gphi3a] + (-i*gV3r-gamloss31)*gphi3a 
          +(-i*gam33-gamloss32)*gdens3*gphi3a 
          - i*0.5*gphi1a*gphi1a +(i*gam13-gamloss132)*gdens1*gphi3a;
      ]]>
    </integrate>
  </sequence>

  <!-- The output to generate -->
  <output format="binary" precision="double">
    <group>
      <sampling>
        <fourier_space>   no    no no</fourier_space>
        <lattice>         16    1  1</lattice>
        <moments>atoms molecules gatoms gmolecules</moments>
        <![CDATA[
          atoms=phi1b*phi1a;
          molecules=phi3b*phi3a;
          gatoms=conj(gphi1a)*gphi1a;
          gmolecules=conj(gphi3a)*gphi3a;
        ]]>
      </sampling>
    </group>
    <group>
      <sampling>
        <fourier_space>    no   no  no</fourier_space>
        <lattice>           0    16   0</lattice>
        <moments>rn_1 rn_2 grn_1 grn_2 excitedn</moments>
        <![CDATA[
    rn_1 = phi1b*phi1a;
    rn_2 = phi3b*phi3a;
    grn_1 = conj(gphi1a)*gphi1a;
    grn_2 = conj(gphi3a)*gphi3a;
    excitedn = g*g/4*phi1b*phi1b*phi1a*phi1a+F*F*phi3b*phi3a
                     - F*g/2*(phi1b*phi1b*phi3a+phi1a*phi1a*phi3b);
        ]]>
      </sampling>
    </group>
    <group>
      <sampling>
        <fourier_space>   no no no</fourier_space>
        <lattice>         4    8  16</lattice>
        <moments>atomsr moleculesr atomsi moleculesi</moments>
        <![CDATA[
          atomsr=phi1a;
          moleculesr=phi3a;
          atomsi=-i*gphi1a;
          moleculesi=-i*gphi3a;
        ]]>
      </sampling>
    </group>
  </output>

</simulation>
\end{xmdsCode}

This simulation is included to highlight the usage of moment groups in evolution.
Here we wish to use various variables integrated over one, two or all the 
transverse dimensions.  This is done in integrate or filter elements by the inclusion of 
a moment group.  Using the first as an example:

\begin{verbatim}        <moment_group>
        <moments>chippy</moments>
        <integrate_dimension>yes yes yes</integrate_dimension>
	    <![CDATA[
		chippy += ~gphi1a*gphi1a;   
		]]>
        </moment_group>         \end{verbatim}  

The syntax is similar to output and filter syntax, but note that the equality must be a "+=" and not simply "=".  This element provides the integral of the modulus squared of the field gphi1a, which can then be used normally in the integration code by the designated name, chippy.  Other variables are actually fields in which any number of transverse dimensions may be integrated, and the others are left.

When \xmdsTag{functions} and \xmdsTag{moment\_group} elements are used, the position of the \xmdsTag{vectors} tag is crucial.  It specifies when the integrate code is to be executed, which will usually need to be after the moment groups are calculated.  


\section{highdim\_vector\_version.xmds} 
\label{sec:highdimVectorVersion}

\begin{xmdsCode}
<?xml version="1.0"?>
<simulation>
  <name>highdim</name>
  <!-- Global system parameters and functionality -->
  <prop_dim>t</prop_dim>
  <error_check>yes</error_check>
  <use_mpi>yes</use_mpi> 
  <use_wisdom>yes</use_wisdom>
  <benchmark>yes</benchmark>
  
  <!-- Global variables for the simulation -->
  <globals>
  <![CDATA[
    const double noise = 0.0;
    const double hbar = 1.05500000000e-34;
    const double M = 1.409539200000000e-25;
    const double omegax = 0.58976353090742;
    const double omegay = 0.58976353090742;
    const double omegaz = 0.58976353090742/30;
    const double U11 = 2.974797272874263e-51;
    const double U13 = -1.417820412490823e-50;
    const double U33 = 2.974797272874263e-51;
    const double inum = 1.0e6;
    const double Uoh11 = U11/hbar;
    const double Uoh13 = U13/hbar;
    const double Uoh33 = U33/hbar;
    const double mu = pow(15*inum*U11*omegax*omegay
              *omegaz/M_PI/4,0.4)*pow(M,0.6)/2;
    const double delta = 1.0e9;
    const double F = 2.0e-2;
    const double g = sqrt(Uoh11*2.0/delta);
    const double loss11=1.0e-2;
    const double loss12=1.6e-22;
    const double loss31=1.0e-2;
    const double loss32=1.6e-22;
    const double loss132=8.0e-17;
    const double chi = F*g*delta;
    const double biggamma = g*g*delta/2;
    const double gam13 = Uoh13/chi;
    const double gam33 = Uoh33/chi;
    const double gameff = (Uoh11-biggamma)/chi;
    const double gamloss11=loss11/2/chi;
    const double gamloss12=loss12/chi;
    const double gamloss31=loss31/2/chi;
    const double gamloss32=loss32/chi;
    const double gamloss132=loss132/chi;
    const double cnoise =  noise/sqrt(2.0); 
  ]]>
  </globals>  

    <argv>
	   <arg>
	       <name>kjoek</name>
	       <type>double</type>
	       <default_value>-1.0e6</default_value>
	   </arg>
	   <arg>
	       <name>joekappamax</name>
	       <type>double</type>
	       <default_value>1.0e2</default_value>
	   </arg>
    </argv>

  <!-- Field to be integrated over -->
  <field>
    <dimensions>x y z</dimensions>
    <lattice>16 16 16</lattice>
<domains>(-1.2e-4,1.2e-4) (-1.2e-4,1.2e-4) (-8.0e-3,8.0e-3)</domains>
    <samples>1 1 1</samples>
    
    <vector>
      <name> vc1 </name>
      <type>double</type>
      <components>Vr(5)</components>
      <fourier_space>no no no</fourier_space>
      <![CDATA[
     Vr(1) = (omegax*omegax*x*x+omegay*omegay*y*y+omegaz*omegaz*z*z);
     Vr(2) = 0.5*M*Vr(1)/hbar/chi -(gameff+gam13/2)/2/(dx*dy*dz);
     Vr(3) = M*Vr(1)/hbar/chi -(gam13/2+gam33)/2/(dx*dy*dz);
     Vr(4) = 0.5*M*Vr(1)/hbar/chi;
     Vr(5) = M*Vr(1)/hbar/chi;
      ]]>
    </vector>
    
    <vector>
      <name> main </name>
      <type>complex</type>
      <components> phi(6) </components>
      <fourier_space>no no no</fourier_space>
      <vectors> vc1 </vectors>
      <![CDATA[
        const double realfn = (mu-0.5*M*Vr(1))/Uoh11/hbar;
        
       for(long j=1; j<7; j++) {
         if (j==1||j==2||j==5)
          phi(j) = realfn>0. ? complex(sqrt(realfn),0) : complex(0,0);
        else
          phi(j) = complex(0,0);
 	}
      ]]>
    </vector>
  </field>

  <!-- The sequence of integrations to perform -->
  <sequence>
    <integrate>
      <algorithm>ARK89EX</algorithm>
      <interval>1e-7</interval>
      <tolerance>1.0e-7</tolerance>
      <lattice>1000</lattice>
      <samples>10 10 1</samples>
      <k_operators>
        <constant>yes</constant>
        <operator_names> L2p L2n L4p L4n </operator_names>
        <![CDATA[
          L2p = complex(0,-hbar/M/2/chi*(kx*kx+ky*ky+kz*kz));
          L2n = complex(0, hbar/M/2/chi*(kx*kx+ky*ky+kz*kz));
          L4p = complex(0,-hbar/M/4/chi*(kx*kx+ky*ky+kz*kz));
          L4n = complex(0, hbar/M/4/chi*(kx*kx+ky*ky+kz*kz));
        ]]>
      </k_operators>
      
        <moment_group>
        <moments>chippy</moments>
        <integrate_dimension>yes yes yes</integrate_dimension>
	    <![CDATA[
		chippy += ~phi(5)*phi(5);   
		]]>
        </moment_group>    
          
        <moment_group>
        <moments>ippy ichippy</moments>
        <integrate_dimension>no no no</integrate_dimension>
	    <![CDATA[
		ippy += phi(1);
		ichippy += phi(5);   
		]]>
        </moment_group>      
        <moment_group>
        <moments>py ic</moments>
        <integrate_dimension>yes yes no</integrate_dimension>
	    <![CDATA[
		py += phi(1);
		ic += phi(5);   
		]]>
        </moment_group>      
        <moment_group>
        <moments>ppy ichi</moments>
        <integrate_dimension>no yes yes</integrate_dimension>
	    <![CDATA[
		ppy += phi(1);
		ichi += phi(5);   
		]]>
        </moment_group>      
           
      <vectors> main vc1 </vectors>
      <![CDATA[
        const complex dens1 = phi(2)*phi(1);
        const complex dens3 = phi(4)*phi(3);
        
     const double gdens1 = (phi(5).re*phi(5).re+phi(5).im*phi(5).im);
     const double gdens3 = (phi(6).re*phi(6).re+phi(6).im*phi(6).im);
               
        dphi_dt(1) = L2p[phi](1) + (-i*Vr(2)-gamloss11
             +(gamloss132/2+gamloss12)/2/(dx*dy*dz))*phi(1) 
             + (-i*gameff-gamloss12)*dens1*phi(1) 
             - (i*gam13+gamloss132)*dens3*phi(1) -i*phi(2)*phi(3) 
             + i*chippy*ippy;
                             
        dphi_dt(2) = L2n[phi](2) + (i*Vr(2)-gamloss11
             +(gamloss132/2+gamloss12)/2/(dx*dy*dz))*phi(2) 
             + (i*gameff-gamloss12)*dens1*phi(2) 
             + (i*gam13-gamloss132)*dens3*phi(2) +i*phi(1)*phi(4);
                             
        dphi_dt(3) = L4p[phi](3) + (-i*Vr(3)-gamloss31
             +(gamloss132/2+gamloss32)/2/(dx*dy*dz))*phi(3) 
             + (-i*gam33-gamloss32)*dens3*phi(3) 
             - i*0.5*phi(1)*phi(1) -(i*gam13+gamloss132)*dens1*phi(3);
                             
        dphi_dt(4) = L4n[phi](4) + (i*Vr(3)-gamloss31+(gamloss132/2
           +gamloss32)/2/(dx*dy*dz))*phi(4) 
           + (i*gam33-gamloss32)*dens3*phi(4) 
           + i*0.5*phi(2)*phi(2) +(i*gam13-gamloss132)*dens1*phi(4);
                             
        dphi_dt(5) = L2p[phi](5) + (-i*Vr(4)-gamloss11)*phi(5) 
           +(-i*gameff-gamloss12)*gdens1*ichippy 
           - (i*gam13+gamloss132)*gdens3*phi(5)-i*conj(phi(5))*phi(6);
                              
        dphi_dt(6) = L4p[phi](6) + (-i*Vr(5)-gamloss31)*phi(6) 
           +(-i*gam33-gamloss32)*gdens3*phi(6) 
           - i*0.5*phi(5)*phi(5) +(i*gam13-gamloss132)*gdens1*phi(6);
      ]]>
    </integrate>
  </sequence>

  <!-- The output to generate -->
  <output format="binary" precision="double">
    <group>
      <sampling>
        <fourier_space>   no    no no</fourier_space>
        <lattice>         16    1  1</lattice>
        <moments>atoms molecules gatoms gmolecules</moments>
        <![CDATA[
          atoms=phi(2)*phi(1);
          molecules=phi(4)*phi(3);
          gatoms=conj(phi(5))*phi(5);
          gmolecules=conj(phi(6))*phi(6);
        ]]>
      </sampling>
    </group>
    <group>
      <sampling>
        <fourier_space>    no   no  no</fourier_space>
        <lattice>           0    16   0</lattice>
        <moments>rn_1 rn_2 grn_1 grn_2 excitedn</moments>
        <![CDATA[
    rn_1 = phi(2)*phi(1);
    rn_2 = phi(4)*phi(3);
    grn_1 = conj(phi(5))*phi(5);
    grn_2 = conj(phi(6))*phi(6);
    excitedn = g*g/4*phi(2)*phi(2)*phi(1)*phi(1)+F*F*phi(4)*phi(3)
     - F*g/2*(phi(2)*phi(2)*phi(3)+phi(1)*phi(1)*phi(4));
        ]]>
      </sampling>
    </group>
    <group>
      <sampling>
        <fourier_space>   no no no</fourier_space>
        <lattice>         4    8  16</lattice>
        <moments>atomsr moleculesr atomsi moleculesi</moments>
        <![CDATA[
          atomsr=phi(1);
          moleculesr=phi(3);
          atomsi=-i*phi(5);
          moleculesi=-i*phi(6);
        ]]>
      </sampling>
    </group>
  </output>

</simulation>
\end{xmdsCode}

This simulation is identical in function to the highdim.xmds example above, but describes the fields as an array of components rather than a list.  This notation may be very valuable when the numbers of component get very large and the equations can be easily described in terms of the index.

WARNING: There is no bounds checking on the index of your field, so be careful when writing your equations in this form.

When using this notation, if XMDS needs to calculate the k-space operator of any of the components of an array, all of them are calculated.  This makes, for example, highdim\_vector\_version.xmds slower than the old version.

IMPORTANT:  For interaction picture algorithms, if a k-space operator is applied to any component of a vector, then it is applied to ALL OF THEM.  This means that highdim\_vector\_version.xmds only solves the correct equations when used with an EX algorithm.