9 <li> \ref sect_ddr </li>
10 <li> \ref sect_elecpotfield </li>
11 <li> \ref sect_rescaling </li>
14 <li> \ref sect_bic </li>
18 <li> \ref sect_dde </li>
19 <li> \ref sect_dpe </li>
22 <li> \ref sect_timestepping </li>
28 <li> \ref sect_overview </li>
29 <li> \ref sect_size_refinement </li>
30 <li> \ref sect_criterion </li>
31 <li> \ref sect_curvature </li>
32 <li> \ref sect_fishpack </li>
33 <li> \ref sect_Conclusions </li>
39 <li> \ref sect_ARCoS_overview </li>
43 <li> \ref sect_sim </li>
44 <li> \ref sect_input </li>
45 <li> \ref sect_parameter </li>
46 <li> \ref sect_needle </li>
47 <li> \ref sect_output </li>
48 <li> \ref sect_structure </li>
49 <li> \ref sect_source </li>
56 \section sect_fluid_eq Fluid model
58 \subsection sect_physical Physical model
60 \subsubsection sect_ddr Drift-diffusion-reaction equations
76 In a fluid model of a streamer, we replace the individual particles in the
77 system by a density function \f${
n}(\mathbf{r},t)\f$.
78 The temporal evolution of
this density
function is governed by the physical
79 processes of the system and
this model takes the form of a
set of partial
80 differential equations (PDEs).
81 The derivation of
this so-called classical streamer model starts from the
82 <a
name=
"continuity"><
i>continuity</
i></a> equation. For particle species \f$i\f$, we have:
84 \frac{\partial n_i(\mathbf{r},t)}{\partial t} + \nabla \cdot \mathbf{j}_i(\mathbf{r},t) = S_i(\mathbf{r},t).
87 Here \f$S_i(\mathbf{r},t)\f$ represents the total of all sources and sinks of
88 species \f$i\f$. \f$\mathbf{j}_i(\mathbf{r},t)\f$ is the term
for the particle
89 current density of
species \f$i\f$.
90 Particles can drift and diffuse as described by the following
91 <a name=
"current_dens"><
i>expression</i></a>
for the particle current density
94 \mathbf{j}_i(\mathbf{r},t) = \mu_i n_i(\mathbf{r},t) \mathbf{E}(\mathbf{r},t) - D_i \nabla n_i (\mathbf{r},t).
97 In the above equation, the first term represents the particle drift due
98 to the electric field, with \f$\mu_i\f$ the mobility coefficient of
species \f$i\f$.
99 The second term represents the diffusion of particles due to the spatial gradient in particle densities with diffusion coefficient \f$D_i\f$. These equations can be derived from the Boltzmann equation
101 On the timescales involved, we consider only
electrons to be mobile,
while ions and neutrals remain stationary, which means that
for heavy
species, the [continuity equation](#continuity) is reduced to
103 \frac{\partial n_i(\mathbf{r},t)}{\partial t} = S_i(\mathbf{r},t).
106 The sources and sinks in the first equation play a very important role in the dynamics of the streamer. In
this model, the sources and sinks correspond to
reactions between the different charged and neutral species present in the gas. The source term due to a single reaction is the product of the densities of the species involved in the reaction and the field-dependent rate coefficient
for that reaction.
108 As an example, the <a name="reaction_imp_ion"><i>impact ionization reaction</i></a>
110 e^- + \mathbf{N}_2 \rightarrow 2 e^- + \mathbf{N}_2^+
114 S_{ionization} = k_{ion}(|\mathbf{E}|) n_e [\mathbf{N}_2].
116 Here \f$n_e\f$ is the local electron density, \f$[\mbox{N}_2]\f$ the density of
117 \f$\mbox{N}_2\f$ and \f$k_{ion}\f$ the reaction coefficient
for impact ionization
118 depending on the magnitude of the local electric field.
119 The value of \f$k\f$ can be determined in different ways, from experiments,
120 theoretical calculations or simulations.
121 The traditional approximation suggested by <a name=
"townsend_ionization"><
i>Townsend</
i> </a>
122 uses an empirical expression
for the impact ionization term, see the
123 [book by Y. P. Raizer,
"Gas Discharge Physics"][11]:
125 \frac{dn_e}{dt} = n_e \mu_e |\mathbf{E}| \alpha_0 e^{-E_0 / |\mathbf{E}|},
127 where \f$\mu_e\f$ is the electron mobility coefficient,
128 \f$\mathbf{E}\f$ is the local electric field and \f$\alpha_0\f$ and \f$E_0\f$ are
129 parameters that can be determined by fitting experimental data.
130 In gases that contain an electronegative admixture, such as \f$\mathbf{O}_2\f$, the process of
131 attachment can provide a sink
for the electron density through the following
134 e^- + \mathbf{O}_2 \rightarrow \mathbf{O} + \mathbf{O}^-
137 e^- + \mathbf{O}_2 + \mathbf{O}_2 \rightarrow \mathbf{O}_2^- + \mathbf{O}_2
139 The first attachment process is dissociative attachment, the second an example
140 of a 3-body attachment (a 3-body attachment can also occur with an oxygen and nitrogen
142 In the
case of the 3-body attachment, the reaction rate scales with the square of the oxygen density:
144 S_{3-body-att} = k_{3-body-att}(|\mathbf{E}|) n_e [\mathbf{O}_2]^2.
146 Further ionization losses can occur via one or more recombination processes,
147 but these typically have a timescale that is much longer than the timescale of
148 streamer development and propagation and are therefore primarily interesting
for
149 the evolution of the charge density after a streamer discharge,
150 as discussed in [PhD Thesis Wormeester(to appear in 2013)][6], Chapter 5.
154 In gases with attachment, detachment may occur, resulting in an additional
156 In gases that contain both nitrogen and oxygen, the photoionization process
157 provides a non-local source of
electrons.
158 Since photoionization is non-local, it can not be modelled by simple reaction
159 equations such as the ones
for impact ionization.
160 Instead, the local contribution of photoionization is calculated by spatially
161 integrating contributions from the entire domain.
162 The commonly used model
for photoionization and the approximations made to make
163 this model suitable
for simulation are discussed in
164 [PhD Thesis Wormeester(to appear in 2013)][6], Chapter 3, section Photoionization.
166 The reaction model
for streamer simulations can be very minimal or very extended,
167 with many species and
reactions, including metastables and various excited states.
168 The complexity of the reaction model depends on the purpose of the simulations.
169 For negative streamers in nitrogen, a model containing no more than three
species
170 (\f$\mbox{e}^-\f$, \f$\mbox{N}_2\f$ and \f$\mbox{N}_2^+\f$) and one reaction
172 is sufficient to simulate the dynamics of the streamer head, see [Montijn et al][1].
173 For more detailed studies of the streamer chemistry, the reaction model should be
174 as complete as possible.
177 \subsubsection sect_elecpotfield Electric potential and field
179 The streamer evolves under the influence of an electric field, which consists
180 of an externally applied electric field and the electric field generated by
182 These space charges are present at the head of the streamer as well as on the edge of the streamer channel.
183 For the further propagation of the streamer, the enhanced electric field in front of the streamer, generated by the space charge in the streamer head is essential.
184 We compute the net <a name=
"charge_dens"><
i>charge density</
i></a>
185 \f$q(\mathbf{r},t)\f$:
187 q(\mathbf{r},t) = e \sum_i q_i n_i(\mathbf{r},t),
189 where
for species \f$i\f$, \f$n_i\f$
190 denotes the density
function of these species and \f$q_i\f$ the charge of
191 a particle in units of the electron charge \f$e\f$.
192 From
this we compute the potential by solving the <a name="Poisson"><i>Poisson</i> equation
194 \nabla^2 \phi(\mathbf{r},t) = \frac{q(\mathbf{r},t)}{\epsilon_0}
196 and the electric field
198 \mathbf{E}(\mathbf{r},t) = -\nabla \phi(\mathbf{r},t).
201 \subsubsection sect_rescaling Rescaling to dimensionless units
203 The classical fluid model
for streamers can be rescaled to dimensionless units and
204 it is with these units that the code used in
this documentation works.
205 We refer the interested reader to the
206 [PhD thesis of Gideon Wormeester (to appear in 2013)][6].
207 From the [Townsend](#townsend_ionization) approximation
for ionization,
208 a characteristic field and length scale
209 emerges: \f$E_0\f$ and \f$l_0 = \alpha_0^{-1}\f$, respectively.
210 The characteristic velocity follows from the drift velocity of
electrons
211 in the characteristic field,
213 E_0: v_0 = \mu_e E_0.
216 The characteristic number density follows from the [Poisson](#Poisson) equation.
217 Values
for \f$\alpha_0\f$, \f$E_0\f$ and \f$\mu_e\f$ were obtained from
218 [PhD Thesis Montijn][7]
219 and are at standard temperature and pressure:
221 \alpha_0 & \simeq & 4332 \quad \mbox{cm}^{-1}\\
222 E_0 & \simeq & 2 \times 10^5 \quad \mbox{V} \mbox{cm}^{-1}\\
223 \mu_e & \simeq & 380 \quad \mbox{~cm}^2 \mbox{V}^{-1} \mbox{s}^{-1}.
226 When we insert these values in the characteristic scales, we obtain the values with which to rescale the equations:
228 l_0 & \simeq & 2.3 \times 10^{-4} \quad \mbox{~cm}\\
229 t_0 & \simeq & 3.0 \times 10^{-12} \quad \mbox{~s}\\
230 n_0 & \simeq & 4.7 \times 10^{14} \quad \mbox{~cm}^{-3}\\
231 D_0 & \simeq & 1.8 \times 10^{4} \quad \mbox{~cm}^2 \mbox{s}^{-1}.
233 We can now make the appropriate substitutions
234 (\f$t^d = t / t_0\f$ and similarly
for the other variables;
235 the superscript \f${\;}^d\f$ will be used to indicate that a variable is
236 in dimensionless form, where
this is not clear from the context.
237 For clarity of reading, the superscript \f${\;}^d\f$ will be omitted where it is clear that variables are
238 dimensionless) to obtain the classical fluid equations in
239 <a name=
"continuity_dimless"><
i>dimensionless continuity</i></a> form:
241 \partial_{t} + \nabla \cdot \mathbf{j}_i = S_i,
243 where \f$t\f$ is the dimensionless time, \f$\mathbf{j}_i\f$ the dimensionless
244 particle density current
for species \f$i\f$ and \f$S_i\f$ the dimensionless source term
246 \f$S_i\f$ is obtained by rewriting reaction equations such as the
247 [impact ionization reaction](#reaction_imp_ion) equation in dimensionless form,
248 where we remark that all rate-coefficients should also be rescaled.
249 The particle density current \f$\mathbf{j}_i\f$ is obtained by rescaling
250 [expression](#current_dens) into
251 <a name=
"current_dens_dimless"><
i>equation</i></a>
253 \mathbf{j}_i = -\mu_i n_i \mathbf{E} - D_i \nabla n_i,
255 where \f$\mathbf{E}\f$ is the dimensionless electric field and \f$n_i\f$,
256 \f$D_i\f$ and \f$\mu_i\f$ are the dimensionless particle density,
257 diffusion coefficient and mobility respectively of species \f$i\f$.
258 We find that in dimensionless units \f$\mu_i\f$ is equal to 1
while for heavy particles
259 \f$\mu_i\f$ is taken as 0, since heavy particles are assumed to be stationary in
261 The [dimensionless current density equation](#current_dens_dimless) can therefore be simplified to
263 \mathbf{j}_e = -n_e \mathbf{E} - D_e \nabla n_e
265 for electrons and \f$ \mathbf{j}_i = 0 \f$
for heavy particles.
266 The expression
for the [charge density equation](#charge_dens) \f$q\f$, is rescaled to
268 q(\mathbf{r},t) = \sum_i q_i n_i(\mathbf{r},t).
270 The [Poisson](#Poisson) equation is rescaled to
274 We remark that although the code
275 described here internally
276 works with the dimensionless equations and variables described in this section,
277 all results are presented in regular units unless otherwise noted.
278 Input parameters for the simulation code are expected to be in dimensionless units.
279 Finally we note that the rescaling to dimensionless units does not change
280 the structure of the equations, it is merely a rescaling to a different
281 set of units, where the dimensionless units yield a set of equations where some
282 constants (such as \f$e\f$, \f$\epsilon_0\f$, \f$\mu_e\f$) become unity.
284 \subsection sect_bic Boundary and initial conditions
287 We consider a cylindrical computational domain with coordinates:
289 (r,z,\theta) \in (0,
L_r) \times (0,
L_z) \times (0,2\pi).
291 Although the code described here is capable of performing full 3D calculations,
292 we assume cylindrical symmetry to greatly simplify the computations.
293 For any spatially dependent function \f$f(r,z,\theta)\f$, we assume:
294 \f$\partial_{\theta} f(r,z,\theta) = 0\f$.
295 Consequently, the coordinate system
for our computations is limited to
296 \f$(0,
L_r) \times (0,
L_z)\f$.
297 We consider a setup with a powered electrode at \f$z = L_z\f$ and a grounded
298 electrode at \f$z = 0\f$.
299 If the powered electrode is a plate, the following boundary conditions are used
for
300 the electric potential \f$\phi(r,z,t)\f$:
304 \forall z \; & \partial_r \phi(0,z,t) & = & 0\\
305 \forall r \; & \phi(
L_r,z,t) & = & 0\\
306 \forall z \; & \phi(r,0,t) & = & 0\\
307 \forall r \; & \phi(r,
L_z,t) & = & \phi_0
311 with \f$\phi_0\f$ the potential applied to the powered electrode.
312 If the powered electrode is a needle protruding from a plate, the needle has the
313 same potential \f$\phi_0\f$ as the plate.
315 <img src=
"http://md-wiki.project.cwi.nl/images/figure2.png" width=
"240px" />
316 \par Figure 1. Schematic of the computational setup.
318 In <a name=
"setup">Figure 1.</a>, the shaded rectangle represents the computational domain
for
319 the fluid equations, the thick horizontal lines the two planar electrodes
320 with the needle and its parameters depicted at the anode.
321 The area between the two planar electrodes is the computational domain
322 for the Poisson equation.
323 The needle is simulated by a single point charge, \f$Q\f$, chosen such
324 that \f$\phi =\phi_0\f$ in the point \f$P\f$, which is the tip of the needle.
325 The calculation assumes cylindrical symmetry around the needle axis
326 represented by the dashed-dotted line.
328 For the density equations, we use homogeneous Neumann conditions on all edges:
330 \partial_r
n(0,z,t) = \partial_r
n(
L_r,z,t) = \partial_z
n(r,0,t) = \partial_z
n(r,
L_z,t) = 0,
332 where we remark that
if the powered electrode is a needle, the computational domain
333 for the density equations is smaller than the computational domain
for the Poisson
334 equation and the \f$L_z\f$ values
for both domains are not equal.
335 This difference is a requirement of the numerical implementation of the needle
336 electrode and is further detailed in section \ref sect_needle.
338 While the boundary conditions mentioned above are the ones used by Wormeester,
339 the code that was used can also handle different choices of boundary conditions:
340 both homogeneous Neumann and homogeneous Dirichlet boundary conditions are available
341 for the top (\f$z =
L_z\f$), bottom (\f$z = 0\f$) and right (\f$r =
L_r\f$)
342 edges of the domain for both the densities and the potential.
343 The Neumann condition on the central axis of the cylindrical domain is required
344 for symmetry reasons.
346 As initial conditions for particle densities, two types of seeds are implemented in
347 the code. A homogeneous
seed, with a constant density over the entire domain and a
348 Gaussian seed of the form
350 n(r,z,0) = n_{max} \mbox{exp}(-\frac{r^2 + (z - z_0)^2}{\sigma^2}).
352 Here \f$z_0\f$ specifies the \f$z\f$-coordinate of the maximum of the
seed
353 (which is located on the symmetry axis with \f$r = 0\f$), where the density is
354 \f$n_{max}\f$. \f$\sigma\f$ is a measure of the radius of the
seed, it is the distance
355 at which the density drops to \f$e^{-1}\f$ of the maximum value.
357 In typical streamer simulations, a seed of
electrons and positive
ions is placed at the tip of the needle to initiate the discharge. Other than these Gaussian seeds and the neutral background gas, initial particle densities are zero with the possible exception of added background ionization, a homogeneous density of negative and positive
ions. The initial distribution of
electrons and
ions is charge neutral at every point of the domain.
359 \subsection sect_numeric Numerical method
361 The physical equations in section \ref sect_fluid_eq are to be solved numerically.
362 The computational code we have applied
for this uses finite volume methods to solve
363 a discretized version of the physical equations.
364 Here we give a basic summary of the numerical technique used.
365 For more details, the reader is referred to the work of [Montijn et al.][1],
366 upon which the current code is based.
368 \subsubsection sect_dde Discretization of density equations
370 The [dimensionless continuity ](#continuity_dimless) and the [dimensionless current density](#current_dens_dimless) equations
371 are discretized using finite volume methods and solved on a uniform rectangular grid with cells:
373 C_{ij} = [(
i - 1) \Delta r,
i \Delta r] \times [(j - 1) \Delta z, j \Delta z]\left(
i = 1 , \cdots , \frac{
L_r}{\Delta r}, j = 1 , \cdots , \frac{
L_z}{\Delta z}\right),
375 where \f$L_r\f$ and \f$L_z\f$ are the \f$r\f$- and \f$z\f$-dimensions of the grid and \f$\Delta r\f$ and \f$\Delta z\f$ the size of a cell in \f$r\f$- and \f$z\f$-direction, respectively. Particle density distributions are represented by their value in the cell center, which can be seen as an average over the cell. For some species \f$n\f$, we use \f$n_{
i,j}\f$ to denote the density at the cell center \f$C_{ij}\f$. For sake of clarity of notation we omit the superscript \f$^d\f$ indicating that variables are in dimensionless units.
377 The discretized continuity equations in cylindrical coordinates, with cylindrical symmetry (\f$\partial_{\theta} f = 0\f$) assumed, have the following form:
380 \frac{d n_{
i,j}}{d t} = & \frac{1}{r_i \Delta r} \Big(r_{
i - \frac{1}{2}} F^a_{
i - \frac{1}{2},j} - r_{
i + \frac{1}{2}} F^a_{
i + \frac{1}{2},j} + r_{
i - \frac{1}{2}} F^d_{
i - \frac{1}{2},j} - r_{
i + \frac{1}{2}} F^d_{
i + \frac{1}{2},j}\Big) + \\
381 & \frac{1}{\Delta z} \Big(F^a_{
i,j - \frac{1}{2}} - F^a_{
i,j + \frac{1}{2}} + F^d_{
i,j - \frac{1}{2}} - F^d_{
i,j + \frac{1}{2}}\Big) + S_{
i,j}.
384 Here \f$F^a\f$ and \f$F^d\f$ represent the advective and diffusive fluxes across the cell boundaries. Since we assume
ions and neutral particles to be stationary, these terms are nonzero only
for electrons. For heavy particles, only the source term \f$S_{ij}\f$ remains.
386 The advective flux, \f$F^a\f$ uses an upwind scheme with flux limiting and is defined as follows:
389 F^a_{
i + \frac{1}{2},j} = & E^+_{r; ~
i + \frac{1}{2},j} \Big[ n_{
i,j} + \psi(P_{
i,j})(n_{
i+1,j} - n_{
i,j}) \Big] \\
390 & E^-_{r; ~
i + \frac{1}{2},j} \Big[ n_{
i + 1,j} + \psi(\frac{1}{P_{
i+1,j}})(n_{
i,j} - n_{
i+1,j}) \Big]
396 F^a_{
i,j + \frac{1}{2}} = & E^+_{z; ~
i,j + \frac{1}{2}} \Big[ n_{
i,j} + \psi(Q_{
i,j})(n_{
i,j+1} - n_{
i,j}) \Big] \\
397 & E^-_{z; ~
i,j + \frac{1}{2}} \Big[ n_{
i,j + 1} + \psi(\frac{1}{Q_{
i,j+1}})(n_{
i,j} - n_{
i,j+1}) \Big],
400 where \f$E^+ = max(-E,0)\f$ and \f$E^- = min(-E,0)\f$ are used to distinguish the upwind directions
for the components of the electric field, \f$E_r\f$ and \f$E_z\f$, and we have
403 P_{
i,j} & = & \frac{n_{
i,j} - n_{
i-1,j}}{n_{
i+1,j} - n_{
i,j}}\\
404 Q_{
i,j} & = & \frac{n_{
i,j} - n_{
i,j-1}}{n_{
i,j+1} - n_{
i,j}}.
407 \f$\psi\f$ is the Koren limiter
function:
409 \psi(x) = max(0, min(1, \frac{1}{3} + \frac{x}{6}, x)).
411 The diffusive flux \f$F^d\f$ is calculated
using a second-order central differences scheme:
414 F^d_{
i + \frac{1}{2},j} & = & \frac{D}{\Delta r}(n_{
i,j} - n_{
i+1,j})\\
415 F^d_{
i,j + \frac{1}{2}} & = & \frac{D}{\Delta z}(n_{
i,j} - n_{
i,j+1})
418 and the reaction term \f$S_{
i,j}\f$ is computed as
420 S_{
i,j} = \sum_{A~\in~{reactions}} \Big[ k_A(|\mathbf{E}|_{
i,j}) \prod_{s~\in~{Spec(A)}} n_{s;
i,j} \Big]
422 where \f$k_A\f$ denotes the field-dependent reaction rate coefficient of
423 reaction \f$A\f$, and \f$Spec(A)\f$ the
set of species that appear as an input
424 for reaction \f$A\f$.
426 \subsubsection sect_dpe Discretization of the Poisson equation
428 We compute the net charge \f$q_{
i,j}\f$ in a cell center by adding up the contributions from the individual charged species:
430 q_{
i,j} = \sum_{s~\in~{species}} n_{s;
i,j} q_s.
432 With
this net charge, the electric potential \f$\phi\f$ can be computed in the cell centers through a second-order central approximation of the dimensionless Poisson equation:
434 q_{
i,j} = \frac{\phi_{
i+1,j} - 2 \phi_{
i,j} + \phi_{
i-1,j}}{\Delta r^2} + \frac{\phi_{
i+1,j} - \phi_{
i-1,j}}{2r_{
i,j} \Delta r} + \frac{\phi_{
i,j+1} - 2 \phi_{
i,j} + \phi_{
i,j-1}}{\Delta z^2}.
436 From the potential we can compute the components of the electric field from \f$\mathbf{E} = - \nabla \phi\f$ in the cell boundaries:
439 E_{r; ~
i + \frac{1}{2},j} & = & \frac{\phi_{
i,j} - \phi_{
i+1,j}}{\Delta r}\\
440 E_{z; ~
i,j + \frac{1}{2}} & = & \frac{\phi_{
i,j} - \phi_{
i,j+1}}{\Delta r}.
443 The electric field strength is determined at the cell center, so we have to compute the field components in the center by averaging the values on the boundaries after which we can compute the field strength:
445 |\mathbf{E}|_{
i,j} = \sqrt{\left(\frac{E_{r;
i - \frac{1}{2},j} + E_{r;
i + \frac{1}{2},j}}{2}\right)^2 + \left(\frac{E_{z;
i,j - \frac{1}{2}} + E_{z;
i,j + \frac{1}{2}}}{2}\right)^2}.
448 \subsection sect_timestepping Time stepping
450 The code uses the
explicit trapezoidal rule, a second order Runge-Kutta method,
451 for the temporal discretization with time step \f$\Delta t\f$.
452 Given some time step \f$t_i =
i \Delta t\f$, density distributions
453 \f$\mathbf{
n}_i(r,z) = \mathbf{
n}(r,z,t_i)\f$ and electric field
454 \f$\mathbf{E}_i(r,z) = \mathbf{E}(r,z,t_i)\f$, the densities and field at the next
455 time step, \f$t_{
i+1}\f$ are calculated by first computing an intermediate result
458 \overline{\mathbf{
n}}_{
i+1} = \mathbf{
n}_i + \Delta t F(\mathbf{
n}_i, \mathbf{E}_i).
460 Using these intermediate densities, the potential can be computed by solving the
461 Poisson equation, after which we obtain the intermediate electric field
462 \f$\overline{\mathbf{E}}_{
i+1}\f$.
463 With
this, we compute the
final values of the densities at \f$t_{
i+1}\f$:
465 \mathbf{
n}_{
i+1} = \mathbf{
n}_i + \frac{\Delta t}{2} F(\mathbf{
n}_i, \mathbf{E}_i) + \frac{\Delta t}{2} F(\overline{\mathbf{
n}}_{
i+1}, \overline{\mathbf{E}}_{
i+1}).
467 Finally, we again compute the potential and electric field, now
using the
final
468 values of the densities.
470 The size of the time step \f$\Delta t\f$ is determined by
using a Courant-Friederichs-Levy (CFL) restriction
for stability of the advection part of the equations:
472 \texttt{max} E_r \frac{\Delta t}{\Delta r} + \texttt{max} E_z \frac{\Delta_t}{\Delta_z} < \nu_a.
474 There are additional restrictions from other diffusion and reaction parts of the
475 equations, but they are dominated by the CFL criterior
for the advection part, see
477 The value of \f$\nu_a\f$ is typically
set to 0.25, which is well below the maximum
478 required
for stability.
479 We refer the interested reader to the [Monograph by Hundsdorfer and Verwer][8].
481 \section sect_refinement Overview of refinement strategies and criteria
483 \subsection sect_overview Overview
485 The \c ARCoS simulation code contains functions
for adaptive grid refinement (also known as adaptive mesh refinement or AMR). Since streamers span different length scales, there is a need to simulate relatively large physical domains
while still having high spatial resolution in areas such as the streamer head. To ensure that such large domains can be simulated without giving up resolution and accuracy, the numerical grid is refined adaptively at each time step. The equations are solved on a coarse grid, after which the solution is analyzed
using refinement criteria to determine the areas where refinement is needed. The equations are then solved on the refined subgrids after which the process is iterated. Grid generation and grid refinement are performed separately
for the density equations and
for the Poisson equation.
487 There are three
main refinement criteria. The first two concern refinement of the density grids: Refinement based on the absolute value of \f$\mathbf{E}\f$ and refinement based on the curvature of densities (both charge density and particle density). The grids used by the [FISHPACK][9] solver use their own refinement scheme where the decision to refine is made
if the difference between the solution on a grid and the solution on a finer grid exceeds a threshold. The [FISHPACK][9] solver is used
for both the Poisson equation that determines the electric potential of the system and the Helmholtz equations
for the photoionization reactions.
489 \subsection sect_size_refinement Size of the refined areas of the density grids
491 All CDR (Convection-Diffusion-Reaction, CDR is the shorthand term
for the density part of the code) refinement criteria are on a per-point basis, which means that the question whether to refine or not is initially answered
for every grid cell. This is inconvenient
for several reasons, primarily due to the computational cost of such a scheme. The regions containing the streamer head will almost always need to be refined, it is not necessary to evaluate
this point-by-point in these regions.
493 To ease
this problem, a minimal refinement area is defined by two parameters:\n
494 \c cdr\_brick\_dr and \c cdr\_brick\_dz, see e.g., file
default.cfg . The refinement module divides the grid it receives (
this can be the coarsest grid covering the entire domain or a refined grid covering only part of the domain, the code and grid structure are recursive) into
"bricks" of these dimensions and searches each brick
for cells that match the refinement criteria. Once such a cell is found, the entire brick containing that cell is refined.
496 For the [FISHPACK][9] module, a different approach is used. The refinement routine scans its input grid, starting at the top (\f$z = z_{min}\f$), going down per
"line" (a
set of cells with equal \f$z\f$d-coordinate). Once it finds a line with points that meet the refinement criterion it searches
for the first line that does not contain any points that meet the criterion. It then refines the smallest rectangular area that contains all the points that meet the criterion. This process is repeated until the bottom (\f$z = z_{max}\f$) of the grid is reached.
498 <img src=
"http://md-wiki.project.cwi.nl/images/refinement_cdr.png" width=
"360px" />
499 \par Figure 2. The nested structure of refined density grids
501 In <a name="refinement_cdr">Figure 2.</a>,
502 the black squares represent grid cells at the coarsest level (level 0),
503 the dark gray cells are the first refined sublevel (level 1). Two rectangular
504 grids are included at
this level, their shared border is indicated by the
506 The light gray cells show grids at a further refined level (level 2).
508 <img src=
"http://md-wiki.project.cwi.nl/images/refinement_poisson.png" width=
"360px" />
509 \par Figure 3. The nested structure of refined Poisson grids
511 The black grid, as shown in <a name=
"refinement_poisson">Figure 3.</a> is the coarsest level, the dark gray cells are the first
512 refined sublevel, the light gray cells show grids at a further refined level.
513 Each grid has at most one subgrid.
515 The tree of grids
for the density equations may contain refined grids that are adjacent to each other. A schematic showing the nested structure of refined density grids is shown in [Figure 2.](#refinement_cdr). The red line in
this figure indicates the shared border between two subgrids. For the Poisson-grids, such a structure is not possible and a grid can have at most one refined child-grid as depicted in [Figure 3.](#refinement_poisson).
517 \subsection sect_criterion The |E| criterion
519 The electric field criterion is the most simple of the three refinement criteria. It is an empirical criterion that is not directly motivated by the underlying numerics. A cell with coordinates \f$(r,z)\f$ qualifies
for refinement
if:
521 |\mathbf{E}(r,z)| > E_c
523 where \f$E_c\f$ is the threshold electric field strength
for refinement. \f$E_c\f$ is a user-determined parameter that is provided in the input file
for a run. Since
this criterion is independent of the grid level or the cell size, once a cell meets the criterion at the coarsest level, it will also
do so at every refined level. Because of
this property, the user can limit the refinement depth that is reached through
this criterion with the \c ref\_level\_eabs input parameter, see e.g., file
default.cfg. Setting \c ref\_level\_eabs to 1,
for example, restricts the refinement from the coarsest level to the first refined level due to the \f$|\mathbf{E}|\f$ criterion.
525 The \f$|\mathbf{E}|\f$ criterion is inflexible in the sense that it requires the user to have advance knowledge of what the field strengths will be. A possible alternative would be to replace the fixed threshold value \f$E_c\f$ by a dimensionless fraction \f$c\f$ and refine if
527 |\mathbf{E}(r,z)| > c E_{max}
529 with \f$E_{max}\f$ the maximum electric field strength in the computational domain. Since the electric field criterion is mostly empirical, picking the right value
for the refinement threshold may be a trial-and-error process.
531 \subsection sect_curvature The curvature criteria
533 There are two criteria that use the curvature of density functions in order to determine which areas to refine. If the curvature is large compared to the cell size, the numerics may become unreliable and it is desirable to work with a finer grid. For a density function \f$u(r,z)\f$ and a cell size \f$\triangle r \times \triangle z\f$ the curvature function \f$C_u(r,z)\f$ is a discretization of the second derivative of \f$u\f$ in cylindrical coordinates \f$(r,z)\f$:
536 C_u(r,z) = & \frac{1}{r + \frac{\triangle r}{2}}\Big[(r + \triangle r)\big(u(r + \triangle r, z) - u(r,z)\big) - r\big(u(r,z) - u(r - \triangle r,z)\big)\Big] + \\
537 & \big[u(r, z + \triangle z) - 2 u(r,z) + u(r, z - \triangle z)\big].
540 Rather than the absolute value of the curvature, the refinement module looks at the curvature relative to the global maximum, \f$Max(u)\f$. The
final criterion then reads:\n
542 \b Refine \f$(r,z)\f$ \b
if \f$\frac{C_u(r,z)}{Max(u)} > C_t\f$\n
544 with \f$C_t\f$ the threshold curvature. This refinement criterion is checked
for two density functions \f$u\f$. The first is the (absolute) charge density function. Here an extra condition applies: the absolute value of the charge needs to exceed a certain threshold value (which is hard-coded) before a cell can qualify for refinement based on this criterion. Secondly, the curvature criterion is applied to the particle density functions. Since only mobile particles require a high spatial resolution, any immobile species are not considered in these criteria (which currently excludes all species other than
electrons). The computational grids for these immobile species are simply the same as the grids used to solve the density equations for
electrons.
546 \subsection sect_fishpack FISHPACK refinement
548 The [FISHPACK][9] module, for the Poisson equation and the photoionization equations, uses a different set of grids than the CDR module and with it a different refinement scheme. Initially, two grids are set up, one coarse and one fine grid (with the fine grid having twice the spatial resolution in each dimension, so 4 times the number of cells). The Poisson/Helmholtz equation is then solved on both grids and the solution of the coarse grid is interpolated onto the fine grid. A grid cell then qualifies for refinement if the absolute difference between the interpolated coarse solution and the fine solution (this difference is called the error) is more than some user-defined threshold. When refinement is needed, a new set of grids is determined using the strategy mentioned earlier and the process is repeated until either the desired accuracy is reached or the maximum number of allowed refinement levels is reached. Since the [FISHPACK][9] module was originally only used to solve the Poisson equation for the electrostatic problem and the value of the electric field
549 is defined on the edge of a cell, a cell that does not meet the error-criterion still qualifies for refinement if its neighbor does meet the error-criterion.
551 One limitation to this scheme is the limited number of grid cells that the [FISHPACK][9] routine can handle. Since [FISHPACK][9] applies a cyclic reduction scheme, the round-off error increases with the number of grid cells. This places a limit on the size of grids that [FISHPACK][9] can solve. Once the refinement module wants to create a grid that is larger than the so-called [FISHPACK][9] limit, the refinement attempt is rejected and the code relaxes the error threshold by a factor of 2 and again determines the area to refine, using the new threshold.
553 To solve the photoionization problem, two Helmholtz equations need to be solved (For details on the implementation of photoionization, the reader is referred to [PhD Thesis Wormeester(to appear in 2013)][6], Chapter 3, section Photoionization
554 and references therein). Each of the so-called "photo-terms" has its own characteristic absorption length, which depends on the gas density and oxygen ratio. The term with the
short absorption length is often dominated by impact ionization in the head of the streamer, while the term with the
long absorption length is the
main contributor of
electrons in front of the streamer head that are required for a positive streamer to propagate.
556 The default behavior of the \c ARCoS code is to treat these two photoionization terms in the same manner as the Poisson problem when it comes to refinement: all user-definable parameters were equal. Since the term with the
short absorption length gives rise to a solution that benefits strongly from high spatial resolution (due to the steep gradients) it will easily trigger the refinement criterion. However, it is this term that is dominated by impact ionization, see [Luque et al.][2],
557 which reduces the relevance of accurate computation of this term. The user can therefore specify the refinement criteria for each of the two photoionization terms separately, providing the user with the means to allow the important,
long absorption length term to benefit from high spatial resolution, while reducing the computational cost incurred by the less important term. However, in tests it was found that tuning the refinement criteria for the photoionization terms has very little effect on
558 computational cost or results.
560 \subsection sect_Conclusions Conclusion
561 The adaptive refinement scheme of \c ARCoS allows for the simulation of large domains while maintaining high spatial resolution in regions that require this. A number of refinement parameters influence both the computational performance and the accuracy of the results, which means that the user has to monitor the results carefully. Since the refinement criteria were setup by
565 for simulations of air and pure nitrogen, application of the code to other gases may require changes to the values of the various thresholds used in the refinement criteria. An example is high-purity oxygen, with a small nitrogen admixture. In such a gas, ionizing photons will have a very
short characteristic absorption length and the calculation of the photoionization terms should be done with high accuracy close to the photon source, primarily the streamer head. However, the limitation of the [FISHPACK][9] refinement method does not permit several smaller, adjacent
566 refined sub-grids, which makes it difficult to properly focus on the streamer head without including too much of the channel.
568 \section sect_software ARCoS software
570 \subsection sect_ARCoS_overview Basic overview and functionality
572 The \c ARCoS simulation software was originally developed by Alejendro Luque as a more flexible version of the adaptive refinement code developed by Carolynne Montijn as described in [PhD Thesis Montijn][7].
573 The original code by Montijn has been written in Fortran90, while \c ARCoS has been written in C.
574 The original [FISHPACK][9] package used for solving the Poisson and Helmholtz equations is written
575 in Fortran77 and was developed by Adams, Swarztrauber and Sweet. The \c ARCoS code is now compiled
576 with [FISH90][10], a modernization of the original [FISHPACK][9], employing Fortran90 to slightly
577 simplify and standardize the interface to some of the routines.
579 \c ARCoS solves the fluid equations for streamers, described in section \ref sect_fluid_eq,
580 on nested Cartesian grids using an adaptive mesh refinement technique.
581 \c ARCoS allows for the simulation of both positive and negative streamers in the
582 electrode configurations plate-plate and needle-plate.
583 The needle-plate electrode geometry is included using a charge simulation method
585 This method replaces the electrode needle by a single point charge, with the
586 location and the size of the charge being updated at every time step to ensure the
587 potential at needle tip remains fixed at the predetermined value.
588 The limitation of this method is that the potential on the rest of the surface
589 of the simulated needle will not be accurate.
590 Consequently, the [continuity](
#continuity) equation is only solved on
591 a smaller grid, not containing the simulated needle.
593 The effect of
this is that \c ARCoS is not well suited
for the study of the
594 inception of streamers, as the area around the tip of the needle is not
596 However, since inception is often affected by the behavior
597 of individual particles, the use of a particle code such as described in
598 [Teunissen et al.][5] and [Li et al.][4].
599 is recommended
for studying streamer inception.
600 The purpose of the \c ARCoS code is to study streamer propagation in the phase
601 after the streamer has formed.
602 Studies performed by [Luque et al.][2] show that the dynamics of
603 streamers in later stages hardly depends on initial conditions.
605 \c ARCoS allows the user full control over the numerical parameters of the
606 simulation: grid size, refinement criteria and CFL numbers can be
set by the user.
607 The kinetic model,
i.e., the list of particle
species, their reactions and
608 initial densities as well as the diffusion and mobility coefficients can be
609 specified via a series of input files, allowing the user to fine-tune the properties
610 of the gas in which the streamer is simulated, see configuration file input/kinetic_example.cfg.
612 The \c ARCoS code can be downloaded from the website
615 \subsection sect_IO Handling the software, input and output
617 \subsubsection sect_sim Starting a simulation
619 Two input parameter files governs all details of the simulation:
621 \li \b Physical \b parameters such as voltage, electrode configuration, size of the gap, etc.
622 \li \b Numerical \b parameters such as grid size, refinement criteria, etc.
623 \li \b Practical \b details like the directory
name, where the output files should be stored and the interval at which output should be generated.
626 structured configuration files, is used
for reading, manipulating and writing these files.
627 The first file, stored as input/
default.cfg, must contain the
default values
for the global variables.
628 This file is a part of the streamer
package distribution.
629 The second file, say input/user_init.cfg, an example is given by input/example_user_init.cfg,
630 has a structure analogously to input/default.cfg, and should contain the parameters which differ
631 from the default values.
632 The program delivers a configuration file, say input/example_user_continue.cfg with the updated parameters
633 from input/user_init.cfg completed with the default values of input/default.cfg.
635 Since the execution time of a single run will take on the order of several days, it is recommended to
636 split the time period into smaller pieces, and restart the execution several times from the point
637 where the previous run stopped.
638 The easiest way to restart the simulation is
639 \li to copy the file input/user_continue.cfg into input/default.cfg, to be sure that equal values for the parameters are used,
640 \li to edit the file input/user_init.cfg, and change the \b \c t_end value, the \b \c restart value and the name of the \b \c load_file. See the end of of this section.
642 The use of the configuration files construction has the following advantages:
643 \li recompilation of the code is not necessary in case of a restart
644 \li the user always has a clear overview of the parameter values used
645 \li results of different or continuing runs can be stored in different output directories,
646 as listed in the configuration file.
647 \li besides the parameter value also comment coupled to a parameter can be changed in the configuration files written by the user. The length of the comment must be restricted by 100 characters
648 \li the order of the parameters in the configuration file is free
649 \li the user has the possibility to control the simulation, many parameter values can be changed.
651 It is easy to resume a simulation by using a set of output data as initial conditions.
652 One has to adapt the parameter file with modified start and end times for the simulation.
653 To start a \c ARCoS simulation use the following command from
654 the directory containing the executables:
657 ./arcos > out.example 2> err.example
661 The \c ARCoS program starts and it will print out the parameter values used:
662 \li in \c out.example
663 \li in input/user_continue.cfg.
665 The program will print some extra information to file \c out.example, e.g., the step size and
666 when a new set of output data has been written, and to which set of file names.
667 Warnings and errors will be printed in file \c err.example. The program can terminate in three different ways:
669 \li The preset end-time is reached.
670 \li The program is terminated by the user.
671 \li The time-step (as determined by the Courant criterion for stability, more details to come) has dropped below a preset threshold. This usually points to some form of instability.
673 In case the simulation runs on a PC or desktop machine, a convenient approach is to set a very
674 large value for the end time and, rather than having the program determine when to terminate,
675 keeping track of the progress of the streamer by checking the output files and manually terminating
676 the program when the desired output is reached, e.g., the streamer has reached the electrode, or,
677 it has started to branch.
678 In other cases, it may be necessary for the program to be able to run for a fixed amount of time.
679 For example when it needs to exit gracefully, which is required for profiling software to work.
680 Also in case of batch jobbing with a limited CPU wall clock time, like on most supercomputers,
681 the end time must be chosen corresponding to the wall clock time.
683 Data files with periodical data controlled by \c $output_dt, stored in directory \c $output\_dir,
684 have names using the format \c variable.C123abc.tsv:
685 \tparam <variable> is a particle species or electric field, e.g. \c charge, \c d_dummyminus, \c d_electrons \n
686 \tparam <123> is the sequence number of that particular output dataset and \n
687 \tparam <abc> specified the subgrid the output belongs to.
689 These files can be used to plot the solution or to restart the simulation. More details on this in the section on output files. To resume the example simulation from output set \c 123, call:
696 \subsubsection sect_input The input files
699 File and directory handling:
700 \li \b \c kin\_input - The filename of the [libconfig](http:
701 input file containing the species, seeds and reactions.
702 By convention, these files have the extension \b .cfg.
703 \li \b \c output\_dir - The directory where the output files will be stored.
704 This can be a relative path to the directory with the executable or an absolute
705 path. When the execution starts, the directory $output_dir must be present and writable.
706 \li \b \c output\_dt - The interval (in dimensionless units) with which an output dataset is to be saved.
707 Lower values mean more frequent output, which gives finer grained time-dependent
708 data at the cost of more disk space.
711 \li \b \c L\_r - Radius of the physical domain in dimensionless units.
712 \li \b \c L\_z - Length of the physical domain. Does not include the needle.
713 \li \b \c has\_photoionization - Whether to enable the photoionization module.\n
714 So, if \b \c has\_photoionization = 1, photoionization is present,\n
715 if \b \c has\_photoionization = 0, execution without photoionization.\n
716 \li \b \c photoionization\_file - Filename of the file containing the photoionization parameters.\n These parameters are further explained in appendix~ref{app:photoionization_parameters}.
718 \todo: WRITE APPENDIX
720 \li \b \c E0\_z, \b \c E0\_y, \b \c E0\_x - Components of the external electric field in dimensionless units. Since the electrodes are located at the top and bottom of the domain \f$(z = L_z\f$ and \f$z = 0)\f$, only \c E0\_z should be nonzero. Positive values of \c E0\_z will generate an electric field in the \f$^z\f$ direction, with the top electrode (at \f$z = L_z\f$) having a negative charge, generating negative streamers and vice verse.
721 \li \b \c pois\_inhom - Whether to use a needle-plane geometry or not.\n
722 If \b \c pois\_inhom = 1 then we have \b \c needle-plane case, \n
723 otherwise if \b \c pois\_inhom = 0 gives the plane-plane case. \n
724 See below for additional remarks regarding the needle-plane geometry.
725 \li \b \c needle\_length and \b \c needle\_radius - The length and radius of the needle in dimensionless units. Only applies when \b \c pois\_inhom = 1.
726 \li \b \c end\_t - Time at which the simulation will stop, in dimensionless units.
727 \li \b \c max\_ntheta - Number of azimuthal grid cells. Default \b \c max\_ntheta = 1. \n
728 \b \c max\_ntheta > 1 will activate the full 3D simulation (without cylindrical symmetry) using a pseudo-spectral method described in more detail in [Luque et al.(2011)][3].
730 Numerical parameters:
731 \li \b \c gridpoints\_r, \b \c gridpoints\_z - Number of gridpoints in \f$r\f$ and \f$z\f$ direction for the density equations at the coarsest level.
732 \li \b \c cdr\_max\_level - Maximum number of refinement levels of the grid for the fluid equations. \n
733 \b \c cdr\_max\_level = 0 means no refinement.
734 \li \b \c cdr\_bnd\_bottom, \b \c cdr\_bnd\_top, \b \c cdr\_bnd\_right - Boundary condition for the density equations at the bottom \f$(z = 0)\f$, top \f$(z = L_z)\f$ and right \f$(r = L_r)\f$ of the domain. \n
735 \b \c cdr\_bnd\_xxx = 1 gives homogeneous Neumann boundary conditions, \n
736 \b \c cdr\_bnd\_xxx = -1 gives homogeneous Dirichlet boundary conditions.
737 \li \b \c pois\_max\_level - Maximum number of refinement levels of the grid for the Poisson equation.
738 \b \c pois\_max\_level = 0 means no refinement.
739 \li \b \c pois\_bnd\_bottom, \b \c pois\_bnd\_top, \b \c pois\_bnd\_right - Boundary condition for the Poisson equation at the bottom \f$(z = 0)\f$, top \f$(z = L_z)\f$ and right \f$(r = L_r)\f$ of the domain. \n
740 \b \c pois\_bnd\_xxx = 1 gives homogeneous Neumann boundary conditions, \n
741 \b \c pois\_bnd\_xxx = -1 gives homogeneous Dirichlet.
743 \subsubsection sect_parameter The parameter file
745 In the parameter file input/default.cfg variables and their values are assigned with the following syntax
754 # <-- this starts a comment-line, which is ignored.
756 # White lines are also ignored.
757 # String-values should be between quotes, numeric values should not
759 # (Note: The above is an approximation)
760 # Scientific notation can be used:
761 pi_times_thousand=3.14E3
764 Non-existent or misspelled parameters are ignored.
765 Parameters are not required, every parameter has a
default value, which can be found
768 The following overview lists the parameters generally meant
for testing purposes and changing them is not required
for streamer simulations. Therefore they are best left at their
default value.
770 File and directory handling:
771 \li \b \c cdr\_output\_margin - Number of margin-cells to be added to the output of the density grids. 2 layers of ghost cells are added on the edge of the computational domain
for the purpose of enforcing boundary conditions. With
this parameter, these can be included in the output. Value must be smaller than or equal to 2. Default \b \c cdr\_output\_margin = 0.
772 \li \b \c pois\_output - Output the grids used in solving the Poisson equation? The grids used
for the Poisson equation are different than those used
for the density equations as detailed in section \ref sect_refinement. \b \c pois\_output = 0. Note that the absolute electric field is already saved as a density grid, so it is not required to
set this parameter to 1 to obtain
this data.
773 \li \b \c pois\_output\_margin - Number of margin-cells to be added to the output of the Poisson grids. Default \b \c pois\_output\_margin = 0. Note : \b \c pois\_output\_margin <= 2.
776 \li \b \c start\_t - Initial time. This parameter is used when resuming simulations and is automatically updated in that
case.
778 Numerical parameters: \
n
780 The adaptive mesh refinement criteria and the parameters related to them are discussed in more detail in section \ref sect_refinement.
782 \li \b \c ref\_threshold\_eabs - Refine grid
if \f$|\mathbf{E}|\f$ exceeds
this value.
783 \li \b \c ref\_level\_eabs - Maximum number of refinement levels of the grid
for the fluid equation due to the \f$|\mathbf{E}|\f$ criterion.
784 \li \b \c ref\_threshold\_charge - Refinement threshold
for the curvature of the charge density.
785 \li \b \c ref\_threshold\_dens - Refinement threshold
for the curvature of the particle densities.
786 \li \b \c cdr\_brick\_r, \b \c cdr\_brick\_z - Size \f$(r\f$ and \f$z)\f$ of the minimal refinement area. See Section \ref sect_size_refinement.
787 \li \b \c pois\_max\_error - An area of the Poisson grid is further refined
if the relative error between two consecutive refined levels exceeds
this value.
788 \li \b \c nu\_a - Courant number based on advection to determine the time step. Note, that \b \c nu\_a < 1 to satisfy CFL stability. In streamer simulations,
this time step restriction will dominate over other parameters (\b \c nu\_d (diffusion) and \b \c nu\_rt (reaction)). More details on the time stepping can be found in \ref sect_timestepping.
790 \subsubsection sect_needle Implementation of the needle-plane configuration
792 As mentioned in section \ref sect_overview , the needle-plane electrode geometry is
793 implemented
using a charge simulation technique, which means that the entire needle is
794 represented by a single point charge located on the axis.
795 The position and strength of
this charge is updated each time step to ensure that the
796 potential at the point that would be the tip of the needle remains fixed.
797 A schematic depiction of
this setup can be seen in [Figure 1.](#setup).
798 This provides a reasonable approximation of the potential in the area below
799 the needle (the needle is always located at the top of the domain), but the potential
800 will be wrong in the areas to the sides of the needle.
801 Consequently, the density equations are only solved from the tip of the needle and downwards.
802 This means that the computational domain
for the density equations is smaller than
803 that
for the Poisson equation.
805 In the specification of the external electric field, \c E0\_z is interpreted as the electric field between the two planar electrodes, far away from the needle. The applied potential is computed as follows:
807 V = E_{0,z} * (
L_z + L_{needle})
810 \subsubsection sect_output Output
812 The code generates a large amount of output files at every \c \b output\_dt units
814 These output files combined are sufficient to resume the simulation
815 from that time step and can also be used
for data analysis.
816 Output files are named \c variable.C123abc.tsv, where
'tsv' stands
for 'Tab-Separated Values'.
817 \c variable is a particle species or electric field.
818 Examples include \c
electrons, \c n2plus and \c eabs
819 (the latter being \f$|\mathbf{E}|\f$). \c 123 is the sequence number or
820 output time step (which is unrelated to the time step of the numerical scheme) of that
821 output
set, it starts at \c 000
for the first
set.
822 This is the initial condition of the system before the simulation starts at \f$t = 0\f$.
823 The alphabetic part of the
filename after the output time step, \c abc denotes the
824 subgrid contained in that file and
this extension is defined recursively.
825 The coarsest grid covering the entire domain has no such alphabetical extension
826 (example: \c electrons.C123.tsv).
827 At every next level, the subgrids at that level are named \c a, \c b, \c c, ...
828 and
this letter is appended to the alphabetical extension.
829 The file \c electrons.C123ba.tsv contains electron densities at the
830 123-th time step (so at \f$t = 123 * output\_dt\f$ )
for the first
831 subgrid of the second subgrid of the
main grid.
832 Beware that since the grid refinement is adaptive, each time step has different subgrids.
834 <img src=
"http://md-wiki.project.cwi.nl/images/grids.png" width=
"360px" />
835 \par Figure 4. Schematic depiction of the naming convention
for output files. Each next letter corresponds to a
new refined level.
837 Data, as shown in <a name=
"output_grids">Figure 4.</a>, is stored as plain text,
838 with each line containing a single number.
839 Data is ordered in columns (with fixed \f$r\f$ coordinate).
840 So to read the data, use the following pseudo-code:
842 for (
i = 0,
i < rmax*zmax;
i++)
846 data[r,z] = read_line_from_file();
850 The dimensions of the grid are not contained in the data files.
851 Instead,
for each subgrid and output-step, two additional files are created:
854 2. _z.C123abc.tsv_\n
\n
856 corresponding to subgrid _abc_ of output step _123_. The structure of these files is identical to that of the regular data files, but instead of particle densities or electric field strengths, these files contain the \f$r\f$ and \f$z\f$ coordinates of the center of the cell corresponding to that line-number. So to determine the coordinates of the \f$n^{th}\f$ line in a regular data file, simply read the \f$n^{th}\f$ line of the corresponding \b r and \b z files.
858 \subsubsection sect_structure File structure
860 <img src=
"http://md-wiki.project.cwi.nl/images/structure.png" width=
"500px" />
861 \par Figure 5. File structure of of ARCoS package.
863 The distribution of the ARCoS
package contains several directories:
865 \li \b \c FISH90 - the directory where FISH90 or FISHPACK should be downloaded
867 \li \b \c output - this directory may be empty. Its name must correspond to the value of \b \c output_dir in file input/default.cfg. Or, the value in input/user_init.cfg if present.
868 \li \b \c python - contains plot files. By means of \b \c plotvar pictures can be made of the output files. Not yet implemented.
869 \li \b \c doc - directory with files for [Doxygen][12] to generate documentation from source code and this \b \c MANUAL. The source of this manual (in MANUAL.md) can be found in its directory \b \c markdown. See also the \b \c doxygen_config_file.
870 \li \b \c arcos_f90 - a part of the \b \c functions in file cdr.c have been replaced by a piece of \b \c FORTRAN90 code in order to accelerate the simulation.
871 \li \b \c src - this directory contains the source files written in \b \c c.
872 \li \b \c include - this directory contains the include files
873 \li \b \c input - this directory contains all input files. Most of them are libconfig configuration files.
876 \subsubsection sect_source Source files
878 The \c ARCoS simulation software was mostly written in \c C and its source code is split up in several files, each dealing with a separate part of the program. Most source files have an associated header file (the source file \b \c example.c has header file \b \c example.h ) containing the type-definitions and preprocessor macros. The function prototypes are aggregated in the header file \c proto.h. Below is a short summary of the important source files and the functionality that is contained within them.
880 \li cdr.c - Functions for solving the convection-diffusion-reaction (CDR) equations, creation, manipulation and refinement of CDR grids and time stepping.
881 \li configuration.c - Module for input/output of parameters. The code uses the library
883 \li cstream.c - Contains some general initialization and termination functions.
884 \li dft.c - Functions related with discrete fourier transformations.
885 \li grid.c - Low-level functions for handling of grids, both CDR and Poisson grids.
886 \li interpol2.c - Interpolation functions for the mapping of one grid to another (for example during refinement).
887 \li main.c - Functions for reading input parameters, starting of the code and the main loop.
888 \li mapper.c - Mapping of one grid tree onto another.
889 \li misc.c - Miscellaneous utilities for allocating memory.
890 \li photo.c - Photoionization functions.
891 \li poisson.c - Functions for solving the Poisson equation, including manipulation of Poisson grids and calling the external [FISPACK][9] solver.
892 \li reaction.c - Functions for computation of reactions between species as part of the density equations.
893 \li react_table.c - Performs initialization of reaction coefficient tables as well as table lookups.
894 \li rt.c - Functions for handling the loading of the input file containing the kinetic model (species, reactions, seeds).
895 \li rz_array.c - Low-level functions for handling Fortran-compatible arrays.
896 \li sprites.c - Routines for the sprites module.
1.8.2
