cutCellPoisson

cutCellPoisson

Works with a VSimPD license.

This is a MultiField FieldUpdater that solves the Poisson equation for the electrostatic potential, with Dirichlet/metal boundaries, using the Trilinos library to execute the large matrix solve. The solve can be performed for arbitrary potential specified on the boundaries (by a previous FieldUpdater), or the boundary values can be set to equipotentials with this updater. In the latter case, it is possible to specify that some of the distinct equipotential surfaces will float, with the potential determined by the charge on each floating conductor.

There are two very different ways to set boundary potentials. Most commonly, the user will specify a number of distinct conducting boundaries, along with the (time-dependent) voltage or charge on each boundary. Each conductor is an equipotential (as a conductor should be), and information about a conductors (such as voltage, charge, and capacitance) can be obtained using a conductorFunc.

Alternatively, the potentials on the boundaries can be taken from an arbitrary Field before each Poisson solve. Although the user must then figure out how to set those field values, this allows arbitrary (non-equipotential) Dirichlet boundary conditions. that specifies the potential at all necessary boundary points. This precludes functionality such as filling the potential within the conducting bodies, calculating the charge on conductors, and calculating capacitance and floating potentials.

cutCellPoisson Common Uses

The cutCellPoisson updater has a number of capabilities. Here are some common uses.

• Finding the electric potential in an arbitrary region with arbitrary Dirichlet boundary conditions specified via a Field.

• Finding the electric potential in a region bounded by arbitrarily-shaped conductors with time-varying voltages.

• Evolving the electric potential in a region with floating conductors, with voltages determined by absorbed/emitted (charged) particles and/or user-specified (external) currents.

• Calculating capacitances.

• Evolving electrostatic systems with simple external circuits and charged particle motion within the sytem. For example, one can observe the current induced by a charged particle flying through a grounded set of parallel plates. (However, currently, all capacitance is simulated spatially; there is no easy way to add external capacitances to the circuit.)

cutCellPoisson Basic Parameters

The cutCellPoisson 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: [phi, [edgeLengths]])

A list of Field objects containing the (scalar) electric potential (at grid nodes, units of Volts) to update, and, optionally, a 3-component vector field to which the grid cell edge lengths (i.e., the lengths not in metal conductors) will be written at initialization time. The edge lengths field may be helpful for diagnostics or post processing.

readFields (required string vector: [rho, [phi]])

A list of Field objects containing one or two elements, the (scalar) charge density (rho) field (at grid nodes, with units of C/m^3 in 3D, C/m^2 in 2D, or C/m in 1D) plus an optional (scalar) electric potential field (units of Volts).

If no second field is given, the boundary potentials will be specified by <Conductor> blocks. In this case, the <Conductor>s will be equipotentials, and their voltages or charges will be specified within the Conductor Block.

If the second field (phi) is given, then the boundary values will be taken from that field, with each boundary shape specified via a ConductorShape Block. This offers maximum flexibility in setting the Dirichlet boundary conditions, but because the (so-called) ‘conductor shapes’ are not necessarily equipotentials, information such as voltage, induced charge, and capacitance cannot be calculated. The second readField can be the same as the first writeField, or it can be a completely different Field.

<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.

cutCellPoisson Parameters (setting boundary values from the second readField)

<ConductorShape> code blocks

Each ConductorShape Block describes the the shape of a conductor. An arbitrary number of <ConductorShape> blocks can be used.

If <ConductorShape>s are specified, a second readField must be given; the potential on the <ConductorShape> boundaries will be set to the values from the second readField. Some other updater must be used to set the values of the second readField.

cutCellPoisson Parameters (setting equipotentials on Conductors)

<Conductor> code blocks

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

If <Conductor>s are specified, only one readField may be given.

cutCellPoisson 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.

cutCellPoisson 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.

cutCellPoisson Parameters (imposing a Neumann (gradient) condition on the simulation boundary)

setGradientBoundary = [xmin xmax ymin ymax zmin zmax]

The entries of the vector xmin, xmax, … and zmax represent the boundaries in the coordinate directions. The entries take either 1 or 0 to specify whether the boundary is Neumann or not. For example, if the upper boundary in the x-direction and the lower boundary in the z-direction are Neumann, the vector is [0 1 0 0 1 0]. The gradient values are specified using the vector ‘boundaryGradient’ as given below.

boundaryGradient = [0 0 0 0 0 0]

The vector entries as the gradient values to be set on the Neumann boundary. At present, the boundary condition takes only ZEROS.

enforceMatrixSymmetry = 0

The parameter specifies whether or not to enforce matrix symmetry. If two of the boundaries at a corner or edge are Neumann, according to the present discretization scheme, the matrix becomes asymmetric. Forcing the matrix to be symmetric will result in an incorrect solution at the corners or edges. ‘enforceMatrixSymmetry’ takes 1 or 0 to enforce or not. Note that the matrix is enforced to be symmetric by default.

cutCellPoisson Usage Notes

• To calculate the electric field from the potential (taking into account the conductor geometry), use the gradBndryUpdater.

• 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.)

cutCellPoisson 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 potential 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 \(- \nabla^2 \phi = \frac{\rho}{\epsilon_0}\) will be tested and printed out (search output for “cutCellPoissonSolver [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 “cutCellPoissonUpdater 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.

Example cutCellPoisson 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 poisson>
  kind = cutCellPoisson
  lowerBounds = [0 0 0]
  upperBounds = [40 80 40]
  writeFields = [phi]
  readFields = [rho]

# Best practice: do not give the same name to the <Conductor> and
# the gridBoundary.
  <Conductor sphereC>
    kind = setVoltage
    calculateCharge = true
    <ConductorShape shape>
      kind = gridBoundary
      gridBoundary = sphere
    </ConductorShape>
    <Expression voltage>
      expression = 10. * sin(628*t)
    </Expression>
  </Conductor>

  <Conductor smallSphereC>
    kind = floating
    initialCharge =  -6e-11
    ptclAbsorbers = [electrons.smallSphereAbsorber, ions.smallSphereAbsorber]
    ptclSources = [electrons.smallSphereAbsorber, ions.emitSmallSphere]
    <ConductorShape shape>
      kind = gridBoundary
      gridBoundary = smallSphere
    </ConductorShape>
    <Expression current>
      # expression should be a string, so if it looks like a number,
      # it should be surrounded by quotes.  This value is in Amps.
      expression = "-6e-11"
    </Expression>
  </Conductor>

  <LinearSolver mySolver>
    kind = iterativeSolver
    convergenceMetric = rhs
    <BaseSolver myBS>
      kind = cg
    </BaseSolver>
    <Preconditioner precond>
      kind = multigrid
      mgDefaults = SA
      numPdeEquations = 1
    </Preconditioner>
    tolerance = 1e-12
    maxIterations = 10000
  </LinearSolver>

</FieldUpdater>

This updater sets the boundary potential to the values given by the phi field (which must be previously set with some other updater…being very careful to set values at every boundary node). Thus it uses <ConductorShape> blocks, not <Conductor> blocks. It solves for the potential and writes it to the same phi field; however, in principle, different phi fields could have been used for reading and writing.

<FieldUpdater poisson>
  kind = cutCellPoisson
  lowerBounds = [0 0 0]
  upperBounds = [40 80 40]
  writeFields = [phi]
  readFields = [rho phi]

  <ConductorShape sphereC>
    kind = gridBoundary
    gridBoundary = sphere
  </ConductorShape>

  <ConductorShape smallSphereC>
    kind = gridBoundary
    gridBoundary = smallSphere
  </ConductorShape>

  <LinearSolver mySolver>
    kind = iterativeSolver
    convergenceMetric = rhs
    <BaseSolver myBS>
      kind = cg
    </BaseSolver>
    <Preconditioner precond>
      kind = multigrid
      mgDefaults = SA
      numPdeEquations = 1
    </Preconditioner>
    tolerance = 1e-12
    maxIterations = 10000
  </LinearSolver>
</FieldUpdater>