implicitCylEmUpdater

implicitCylEmUpdater

Works with a VSimPD license.

implicitCylEmUpdater is an implicit method that solves for the change in electric field \(\Delta E_{\theta}\). The method works by solving Ampere’s equation for \(\Delta E_{\theta}\). We then update the total electric field via \(E^{new}_{\theta} = E^{old}_{\theta} + \Delta E_{\theta}\). Then Faraday’s law is used to update the magnetic field. Of course, these steps must be centered in time to preserve 2nd order accuracy. This updater only works in a cylindrical 2D geometry and therefore solves for \(E_{\theta}\), \(B_r\), and \(B_z\).

This updater allows for the inclusion of dielectrics. The only allowed simulation boundary conditions are Dirichlet in \(\Delta E_{\theta}\) which is the same as setting the electric field on the boundaries. Under most conditions, conducting boundary conditions implies \(E=0\), though the user is free to choose other values.

The following equation is solved:

\[\begin{split}\Delta E_{\theta} - \frac{(\Delta t)^2}{4\epsilon \mu_0}\nabla^2 \Delta E_{\theta} = rhs \\ rhs = \frac{\Delta t}{\epsilon \mu_0} \left[-\mu_0 J_{\theta} + \left(\nabla \times \vec{B} \right)_{\theta} + \frac{\Delta t}{2}\nabla^2 E_{\theta} \right]\end{split}\]

This is a Poisson-like equation for \(\Delta E_{\theta}\) and therefore is solved via an iterative matrix inversion method. Once \(\Delta E_{\theta}\) is computed, the total electric and magnetic fields are updated.

implicitCylEmUpdater Basic Parameters

The implicitCylEmUpdater kind takes the lowerBounds and upperBounds parameters of FieldUpdater (typically set to [0, 0, 0] and [NX+1, NY+1, NZ+1], the bounds of the entire simulation including the uppermost nodes; however if the simulation is periodic in, e.g., the y direction, then the upper bound should typically be NY, not NY+1), as well as the following parameters:

writeFields (required string vector: [dE]): A list of Field objects containing the theta component of the change in electric field (at grid edges, units of V/m) to update.

readFields (required string vector: [rhs])

A list of Field objects containing the right hand side (rhs) of the Poisson-like equation. (at grid edges, with units of V/m).

The boundary potentials will be specified by <Conductor> blocks. In this case, the <Conductor>s will be equipotentials, and the electric field within the Conductor Block.

<LinearSolver> (code block, required): A code block describing the linear solver used to solve the Poisson equation: cf. LinearSolver. See the example below: solvers for symmetric matrices are appropriate for this problem, and for large (parallel) problems an iterativeSolver with a multigrid preconditioner may offer the best performance.

implicitCylEmUpdater Parameters (setting equipotentials on Conductors)

<Conductor> code blocks: Each Conductor Block describes the the shape of a conductor, and how the electric field is determined on that conductor. An arbitrary number of <Conductor> blocks can be used.

implicitCylEmUpdater Parameters (setting a Dirichlet boundary condition)

<Conductor> code blocks.: The Dirichlet boundary condition could be imposed on a boundary using the Conductor block. Inside the conductor block, a <ConductorShape> block could be used to select the boundary.

implicitCylEmUpdater Parameters (defining a Dielectric shape within the simulation domain.)

<DielectricShape > code blocks.: The ‘DielectricShape’ block could be used to define a dielectric region within the simulation domain.

enforceMatrixSymmetry = 0: The parameter specifies whether or not to enforce matrix symmetry. ‘enforceMatrixSymmetry’ takes 1 or 0 to enforce or not. Note that the matrix is enforced to be symmetric by default.

implicitCylEmUpdater Usage Notes

When particles are present near conducting surfaces and accuracy is important, the nodal electric field should be extrapolated (to nodes just inside the metal regions) using the deyMittraFieldExtrapolator, so that particles experience correct fields near metal surfaces.
It is best not to let <Conductor>s overlap. If they do overlap, the updater will issue a WARNING and do its best to yield a reasonable solution. However, the updater may calculate edge lengths inaccurately. Of course, this situation is actually unphysical if the <Conductor>s have different voltages; if they have the same voltage, it’s usually best to replace them with a single <Conductor>.
It is okay to have two conductors get very close to each other (e.g., much closer than a grid cell size). In general (of course) the simulation cannot be expected to resolve features smaller than a grid cell; however, cases like two large spheres (or parallel plates) that almost touch should be handled correctly if the spheres (or plates) are different <Conductor>s.
For performing repeated poisson solves with non-zero space charge, setting convergenceMetric=rhs in the iterativeSolver may reduce the solve time if the previous solution is similiar to the next solution. (Using rhs with zero space charge is not recommended, because the right-hand-side will be zero except near conducting boundaries with non-zero voltage, making it a poor choice for normalizing the residual error.)

implicitCylEmUpdater Advanced Parameters

The following parameters have default values that are appropriate for the majority of simulations.

exactRestart (default false)

If true, then restarting a simulation should produce exactly the same results as if the simulation hadn’t been halted and restarted. If false, then (at each timestep) the previous solution will be used as a guess for the next solution – this reduces the time to solution (if convergenceMetric=rhs in the iterativeSolver), but prevents the restart from being exact (although the previous solution is nearly reconstructed, it cannot always be exactly reconstructed). Restarting with exactRestart=false still yields an equally valid solution – but a solution that differs (by machine-precision levels, or by the solver tolerance) from the non-restarted solution.

Important: to work, this option must be set to true in the original run and also in the restarted simulation.

shiftNodesNearMetalToMetal (optional boolean, default true)

This option specifies how to deal with the difficult case of tiny edges, when the boundary surface cuts through a cell edge extremely close to the cell node (i.e., almost all of the edge is in metal). An edge is tiny (in this sense) if its fraction in dielectric is less than minEdgeFrac, which is by default the square root of the machine precision. Tiny fractional edges are a problem because the potential difference across them is inevitably tiny, so that subtracting the potentials (e.g., to find the electric field) can result in extreme precision loss, yielding an essentially random value for the potential difference (hence for the electric field).

If true, then tiny edges are replaced with all-metal edges, effectively shifting the metal boundary slightly so that it goes through the node. This is the recommended option; it’s usually the most robust and convenient option. However, it does change the metallicity of some nodes and edges: i.e., some nodes that are technically (but just barely) in dielectric will be considered to be in metal.

If false, then the metal boundary will be shifted slightly away from the node until all edge fractions are greater than minEdgeFrac. This option should be used when it’s important not to change the metallicity of any node.

The effect of shifting the metal boundary introduces some error to the solution, but that error (especially in double precision) will typically be orders of magnitude lower than the finite difference error due to finite cell size.

minEdgeFrac (optional float, default: square root of machine
precision): the minimum edge fraction allowed in dielectric. Any edges with smaller non-zero fractions in dielectric will either be reduced to zero, or increased to at least this value, depending on the samp:shiftNodesNearMetalToMetal attribute.

messagePotentialField (optional float, default true): whether to message the field (the first writeField) across MPI ranks, setting guard cell values. If true, then the UpdateStep containing this updater does not need to mention the potential field, though it won’t hurt to do it twice. If this is false, then the accuracy test (cf. accuracyTestStepSlice) may give incorrectly high values on domain boundaries.

asymmetryWarningLevel (optional float, default 0.005): When using floating potentials, capacitance and conductance matrices must be formed; they should be symmetric, but sometimes can be asymmetric, either due to finite precision truncation or to input file mistakes (e.g., letting conductors overlap). In these cases, a warning will be issued if the asymmetry (relative to the maximum element of the matrix) surpasses this level; in any case, the matrix will be symmetrized. Small asymmetries arise naturally from finite precision arithmetic; however, large asymmetries may indicate a problem interpreting the given geometry correctly (e.g., due to overlapping conductors).

accuracyTestStepSlice (optional vector of integers, default

[20, 2147483647, 1000])

a python-like ‘slice’ specifying the timesteps [start, stop, step] at which the accuracy of the Poisson solve will be tested and printed out (search output for “implicitCylEmUpdaterSolver [updater] accuracy”). This is a relatively time-consuming operation, so it should be done infrequently. The default will test accuracy at steps 20, 1020, 2020, 3020, etc.

This is sometimes useful because Trilinos solves a transformed matrix with a transformed right-hand side vector; this error estimate shows the error of the solution in a more natural basis. Also, sometimes Trilinos thinks that it hasn’t done a good enough job minimizing the residual (sometimes this happens, e.g., at t=0 when the potential is zero everywhere), and this usually clarifies whether the solve has actually succeeded or not.

For this to work, messagePotentialField must be true.

detailedAccuracyEstimate (optional boolean, default: false)

If false, a basic error estimate is printed out at steps determined by accuracyTestStepSlice. If true, the estimate will include a comma-separated-value table listing various norms and ways of measuring error (for each individual MPI rank, in the appropriate comms file, as well as over all ranks). These quantities described in terms of:

rho: the charge density (at a node), in C
phi: the potential (at a node), in V
laplPhi: the laplacian of the potential (at a node), in V per square m; in finite difference, the divergence involves the nodal volume V and dual face areas A; and the gradient uses edge lengths L, so we write laplPhi = - V^{-1} div A L^{-1} grad phi.
max_nn|phi A/Vdx|: for a given node, the maximum value, over the node and its nearest neighbors, of the magnitude of phi times A/(V L), where A, V, and L are the geometric factors involved in calculating the finite difference Laplacian.
avg|x|: the average value of the magnitude of x over all nodes
max|x|: the maximum value of the magnitude of x over all nodes

Norms and errors are reported for all interior grid nodes, as well as the subset of nodes that don’t touch any cell edges cut by a conducting boundary; norms and errors for individual ranks are written to the appropriate comms file along with the norms and errors for the entire simulation.

timingStepSlice (optional vector of integers, default
[10, 2147483647, 100]): a python-like ‘slice’ specifying the timesteps [start, stop, step] at which the timing for this updater will be printed out (search output for “implicitCylEmUpdaterUpdater timing” and “time per Trilinos matrix solve”).

printEdgeLengths (optional string, default 'none'): whether to print (selected) cell edge lengths to the per-rank comms file: can be ‘none’, ‘modifiedOnly’, ‘fractionalOnly’, or ‘all’. If ‘modifiedOnly’, only edge lengths that were modified (usually because they were extremely small; cf. samp:minEdgeFrac) will be printed.

printMatrices (optional vector of strings, default empty): a list of Trilinos matrices/vectors to write to file in Matrix Market format: can be ‘M’ (interior Laplacian), ‘sqrtLminV’, ‘sqrtLminOverV’, ‘PMPbar’ (laplacian operating on boundary values), ‘phi1V_0’, ‘phi1V_1’, …, ‘epsEdotdA_0’, ‘epsEdotdA_1’, …. where ‘_0’ and ‘_1’ (etc.) refer to the first and second (etc.) <Conductor>s or <ConductorShape>s. Note that in the Matrix Market format that Trilinos outputs, indices are one-based, although Trilinos uses zero-based indices.

rhsDumpSteps (optional vector of integers, default empty): a list of timesteps at which the right-hand side of the (transformed) Poisson equation should be written to disk.

indexField (optional string): The name of a vector <Field> into which to put the Trilinos vector index at each field location.

electricField (optional string): Name of electric field when dielectrics have conductivity. If this is used then conductivity must be specified in a DielectricShape block.

Example implicitCylEmUpdater Blocks

This updater sets the voltage on one conductor to a function of time, while the other conductor floats—with a prescribed current, as well as current from electrons and ions. Thus it uses <Conductor> blocks.

<FieldUpdater iemSolve>
  kind = implicitCylEmUpdater
  lowerBounds = [0 0 0]
  upperBounds = [885  833  13]
  readFields = [rhs]
  readComponents = [2]
  writeFields = [dE]
  writeComponents = [2]
  writeEquationToFile = False
  detailedAccuracyEstimate = 1
  accuracyTestStepSlice = [0   10000000   100000]
  timingStepSlice = [0   10000000   100000]
  verbosity = 255
  exactRestart = 1
  enforceMatrixSymmetry = 0
      electricField = E #(optional, in this case "E" must
                        # be a edge electric field that is updated by "dE")

  <DielectricShape DielectricWindowShape>
    kind = dielecShape
    gridBoundary = DielectricWindow
    permittivity = 3.7
    isFlipped = 0
        conductivity = 2.0 # optional block used to compute the
                           # electric field when the dielectric
                           # has a conductivity. "electricField"
                                               # must point to an edge electric field above.
  </DielectricShape>

  <Conductor lowerZ>
    kind = setVoltage
    calculateCharge = 0

    <ConductorShape shape>
      kind = slabCndctrShape
      lowerBounds = [0 0 0]
      upperBounds = [0 833 13]
      cellsInBoundsAreMetal = 1
    </ConductorShape>

    <Expression voltage>
      expression = 0.0
    </Expression>

  </Conductor>

  <Conductor upperZ>
    kind = setVoltage
    calculateCharge = 0

    <ConductorShape shape>
      kind = slabCndctrShape
      lowerBounds = [884 0 0]
      upperBounds = [884 833 13]
      cellsInBoundsAreMetal = 1
    </ConductorShape>

    <Expression voltage>
      expression = 0.0
    </Expression>

  </Conductor>

  <Conductor upperR>
    kind = setVoltage
    calculateCharge = 0

    <ConductorShape shape>
      kind = slabCndctrShape
      lowerBounds = [1 832 0]
      upperBounds = [884 832 13]
      cellsInBoundsAreMetal = 1
    </ConductorShape>

    <Expression voltage>
      expression = 0.0
    </Expression>

  </Conductor>

  <Conductor lowerR>
    kind = setVoltage
    calculateCharge = 0

    <ConductorShape shape>
      kind = slabCndctrShape
      lowerBounds = [1 0 0]
      upperBounds = [884 0 13]
      cellsInBoundsAreMetal = 1
    </ConductorShape>

    <Expression voltage>
      expression = 0.0
    </Expression>

  </Conductor>

  <LinearSolver linearSolver>
    kind = iterativeSolver
    verbosity = 255
    output = last
    maxIterations = 200
    tolerance = 1e-08
    convergenceMetric = r0

    <BaseSolver myBS>
      kind = gmres
      kspace = 30
      orthog = classic
    </BaseSolver>

    <Preconditioner precond>
      kind = multigrid
      mgOutput = 0
      mgDefaults = DD
      maxLevels = 10
      smootherType = symmetricGaussSeidel
      smootherSweeps = 3
      smootherPrePost = both
      coarseType = Jacobi
      dampingFactor = 1.333
      threshold = 0.0
      increaseDecrease = increasing
    </Preconditioner>

  </LinearSolver>

</FieldUpdater>