Preface: This is in response to the request of the ILC GDE for a
discussion of the strengths of the VORPAL parallel EM framework
for the computations needed for the ILC.  As background, we note
that there are two codes proposed for such work, the VORPAL
framework, which can incorporate algorithms for arbitrary
hexahedral meshes, and SLAC finite element suite (Omega3p, Tau3p,
etc.) which work on unstructured meshes.

     Applicability of VORPAL to Electromagnetics Calculation

Summary: Recent algorithm additions have made the VORPAL
framework [1] well suited for complex electromagnetic
calculations such as those needed for computational Grand
Challenges like ILC cavity modeling.  Those additions include
embedded boundary methods and python-scripted geometry
generation.  These new capabilities augment the previously
demonstrated features of VORPAL: fully self-consistent particle
modeling, scalability to 1000s of processors, and modest memory
requirements compared to global algorithms.  Object-oriented
design allows developers to quickly add new capability to VORPAL,
for example perfect dispersion algorithms and implicit
time-stepping algorithms. VORPAL is presently installed at and
available to researchers at FNAL, LBNL, JLAB, BNL, and ANL.

The VORPAL framework computes the evolution for partial
differential equations on a hexahedral mesh.  It is a general
framework that allows for particles, fields, and their
interactions.  Of most interest here is the finite-difference
(FD), time-domain (TD) electromagnetics including particle
effects through the Particle-In-Cell methods thus allowing the
computation of the propagation of beams fully self-consistently.
That is, particles respond to both the electromagnetic fields
powered by the klystrons as well as those produced by the beam
particles. VORPAL has had great success in the laser-wake-field
research, and it now has demonstrated capability to handle
complex boundary conditions needed for the ILC.  Recent
accomplishments have been documented in multiple places,
including the proceedings of the 2006 SciDAC meeting, as well as
recent invited talks at FNAL and the just held International
Computational Accelerator Physics Workshop.  VORPAL
visualizations have appeared on the cover of Nature and of the
proceedings of SciDAC 2006.

The advantages of the algorithms implemented within the VORPAL
framework are their speed, scalability, and minimal memory usage.
The speed comes from the fact that all updates can be written as
local operations on grid variables, involving only the nearest
neighbors.  (There are no global linear equations to solve.) The
scalability then derives from this, as such updates require the
transfer of only the data on the surfaces of the domains assigned
to each processor.  To take full advantage of this, VORPAL has
within it a set-theory based approach to setting up
communications that allows an arbitrary covering decomposition.
The result has been demonstrated scalability up the the largest
supercomputers available at the time of the last tests, namely to
4000 processors on Seaborg, the NERSC POWER III supercomputer.
Finally, as to memory usage, the algorithms and structures within
VORPAL ensure that almost no more data than exists with the
fields and particles is needed. Thus, it is possible to run a
problem with 1.5B degrees of freedom with no more than 24 GB of
memory [a].

The FD methods have made important strides in the last decade
through the development of embedded boundaries (aka cut cells or
conformal boundaries). Important early work was that by Dey and
Mittra [2] and this has been refined in the work of Zagorodnov et
al [3].  In these papers it was shown that one could obtain high
accuracy in the presence of curved boundaries. VORPAL has (to the
best of my knowledge) the only parallel implementation of
conformal boundaries in a finite difference computational EM
application.  As noted in our presentation at SciDAC 2006, this
has allowed the use of Richardson extrapolation to increase all
computations by one higher order.  Thus, the implemented
second-order algorithm produces results with third-order error.
This allows one to obtain accuracy to on part on 10^5 with only
about 50 cells per characteristic dimension.

The above results rely on the ability to generate meshes, and
here is another place where embedded boundary methods have a real
advantage.  With embedded boundaries, one need not compute the
mesh throughout a volume, only at the surface, where one
essentially computes only the intersections of the boundary with
the mesh and from that computes fractional lengths, areas, and
volumes of cells.  Consequently, mesh generation parallelizes
trivially, and mesh generation is very fast - typically minutes
at most on the largest problems - sufficiently fast that for
convenience, we regenerate the mesh at each run because the
computation is so minor.

We currently generate meshes from a python-coded logical function
that returns whether a point is inside the cavity.  This implies
that cavities can be defined through a scripting process.  (We
have an STL-CAD implementation mostly implemented but not fully
tested.)

The above two capabilities imply the ability to do error and
parameter studies rapidly.  For example, to add a deformity to a
cavity, one modifies the boundary function through addition of an
error term and reruns the simulation.  There is no need to go
through the cycle of CAD-description-to-mesh-generation-to-
code-importation to do a calculation.  Streamlining this cycle
becomes important for the problem in which one would like to
infer the imperfections from the measured spectrum.

There has been considerable work in algorithms for FD
electromagnetics, including the absolutely stable methods for
time integration.  These include the work of deRaedt [4], Zheng
et al [5], with the latter recently followed up by Fornberg [6].
The first of these is fully explicit (no linear solve), while the
latter are methods that require at each stage of the algorithm
only to solve a set of tridiagonal systems.  Such solves converge
in a single iteration with only 6 times the number of numerical
operations as an explicit update.  Our own work [7] has provided
a modification of the Zheng et al and Fornberg methods that is
charge conserving, so that Poisson's equation is satisfied. Thus,
we can now have absolute stability with implicit methods that
require minimal work.  This will be important in SciDAC-2 as we
model systems in which the needed transverse resolution, dictated
by the beam size, is much greater than the longitudinal.

Yet another recent advance has been in the development of perfect
dispersion algorithms [8], in which the dispersion of
perturbations aligned with the beam perfectly follows the vacuum
result of omega = kc.  This is important for the modeling of
relativistic beams in cavities, as the emission of radiation,
hence the wake fields, depends sensitively on phase velocity of
light relative to the velocity of the beam.  Most numerical
methods produce short wavelength subluminal waves, which are then
artificially radiated by a highly relativistic beam.  This
"numerical Cerenkov" radiation results in unphysical radiation
generated by the beam and ultimately traveling down the pipe or
out the HOM couplers.  Superluminal waves can be problematic as
well, as the interaction between beam and the mode changes,
though not as dramatically.  The perfect dispersion algorithms
eliminate these unphysical effects completely.  Finally, we have
now shown [9] how to modify our more rapid new algorithm [7] to
obtain perfect dispersion and charge conservation simultaneously.
We believe that this algorithm will be crucial for end-to-end,
self-consistent modeling of the ILC.

Higher-order algorithms are available as well.  A recent,
promising approach to combining those with conformal boundaries
is presented in the compact stencil work of [10].

While VORPAL has so far been used in only the time domain,
implementing a frequency-domain model is straight forward. VORPAL
already has implementations for the Bowers [11] implicit advance
through use of matrix methods that call out to the Trilinos
parallel linear algebra library. While the code is not in place
for conformal boundaries (it is present for planar boundaries),
this could be done rapidly in SciDAC-2, were it declared a
priority, and this would provide another check on the
frequency-domain computations currently being carried out by the
SLAC electromagnetics group.

Additionally, the VORPAL framework was built very generally, to
be able to incorporate any algorithms that can exist on a
hexahedral mesh.  Thus, finite-element algorithms can be
implemented within  the VORPAL framework as well, should the
case be made.

Finally, VORPAL is presently installed at many of the national
laboratories, including FNAL, LBNL, JLAB, BNL, and ANL, and at
several universities, including the Universities of Texas and
California.  Scientists at these labs have full access to VORPAL
to assist them with their research.

[1] C. Nieter and J. R. Cary J. Comp. Phys. 196, 538 (2004).

[2] S. Dey and R. Mittra, IEEE Microwave and Guided Wave Lett. v.
7, 273 (1997).

[3] I. A. Zagorodnovn, R. Schuhmann and T. Weiland, Int. J.
Numer.  Model: Elec net, Devices, Fields v. 16, 127 (2003).

[4] J. S. Kole, M. T. Figge, and H. De Raedt, Phys. Rev. E v. 64,
066705 (2001).

[5] F. Zheng, Z. Chen, and J. Zhang, IEEE Trans. Microwave Theory
Techniques v. 48, 1550 (2000).

[6] J. Lee and B. Fornberg, J. Comp. App. Math. v. 66, 497
(2003).

[7] D. Smithe, J. R. Cary, J. A. Carlsson, in draft.

[8] I. Zagorodnov and T. Weiland, Phys. Rev. ST/AB v. 8, 042001
(2005).

[9] D. Smithe, J. R. Cary, J. A. Carlsson, in draft.

[10] A. Yefet, E. Turkel, App. Num. Math. v. 33, 125 (2000).

[11] K. Bowers, Ph.D. thesis, U.  Calif. (2001).