1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
|
\input{newcommands}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mesh refinement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[htb]
\centering
\includegraphics[width=15cm]{ICNSP_2011_Vay_fig1.png}
\caption{Sketches of the implementation of mesh refinement in WarpX with the electrostatic (left) and electromagnetic (right) solvers. In both cases, the charge/current from particles are deposited at the finest levels first, then interpolated recursively to coarser levels. In the electrostatic case, the potential is calculated first at the coarsest level $L_0$, the solution interpolated to the boundaries of the refined patch $r$ at the next level $L_{1}$ and the potential calculated at $L_1$. The procedure is repeated iteratively up to the highest level. In the electromagnetic case, the fields are computed independently on each grid and patch without interpolation at boundaries. Patches are terminated by absorbing layers (PML) to prevent the reflection of electromagnetic waves. Additional coarse patch $c$ and fine grid $a$ are needed so that the full solution is obtained by substitution on $a$ as $F_{n+1}(a)=F_{n+1}(r)+I[F_n( s )-F_{n+1}( c )]$ where $F$ is the field, and $I$ is a coarse-to-fine interpolation operator. In both cases, the field solution at a given level $L_n$ is unaffected by the solution at higher levels $L_{n+1}$ and up, allowing for mitigation of some spurious effects (see text) by providing a transition zone via extension of the patches by a few cells beyond the desired refined area (red \& orange rectangles) in which the field is interpolated onto particles from the coarser parent level only.}
\label{fig:ESAMR}
\end{figure}
The mesh refinement methods that have been implemented in WarpX were developed following the following principles: i) avoidance of spurious effects from mesh refinement, or minimization of such effects; ii) user controllability of the spurious effects' relative magnitude; iii) simplicity of implementation. The two main generic issues that were identified are: a) spurious self-force on macroparticles close to the mesh refinement interface \cite{Vaylpb2002,Colellajcp2010}; b) reflection (and possible amplification) of short wavelength electromagnetic waves at the mesh refinement interface \cite{Vayjcp01}. The two effects are due to the loss of translation invariance introduced by the asymmetry of the grid on each side of the mesh refinement interface.
In addition, for some implementations where the field that is computed at a given level is affected by the solution at finer levels, there are cases where the procedure violates the integral of Gauss' Law around the refined patch, leading to long range errors \cite{Vaylpb2002,Colellajcp2010}. As will be shown below, in the procedure that has been developed in WarpX, the field at a given refinement level is not affected by the solution at finer levels, and is thus not affected by this type of error.
\subsection{Electrostatic}
A cornerstone of the Particle-In-Cell method is that assuming a particle lying in a hypothetical infinite grid, then if the grid is regular and symmetrical, and if the order of field gathering matches the order of charge (or current) deposition, then there is no self-force of the particle acting on itself: a) anywhere if using the so-called ``momentum conserving'' gathering scheme; b) on average within one cell if using the ``energy conserving'' gathering scheme \cite{Birdsalllangdon}. A breaking of the regularity and/or symmetry in the grid, whether it is from the use of irregular meshes or mesh refinement, and whether one uses finite difference, finite volume or finite elements, results in a net spurious self-force (which does not average to zero over one cell) for a macroparticle close to the point of irregularity (mesh refinement interface for the current purpose) \cite{Vaylpb2002,Colellajcp2010}.
A sketch of the implementation of mesh refinement in WarpX is given in Figure~\ref{fig:ESAMR} (left). Given the solution of the electric potential at a refinement level $L_n$, it is interpolated onto the boundaries of the grid patch(es) at the next refined level $L_{n+1}$. The electric potential is then computed at level $L_{n+1}$ by solving the Poisson equation. This procedure necessitates the knowledge of the charge density at every level of refinement. For efficiency, the macroparticle charge is deposited on the highest level patch that contains them, and the charge density of each patch is added recursively to lower levels, down to the lowest.
\begin{figure}[htb]
\centering
\includegraphics[width=15cm]{ICNSP_2011_Vay_fig2.png}
\caption{Position history of one charged particle attracted by its image induced by a nearby metallic (dirichlet) boundary. The particle is initialized at rest. Without refinement patch (reference case), the particle is accelerated by its image, is reflected specularly at the wall, then decelerates until it reaches its initial position at rest. If the particle is initialized inside a refinement patch, the particle is initially accelerated toward the wall but is spuriously reflected before it reaches the boundary of the patch whether using the method implemented in WarpX or the MC method. Providing a surrounding transition region 2 or 4 cells wide in which the potential is interpolated from the parent coarse solution reduces significantly the effect of the spurious self-force. }
\label{fig:ESselfforce}
\end{figure}
The presence of the self-force is illustrated on a simple test case that was introduced in \cite{Vaylpb2002} and also used in \cite{Colellajcp2010}: a single macroparticle is initialized at rest within a single refinement patch four cells away from the patch refinement boundary. The patch at level $L_1$ has $32\times32$ cells and is centered relative to the lowest $64\times64$ grid at level $L_0$ (``main grid''), while the macroparticle is centered in one direction but not in the other. The boundaries of the main grid are perfectly conducting, so that the macroparticle is attracted to the closest wall by its image. Specular reflection is applied when the particle reaches the boundary so that the motion is cyclic. The test was performed with WarpX using either linear or quadratic interpolation when gathering the main grid solution onto the refined patch boundary. It was also performed using another method from P. McCorquodale et al (labeled ``MC'' in this paper) based on the algorithm given in \cite{Mccorquodalejcp2004}, which employs a more elaborate procedure involving two-ways interpolations between the main grid and the refined patch. A reference case was also run using a single $128\times128$ grid with no refined patch, in which it is observed that the particle propagates toward the closest boundary at an accelerated pace, is reflected specularly at the boundary, then slows down until it reaches its initial position at zero velocity. The particle position histories are shown for the various cases in Fig. \ref{fig:ESselfforce}. In all the cases using the refinement patch, the particle was spuriously reflected near the patch boundary and was effectively trapped in the patch. We notice that linear interpolation performs better than quadratic, and that the simple method implemented in WarpX performs better than the other proposed method for this test (see discussion below).
\begin{figure}[htb]
\centering
\includegraphics[width=15cm]{ICNSP_2011_Vay_fig3.png}
\caption{(left) Maps of the magnitude of the spurious self-force $\epsilon$ in arbitrary units within one quarter of the refined patch, defined as $\epsilon=\sqrt{(E_x-E_x^{ref})^2+(E_y-E_y^{ref})^2}$, where $E_x$ and $E_y$ are the electric field components within the patch experienced by one particle at a given location and $E_x^{ref}$ and $E_y^{ref}$ are the electric field from a reference solution. The map is given for the WarpX and the MC mesh refinement algorithms and for linear and quadratic interpolation at the patch refinement boundary. (right) Lineouts of the maximum (taken over neighboring cells) of the spurious self-force. Close to the interface boundary (x=0), the spurious self-force decreases at a rate close to one order of magnitude per cell (red line), then at about one order of magnitude per six cells (green line).}
\label{fig:ESselfforcemap}
\end{figure}
The magnitude of the spurious self-force as a function of the macroparticle position was mapped and is shown in Fig. \ref{fig:ESselfforcemap} for the WarpX and MC algorithms using linear or quadratic interpolations between grid levels. It is observed that the magnitude of the spurious self-force decreases rapidly with the distance between the particle and the refined patch boundary, at a rate approaching one order of magnitude per cell for the four cells closest to the boundary and about one order of magnitude per six cells beyond. The method implemented in WarpX offers a weaker spurious force on average and especially at the cells that are the closest to the coarse-fine interface where it is the largest and thus matters most.
We notice that the magnitude of the spurious self-force depends strongly on the distance to the edge of the patch and to the nodes of the underlying coarse grid, but weakly on the order of deposition and size of the patch.
A method was devised and implemented in WarpX for reducing the magnitude of spurious self-forces near the coarse-fine boundaries as follows. Noting that the coarse grid solution is unaffected by the presence of the patch and is thus free of self-force, extra ``transition'' cells are added around the ``effective'' refined area.
Within the effective area, the particles gather the potential in the fine grid. In the extra transition cells surrounding the refinement patch, the force is gathered directly from the coarse grid (an option, which has not yet been implemented, would be to interpolate between the coarse and fine grid field solutions within the transition zone so as to provide continuity of the force experienced by the particles at the interface). The number of cells allocated in the transition zones is controllable by the user in WarpX, giving the opportunity to check whether the spurious self-force is affecting the calculation by repeating it using different thicknesses of the transition zones. The control of the spurious force using the transition zone is illustrated in Fig.~\ref{fig:ESselfforce}, where the calculation with WarpX using linear interpolation at the patch interface was repeated using either two or four cells transition regions (measured in refined patch cell units). Using two extra cells allowed for the particle to be free of spurious trapping within the refined area and follow a trajectory that is close to the reference one, and using four extra cells improved further to the point where the resulting trajectory becomes undistinguishable from the reference one.
We note that an alternative method was devised for reducing the magnitude of self-force near the coarse-fine boundaries for the MC method, by using a special deposition procedure near the interface \cite{Colellajcp2010}.
%\begin{figure}[htb]
% \centering
% \includegraphics[width=15cm]{ICNSP_2011_Vay_fig4.png}
% \caption{Snapshot from a 3D self-consistent simulation of the injector in the High Current Experiment shows the beam emerging from the source at low energy (blue) and being accelerated (green-yellow-orange) and transported in a four quadrupole front end. The automatic layout of the mesh refinement patches from a 2D axisymmetric simulation of the source area shows 2 levels of refinement, concentrating the finer meshes around the emitter (white curve surface) and the beam edge (dark blue).}
% \label{fig:ESHCX}
%\end{figure}
%Automatic remeshing has been implemented in WarpX following the procedure described in \cite{Vaynim2005}, refining on criteria based on measures of local charge density magnitude and gradients. AMR WarpX simulations were applied to the modeling of the front end injector of the High Current Experiment (HCX) \cite{Prostprstab2005}, and provided the first numerically converged estimates of phase space beam distorsions, which directly affects beam quality \cite{Vaypop04}. Fig.~\ref{fig:ESHCX} shows snapshots from 2D axisymmetric simulation of the souce area illustrating the automatic placement of refined patches, and 3D simulation of the full injector showing the beam generation, acceleration and transport.
\subsection{Electromagnetic}
The method that is used for electrostatic mesh refinement is not directly applicable to electromagnetic calculations. As was shown in section 3.4 of \cite{Vayjcp01}, refinement schemes relying solely on interpolation between coarse and fine patches lead to the reflection with amplification of the short wavelength modes that fall below the cutoff of the Nyquist frequency of the coarse grid. Unless these modes are damped heavily or prevented from occurring at their source, they may affect particle motion and their effect can escalate if trapped within a patch, via multiple successive reflections with amplification.
To circumvent this issue, an additional coarse patch (with the same resolution as the parent grid) is added, as shown in Fig.~\ref{fig:ESAMR}-right and described in \cite{Vaycpc04}. Both the fine and the coarse grid patches are terminated by Perfectly Matched Layers, reducing wave reflection by orders of magnitude, controllable by the user \cite{Berengerjcp96,Vayjcp02}. The source current resulting from the motion of charged macroparticles within the refined region is accumulated on the fine patch and is then interpolated onto the coarse patch and added onto the parent grid. The process is repeated recursively from the finest level down to the coarsest. The Maxwell equations are then solved for one time interval on the entire set of grids, by default for one time step using the time step of the finest grid. The field on the coarse and fine patches only contain the contributions from the particles that have evolved within the refined area but not from the current sources outside the area. The total contribution of the field from sources within and outside the refined area is obtained by adding the field from the refined grid $F(r)$, and adding an interpolation $I$ of the difference between the relevant subset $s$ of the field in the parent grid $F(s)$ and the field of the coarse grid $F( c )$, on an auxiliary grid $a$, i.e. $F(a)=F(r)+I[F(s)-F( c )]$. The field on the parent grid subset $F(s)$ contains contributions from sources from both within and outside of the refined area. Thus, in effect, there is substitution of the coarse field resulting from sources within the patch area by its fine resolution counterpart. The operation is carried out recursively starting at the coarsest level up to the finest.
An option has been implemented in which various grid levels are pushed with different time steps, given as a fixed fraction of the individual grid Courant conditions (assuming same cell aspect ratio for all grids and refinement by integer factors). In this case, the fields from the coarse levels, which are advanced less often, are interpolated in time.
The substitution method has two potential drawbacks due to the inexact cancellation between the coarse and fine patches of : (i) the remnants of ghost fixed charges created by the particles entering and leaving the patches (this effect is due to the use of the electromagnetic solver and is different from the spurious self-force that was described for the electrostatic case); (ii) if using a Maxwell solver with a low-order stencil, the electromagnetic waves traveling on each patch at slightly different velocity due to numerical dispersion.
The first issue results in an effective spurious multipole field whose magnitude decreases very rapidly with the distance to the patch boundary, similarly to the spurious self-force in the electrostatic case. Hence, adding a few extra transition cells surrounding the patches mitigates this effect very effectively.
%[Add hyperbolic correction?]
The tunability of WarpX's electromagnetic finite-difference and pseudo-spectral solvers provides the means to optimize the numerical dispersion so as to minimize the second effect for a given application, which has been demonstrated on the laser-plasma interaction test case presented in \cite{Vaycpc04}.
Both effects and their mitigation are described in more detail in \cite{Vaycpc04}.
Caustics are supported anywhere on the grid with an accuracy that is set by the local resolution, and will be adequately resolved if the grid resolution supports the necessary modes from their sources to the points of wavefront crossing. The mesh refinement method that is implemented in WarpX has the potential to provide higher efficiency than the standard use of fixed gridding, by offering a path toward adaptive gridding following wavefronts.
%\begin{figure}[htb]
% \centering
% \includegraphics[width=13cm]{ICNSP_2011_Vay_fig5.png}
% \caption{Electron density $n_e$ (normalized to the density of the injected plasma) from WarpX simulations in 2-1/2D for a), b), c) and 3D for d) of a rigid beam (thin light-blue outline) propagating through a neutral plasma, for grid sizes of a) $128\times320$, b) $512\times1280$, c) $128\times320$ (main grid, red box) + $128\times640$ (patch 1, orange box) + $128\times1280$ (patch 2, yellow box), such that the resolution of patch 2 matched the resolution of the grid used for b), d) grid size of $64\times64\times160$ (main grid, red box) + $64\times64\times320$ (patch 1, orange box) + $64\times64\times640$ (patch 2, yellow box). For c) and d), the number and weight of injected plasma macroparticles was adjusted to keep the number of macroparticles per cell constant in each grid at injection in front of the beam.}
% \label{fig:EMplasma}
%\end{figure}
%As a test to the electromagnetic PIC implementation, WarpX simulations of wave excitations by a beam propagating through plasma, as described in \cite{Kaganovichpop2004}, were conducted. In these simulations, a hard-edged, elliptical, rigid beam propagates at constant velocity $v_z = 0.5c$ where $c$ is the speed of light through an initially cold neutral plasma of initial density $n_0$. The beam has a flat-top density profile of $n_b = n_0/2$, and an elliptical shape of length $l = 15 c/\omega_p$ and diameter $d = l/10$, where $\omega_p$ is the electron plasma frequency. It is shown in \cite{Kaganovichpop2004} that waves with a wavenumber of approximately $2\omega_p/v_z$ are generated in the plasma by the beam's electrostatic field, and have larger amplitude inside the beam, due to their interaction with the beam's sharp edges.
%Resolving the beam edge and the small structures developing in the wake inside the beam forces small cell sizes. The resolution that is needed for macroscopic convergence was explored in 2-1/2D in a series of four runs where the number of grid cells was varied from $64\times160$ to $512\times1280$ by incremental factors of 2. Third order spline interpolation was used for the beam and plasma macroparticle current deposition and force gathering. The details of the plasma wake were very similar between the two highest resolution cases, indicating that macroscopic convergence was reached. The results from the runs using $128\times320$ and $512\times1280$ grids are shown in Fig.~\ref{fig:EMplasma}. The result from the highest resolution run serves as the reference for subsequent calculations with mesh refinement.
%A run was conducted where the main grid had $128\times320$ cells and was complemented by two refinement patches (with successive refinement factors of 2 in each direction), such that the resolution in the central patch matched the resolution of the case of reference.
%The number and weight of the injected plasma macroparticles was varied, such that the number of macroparticles per cell in each grid at injection was constant. Results are plotted in Fig.~\ref{fig:EMplasma} (bottom-left) showing a good reproduction of the fine scale structures within the central fine patch in good agreement with the reference case.
%Lastly, a three-dimensional simulation with mesh refinement of the same physical setup was conducted. The grid setup and 3D isosurfaces of the plasma electron density as the beam enters the plasma are shown in Fig.~\ref{fig:EMplasma} (bottom-right). As expected, structures similar to the ones observed in 2D are present within the beam envelope. The speedup achieved by the use of mesh refinement was estimated to be approximately one order of magnitude in 3D.
|
