Domain Decomposition Framework for Maxwell Finite Element Solvers and Application to PIC

The charge and current density is defined using the PSDF conventional definition, $\rho(t,\mathbf{r})=q\int_{\Omega}f(t,\mathbf{r},\mathbf{v})d\mathbf{v}$ and $\mathbf{J}_{\rho}(t,\mathbf{r})=q\int_{\Omega}\mathbf{v}f(t,\mathbf{r},\mathbf{v})d\mathbf{v}$. The PSDF is approximated by using $N_{p}$ samples, or particles, with shape $S(\mathbf{r})$ such that the charge and current density are defined as $\rho(t,\mathbf{r})=q\sum_{p=1}^{N_{p}}S(\mathbf{r}-\mathbf{r}_{p}(t))$ and $\mathbf{J}_{\rho}(t,\mathbf{r})=q\sum_{p=1}^{N_{p}}\mathbf{v}(t)S(\mathbf{r}-\mathbf{r}_{p}(t))$ where $\mathbf{r}_{p}(t)$ and $\mathbf{v}_{p}(t)$ are the position and velocity of particle $p$. In this work, the shape functions are delta functions, though other shape functions are also valid and can be used to improve noise quality. The proper function to measure the charge density is the 0-form and the 1-form for the current density. Different schemes to evolve Newton’s equations can be used to preserve certain quantities like phase. In this work, to maintain accuracy of the time evolution at larger time step sizes, the Adams-Bashforth scheme is used. Lastly, the time integral of the current density

is used to allow the particle current to evolve consistently with the fields. Otherwise, there is a disconnect between the continuity equation and Gauss’ law.

The first formulation presented uses a Neumann-like boundary condition to enforce the continuity conditions of the fields at the interior boundaries,

This is accomplished by ensuring $\mathbf{x}^{n}_{i}=\mathbf{x}^{n}_{j}|_{\Gamma_{ij}}$ and defining

The system of equations to be solved in each sub-region

where

It is noted that the time history of the Lagrange multiplier is not explicitly stored, unlike the other currents in the system. Let the unknowns defined on $\Gamma_{c}$ be denoted by the subscript $c$ and all other unknowns in $\Omega_{i}$ denoted by subscript $r$, such that (7) can be written as

where $\bar{\bar{A}}=.5\bar{\bar{A}}_{1}+.25\Delta_{t}\bar{\bar{A}}_{0}$. The first set of equations can be rewritten as

which can be substituted into the second set of equations

Define a signed Boolean matrix $\bar{\bar{B}}_{r}^{(i)}$ that maps unknowns in $\Omega_{i}$, excluding those on $\Gamma_{c}$, to a set of unique global unknowns defined on $\cup\Gamma_{ij}$. The definition of $\bar{\bar{B}}_{r}^{(i)}$ is such that $\sum_{i=1}^{N_{v}}\bar{\bar{B}}_{r}^{(i)}\bar{x}^{(i)}=0$, which enforces (5a) and (5b). Furthermore, the transpose $\bar{\bar{B}}_{r}^{(i),T}$ forms a map from the global set of boundary unknowns to those of a particular volume. A similar unsigned Boolean matrix can be defined for unknowns that lie on $\Gamma_{c}$ within $\Omega_{i}$, $\bar{\bar{B}}_{c}$, which maps to a set of unique global unknowns. Using $\bar{\bar{B}}_{r}$ on (10) and summing over all volumes results in

where

and $\bar{X}_{c}$ the degree of freedom vector for corner unknowns.

The boundary condition (5a) is enforced at the corners of $\Omega_{i}$ as $\bar{\bar{B}}_{c}^{(i)}$ is defined so that $\sum_{i=1}^{N_{v}}\bar{\bar{B}}_{c}^{(i),T}\lambda^{(i)}_{c}=0$. Using that definition, (11) becomes

where

As in the vector wave equation domain decomposition formulations, $\bar{\bar{K}}_{cc}$ is sparse and generally much smaller than $\bar{\bar{K}}_{rr}$. Therefore, if solving the boundary problem directly, it is possible to define

and then use substitution into (12) to yield

where $\bar{\bar{\mathcal{K}}}_{rr}=\left(\bar{\bar{K}}_{rr}+\bar{\bar{K}}_{rc}\bar{\bar{K}}_{cc}^{-1}\bar{\bar{K}}_{cr}\right)$. If using an iterative solver, (12) and (14) can be solved as a coupled system instead.

In the previous formulation, a unique Lagrange multiplier was defined for each edge and face on $\Gamma$. This formulation enforces field continuity though a Neumann-like boundary condition. A more robust enforcement of field continuity is to use a Robin boundary condition [jin2015finite]. The boundaries between $\Omega_{i}$ and $\Omega_{j}$, $\Gamma_{ij}$ and $\Gamma_{ji}$, are distinct and have a unique Lagrange multiplier defined for each subdomain. Instead of (6), the boundary condition on the Lagrange multipliers becomes

The edges on $\Gamma_{ij}$, as well as the interior unknowns and unknowns on $\partial\Omega_{i}$ are included in a vector $\bar{X}^{(i)}_{I}$. The Lagrange multipliers for the corner edges remain the same, with a unique Lagrange multiplier for each feature. The faces on an $\Gamma_{ij}$ are grouped together with the corner edges, which we will denote as $\bar{X}_{N}$. The system of equations solved in each subdomain is

where given

$\mathcal{R}^{(i)}_{I}$ contains $\mathcal{R}$ associated with interior unknowns, external unknowns, and unknowns of edges on $\Gamma_{ij}$ and $\mathcal{R}^{(i)}_{N}$ contain the remaining unknowns in subdomain $\Omega_{i}$. To derive the set of equations to solve for the dual unknowns, first consider the summation of (18) for two volumes

Define an unsigned Boolean matrix such that $\bar{\bar{B}}_{I}^{(i)}\bar{X}_{I}=\bar{e}_{I}^{(i)}$, the unknowns for the electric field on $\Gamma_{ij}$. Also, define the map from the boundary edges on $\Gamma_{ij}$ to its counterpart on $\Gamma_{ji}$, $\bar{T}_{ij}$. Testing (21) with curl-conforming basis functions $\mathbf{W}^{e}$, generalizing the case to any number of subdomains, and using $\bar{\bar{B}}_{I}^{(i)}$ and $\bar{\bar{T}}_{ji}$ restrict and map the unknowns of $\Omega_{j}$ to $\Omega_{i}$ yields

where $\sigma(i)$ are all subdomains $\Omega_{j}$ which share a boundary with $\Omega_{i}$ and the matrix $\bar{\bar{M}}_{T}^{(i)}\langle\hat{n}\times\mathbf{W}^{e},\hat{n}\times\mathbf{W}^{e}\rangle|_{\Gamma_{ij}}$.

Rewrite the first set of equations in (19) as

Now, (23) can be plugged into (22) to yield

In order to complete the definition of the equation to solve for $\Lambda_{I}$ over $\cup\Gamma_{ij}$, we first define the equation for $\bar{X}_{N}$, the unknowns associated with faces on $\Gamma_{ij}$ and edges on $\Gamma_{c}$.

The derivation of the equation to solve for $\bar{X}_{N}$ follows a similar path as $\bar{X}_{c}$ from MFET-DP1. First, substitute (23) into the second set of equations in (22) to yield

The Boolean matrix $\bar{\bar{B}}_{N}^{(i)}$ is defined as

to enforce (5b) for the faces on $\Omega_{i}$ and (5a) edges on corners where $\bar{\bar{B}}_{r,f}^{(i)}$ is the portion of $\bar{\bar{B}}_{r}^{(i)}$ that acts on the face degrees of freedom. Using a similar process as (14), the equation for

where

The global system of equations for the edge Lagrange multipliers Define a Boolean matrix $\bar{\bar{B}}_{L}^{(i)}$ that maps the local Lagrange multipliers to a global vector. The Boolean matrix is applied to (24) and discretized in time to define the equation

where

It can be solved as a coupled system with (27) or in a Schur complement system similar to (17) using

where $\bar{\bar{\mathcal{K}}}_{II}=\left(\bar{\bar{K}}_{II}-\bar{\bar{K}}_{IN}\bar{\bar{K}}_{NN}^{-1}\bar{\bar{K}}_{NI}\right)$. However, $\bar{\bar{K}}_{NN}$ is much larger than its MFETI-DP1 counterpart $\bar{\bar{K}}_{cc}$ and will take more time and memory to compute and store its inverse.

The proposed domain decomposition formulations reduces the cost of solving for the fields by defining a set of equations on fictitious interior boundaries of the true problem. The definitions of these equations uses small matrix inverses, and therefore are generally not sparse. Furthermore, those matrix inverses have to be stored in memory. At the very least, one matrix inverse for each subdomain is needed, making the memory storage more than solving for the original system. However, it is expected that for direct inverses, inverting either the coupled equations of (12) and (14) or (29) and (27) would be more efficient than solving the original problem. If storing the inverse of the matrix associated with $\Gamma_{c}$ is possible, (17) and (31) are even smaller in size. The following discussion uses a sample problem to further probe properties of MFETI-DP1 and MFETI-DP2. In particular, we examine the performance of the method with a sample geometry with several configurations of subdivision.

Define $\Omega_{rect}$ as a rectangular region in free space with dimension .25 m $\times$ .2 m $\times$ .3 m discretized with average edge length of 5.5 cm. The region is the split into subdomains in several configurations listed in Table I so that MFETI-DP1 and MFETI-DP2 can be applied, where $\alpha$ is the ratio of interior unknowns to the number of boundary unknowns. This table is provided to provide the connection between the number of subdomains and $\alpha$ for the test example, as $\alpha$ will be more extensible to other geometries than a set number of subdomains. Figure 2 shows the condition number of the matrix formed by the coupled equations (12) and (14) for MFETI-DP1, (29) and (27) for MFETI-DP2, as well as the Schur complement forms $\bar{\bar{\mathcal{K}}}_{rr}$ and $\bar{\bar{\mathcal{K}}}_{II}$ with different numbers of subdivisions compared to the condition number of of the monolithic system matrix $\bar{\bar{A}}$. Solving the coupled equations results in condition number that are generally at or worse than the monolithic solve. Using the Schur complement results in condition numbers comparable to the monolithic solve for MFETI-DP1 and several orders of magnitude smaller for MFETI-DP2.

Typically, MFETI-DP1 and MFETI-DP2 when used with iterative solvers converge faster than the monolithic solve. Consider a domain $\Omega_{rect}$ is illuminated by an incident field defined by

The bandwidth of the modulated Gaussian is $\sigma=3(\pi f_{bw}$, where the bandwidth $f_{bw}=f_{max}-f_{0}$ with maximum frequency $f_{max}=500$MHz and center frequency $f_{0}=300$MHz. The time step size is $0.333$ ns which is ten times larger than what could be taken by a similar FDTD solver or leapfrog method for MFEM. Figure 3 shows the iteration counts for the coupled solves for MFETI-DP1 and MFETI-DP2 compared to the monolithic solve using GMRES with a tolerance of $10^{-6}$ and restarted after 20 iterations when $\Omega_{rect}$ is illuminated by an incident field. With a diagonal preconditioner, both domain decomposition approaches perform better than the monolithic solve, with MFETI-DP2 performing best overall.

Next, as alluded to earlier, the principal advantage of Newmark-$\beta$ schemes is the unconditional stability. As a result, it is important to verify that the both MFETI-DP1 and MFETI-DP2 retain this feature. To ascertain this, we examine the eigenspectra of these formulations to understand the stability of the time marching scheme. For the time marching scheme used in this work to be stable, the matrix used to evolve the system in time must have eigenvalues $\real(\lambda)>=0$ for the system to be unconditionally stable. Figure 4 shows that for solving both the coupled system and using the Schur complement approach, the domain decomposition methods will be stable. From the insert in Figure 4, the effect of using the schur complement shifts the eigenvalues of $\bar{\bar{\mathcal{K}}}_{rr}$ and $\bar{\bar{\mathcal{K}}}_{II}$ away from the imaginary axis, though MFETI-DP2 is shifted further away than MFETI-DP1.

The first example shows performance of the methods and difference in the fields. In this example, a $0.25m\times 0.2m\times 1m$ box is illuminated an incident field defined by

The bandwidth of the modulated Gaussian is $\sigma=3(\pi f_{bw}$, where the bandwidth $f_{bw}=f_{max}-f_{0}$ with maximum frequency $f_{max}=595$MHz and center frequency $f_{0}=295$MHz. The time step size $\delta_{t}=(150f_{0})^{-1}$ and the experiment is run for 1500 time steps. The boundary conditions on the external boundaries are robin boundary conditions. Three approaches are compared: monolithic (no DDM), MFETI-DP1, and MFETI-DP2 using three levels of discretization with average edge lengths of $h_{1}=5.20$cm, $h_{2}=3.63$cm, and $h_{3}=2.59$cm. The DDM meshes were split using METIS, such that no volumes have a regular shape. Figure 5, shows that the fields obtained through domain decomposition for edge lengths $h_{1}$, $h_{2}$, and $h_{3}$ respectively. Regardless of the number of subdivisions, the methods remain accurate to near machine precision with respect to the monolithic solve. This implies that using either domain decomposition formulation would cause the same acceleration on particles as the fields that are similar to each other within machine precision.

The next result demonstrates the effect of MFETI-DP on simulation time. As seen in Figure 6, both formulations are faster than the original monolithic solve for any number of subdivisions. However, there is an optimal ratio of interior volume unknowns to boundary unknowns for the fastest runtime. For each case, there is a crossover point at which further subdividing $\Omega$ incurs a cost in solving the boundary problem that is greater than the savings by reducing the unknowns in each volume.

The next two examples demonstrate MFETI-DP1 and MFETI-DP2 with particles. In these experiments, the fields that are obtained through domain decomposition are used to push the particles, which are not subject to any decomposition scheme. The first example is a particle beam traveling through a cylindrical PEC cavity. This case is a standard test to show charge conservation as errors readily appear as striations in the particle distribution. Furthermore, approximate solutions can be used to further verify the result. The first test is done on a cavity that is 10 cm in length and 0.02 cm in radius disctretized with an average edge length of 5 mm resulting in roughly 21,000 unknowns. The ten particles at each time step enter at the bottom face of the cavity with a radius of 8 mm. The beam voltage and current are 500 V and 50 mA respectively, causing the an initial velocity $1.326\times 10^{4}$ km/s. The geometry was subdivided in four configurations, the details of which are listed in Table II. The experiment was run for 1000 time steps at a time step size of 1 ns, which for the monolithic case took 1,342 seconds. As shown in figure 7, all the simulations maintained charge conservation over the course of the run. There is no notable difference in the simulations, regardless of the number of subdivisions used. Seven particles were tracked over the course of the run to track the particle trajectories with the domain decomposition. Because the fields are virtually equivalent, the particle trajectories are consistent across the different geometries. Figure 8 shows the error in the tracked particles’ path with rep sect to the path in the monolithic solve.

The next result was a particle beam in a larger PEC cavity. In this case, the PEC cylinder was 40 cm in length and 10 cm in radius. The geometry has an average edge length of 7.38 mm and a total of 218,059 unknowns. The beam voltage and current are 3.2 kV and 50 mA respectively, resulting in an initial velocity of $3.35\times 10^{4}$ km/s. The experiment was run for both MFETI-DP1 and MFETI-DP2 using the iterative solver GMRES using a tolerance of $10^{-6}$. The particle distribution is compared to a well-validated axialsymmetric EM-FDTDPIC code, XOOPIC [32] and shows good agreement in Figure 9.

The last example presented is the case of an adiabatic expansion of a plasma ball using MFETI-DP2. This example takes a considerable number or degrees of freedom to resolve the debye length and allow the particles to expand for an extended period of time. However, this can be made more tractable by using either domain decomposition schemes and unconditionally stable time marching to evolve the system with larger time steps than a similar EM-FDTDPIC formulation would allow. Additionally, this experiment has analytic solutions [33] that can be used for validation. A charge neutral, Gaussian distribution of 24000 particles is set at the center of a spherical region 24 cm in diameter. The $Sr^{+}$ ions are at 1K and the electrons are at 100K with a density of $5\times 10^{8}$ particles per cubic meter. The outer boundary has an impedance boundary condition and the boundary is sufficiently far from the particles that none leave during the simulation. The average edge length is 1 cm with a total of 681,214 unknowns. Figure 10 shows good agreement of the particle distribution with the analytic result at 150 ns, 300 ns, and 800 ns.