|
MPI-AMRVAC
3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
|
The test particle module in MPI-AMRVAC provides the possibility to track test particles, i.e. particles evolved according to the HD/MHD fields in a simulation without any feedback from the particles to these fields. Test particles can be added on top of an (M)HD state and evolved concurrently to the fluid. Alternatively, test particles can be evolved in a static snapshot, i.e. without evolving the underlying fluid quantities. The test particle module can also be used to sample the fluid quantities at specific locations (which may differ from the computational grid points) and output the results separately from the usual .dat and .vtu files.
The particle module can work in four different modes, depending on the user's choice:
In the next sections, the details of each mode are illustrated.
This mode is used to track fluid parcels moving through the simulation domain with the fluid speed. At each step, a simple equation of motion is solved for each particle,
\[ \frac{d\mathbf{x}}{dt} = \mathbf{v}, \]
where \(\mathbf{v}\) is the local fluid velocity, linearly interpolated at the particle position \(\mathbf{x}\). The equation of motion above is solved by means of a fourth-order Runge-Kutta integrator with adaptive time-stepping.
As the Lagrangian motion of the fluid particles is computed, it is possible to track various fluid quantities at the particle locations, e.g. the fluid density, pressure, etc. See the "Payloads" section below for further details.
A charged particle of charge \(q\) and mass \(m\) in electromagnetic fields evolves in time according to Newton's equations of motion
\[ \frac{d\mathbf{x}}{dt} = \mathbf{v}, \]
\[ \frac{d\mathbf{v}}{dt} = \frac{q}{m}\left(\mathbf{E}+\mathbf{v}\times\mathbf{B}\right), \]
where \(\mathbf{x}\) is the particle position, \(\mathbf{v}\) is the velocity, and the electric and magnetic fields \(\mathbf{E}\) and \(\mathbf{B}\) provide the accelerating Lorentz force for the particle. The typical particle trajectory is a superposition of a parallel motion along magnetic field lines, a gyro-motion around magnetic field lines, and various "drift" mechanisms that allow the particle to cross magnetic field lines. As the magnetic field does no work on charged particles, a change in the particle kinetic energy is only associated to the presence of electric fields.
The equations above are valid for a particle travelling at speeds much smaller than the speed of light in vacuum \(c\). When \(v\rightarrow c\), it is necessary to adopt the relativistic equations of motion
\[ \frac{d\mathbf{x}}{dt} = \frac{\mathbf{u}}{\gamma}, \]
\[ \frac{d\mathbf{u}}{dt} = \frac{q}{m}\left(\mathbf{E}+\frac{\mathbf{u}}{\gamma}\times\mathbf{B}\right), \]
where \(\mathbf{u}=\gamma\mathbf{v}\) is the normalised particle momentum, with \(\gamma=\sqrt{1+u^2/c^2}=1/\sqrt{1-v^2/c^2}\) the Lorentz factor.
For the relativistic case, the solution of the equations for the Lorentz dynamics can be carried out numerically in several fashions. Below we list the options available in MPI-AMRVAC:
In this mode, a "reduced" set of equations is employed to resolve the dynamics of charged particles, according to the so-called "guiding centre approximation" (GCA) paradigm. This neglects the particle gyro-motion, and may present advantages in certain physical scenarios. The equations of motion in this case are
\[ \begin{aligned} \frac{d\mathbf{R}}{dt} = & \mathbf{v}_{\|} + \mathbf{v}_E + \frac{\hat{\mathbf{b}}}{B}\times\frac{\mu}{q}\nabla B \\ & + \frac{m\hat{\mathbf{b}}}{qB} \times\left\{ v_{\|}^2\left(\hat{\mathbf{b}}\cdot\nabla\right)\hat{\mathbf{b}} + v_{\|}\left(\mathbf{v}_E\cdot\nabla\right)\hat{\mathbf{b}} + v_{\|}\left(\hat{\mathbf{b}}\cdot\nabla\right)\mathbf{v}_E + \left(\mathbf{v}_E\cdot\nabla\right)\mathbf{v}_E \right\}, \end{aligned} \]
\[ \frac{dv_{\|}}{dt} = \frac{q}{m}E_{\|} - \frac{\mu}{m}\hat{\mathbf{b}}\cdot\nabla B + \mathbf{v}_E\cdot\left(v_{\|}\left(\hat{\mathbf{b}}\cdot\nabla\right)\hat{\mathbf{b}}+\left(\mathbf{v}_E\cdot\nabla\right)\hat{\mathbf{b}}\right), \]
\[ \frac{d\mu}{dt} = 0, \]
where \(\mathbf{R}\) is the guiding centre position and \(\mathbf{v}_{\|}=(\mathbf{v}\cdot\hat{\mathbf{b}})\hat{\mathbf{b}}\) is the particle velocity in the direction parallel to the magnetic field. The unit vector \(\hat{\mathbf{b}}=\mathbf{B}/B\) also defines the drift velocity \(\mathbf{v}_E=\mathbf{E}\times\hat{\mathbf{b}}/B\) and the parallel electric field \(E_{\|}=\mathbf{E}\cdot\hat{\mathbf{b}}\). The conservation of the adiabatic invariant \(\mu=mv_\perp^2/(2B)\), with \(v_\perp\) the particle velocity perpendicular to \(\mathbf{B}\), is an assumption of this paradigm and may not be valid in general.
For particles travelling at velocities close to \(c\), the special-relativistic formulation of the equations above reads
\[ \begin{aligned} \frac{d\mathbf{R}}{dt} = & \frac{\mathbf{u}_{\|}}{\gamma} + \mathbf{v}_E + \frac{\kappa^2\hat{\mathbf{b}}}{B}\times\left\{\frac{\mu_r}{q\gamma}\nabla (B/\kappa) + \frac{u_{\|}}{\gamma}E_{\|}\mathbf{v}_E\right\} \\ & + \frac{m\kappa^2\hat{\mathbf{b}}}{qB}\times\left\{\frac{u_{\|}^2}{\gamma}\left(\hat{\mathbf{b}}\cdot\nabla\right)\hat{\mathbf{b}} + u_{\|}\left(\mathbf{v}_E\cdot\nabla\right)\hat{\mathbf{b}} + u_{\|}\left(\hat{\mathbf{b}}\cdot\nabla\right)\mathbf{v}_E + \gamma\left(\mathbf{v}_E\cdot\nabla\right)\mathbf{v}_E\right\}, \end{aligned} \]
\[ \frac{du_{\|}}{dt} = \frac{q}{m}E_{\|} - \frac{\mu_r}{m\gamma}\hat{\mathbf{b}}\cdot\nabla \left( B/\kappa\right) + \mathbf{v}_E\cdot\left(u_{\|}\left(\hat{\mathbf{b}}\cdot\nabla\right)\hat{\mathbf{b}}+\gamma\left(\mathbf{v}_E\cdot\nabla\right)\hat{\mathbf{b}}\right), \]
\[ \frac{d\mu_r}{dt} = 0, \]
where \(\mathbf{u}_{\|}=\mathbf{v}_{\|}\gamma\), and \(\kappa=1/\sqrt{1-v_E^2/c^2}\) is the Lorentz factor corresponding to the frame moving at a speed equal to \(\mathbf{v}_E\). These equations are obtained by averaging over the particle gyro-motion, and therefore assume the conservation of the relativistic magnetic moment \(\mu_r=m\gamma^2 v_\perp^2/(2B)\). Because of the averaging, it can be shown that the expansion of the guiding centre velocity leads to the definition of the Lorentz factor \(\gamma = \kappa\sqrt{1+(u_{\|}^2+2\mu_r B/m)/c^2} = \kappa\sqrt{(1+2\mu_r B/(mc^2))/(1-\kappa^2 v_{\|}^2/c^2)}\) up to first order. This implies that the dominant drift mechanism in the guiding centre motion is represented by \(\mathbf{v}_E\). For both the Newtonian and relativistic cases, the assumption of slowly time-varying electromagnetic fields has been made, in order to simplify the equations of motion. This is justified by the reasonable assumption that the particle dynamics takes place on much faster time scales than that of the typical MHD evolution.
The solution of the equations above is carried out in MPI-AMRVAC with a fourth-order Runge-Kutta numerical integrator with adaptive time stepping.
With this mode, the user is allowed to sample the fluid quantities at an arbitrary number of points in space, which may or may not be spatially distinguished from the numerical grid points. All quantities are retrieved at such locations via linear interpolation. By default, the quantities used in a standard HD/MHD run are computed at the sampling points (density, momentum, pressure/energy/temperature, and magnetic field components in the MHD case). The user is also free to define their own list of additional quantities to be sampled (e.g. the current).
All particle-related computations can be activated by setting hd_particles=.true. in the hd_list (if running in HD) or mhd_particles=.true. in the mhd_list (if running in MHD) of the .par file. If the particle calculations are switched on, additional parameters can be specified in the particles_list. Below is a description of these parameters and their role in the particle modules.
Below, we describe the essential steps needed to correctly set up a particle simulation.
physics_type_particles, which can be set to 'advect', 'Lorentz', 'GCA', or 'sample'.num_particles (an integer).Particle position and velocity: for all choices of physics_type_particles, the particles must be initialised at the desired locations in space. For the 'Lorentz' and 'GCA' modes, a velocity distribution should also be specified.
The user can specify their own particle initialisation setup in the mod_usr.t file. Then, the user must associate the pointer usr_create_particles with their particle initialisation subroutine (named e.g. generate_particles). The format for such a routine must be:
Here, x is the particle position, v is the particle velocity, q is the particle charge, and m is the particle mass. The follow variable tells the routine whether a certain particle should be tracked individually during the simulation (see section "Output and visualisation" below). While the position will be used in all modes, information on the velocity, charge, and mass will only be employed if physics_type_particles is 'Lorentz' or 'GCA'. If the user does not provide their own particle initialisation subroutine (and/or if the latter is not associated with the usr_create_particles pointer), the code will by default initialise all particles at randomly generated locations inside the domain, with zero velocity, mass, and charge, and with follow=.false..
By default, the code generates the particles on the processor #0 and then distributes them among the processors depending on which processor handles which spatial region. After the particles generated with the initial position and velocity have been distributed, the user has the option to perform additional operations such as modifying the particle velocity based on local MHD properties, discard particles that are found outside of specific regions of interest, etc. For example, the user may choose to retain only particles belonging to regions wheret he temperature is above a chosen threshold (an operation that cannot be performed at generation time, since processor #0 where particles are generated does not possess information on the spatial regions where the particles will be sent). This can be done by associating the pointer usr_check_particle in the user file. The pointer must be associated with a subroutine (named e.g. particle_modification) which must have the follwing format:
Here, the user is free to apply modifications to the particle velocity, mass, charge, and to the follow parameter, but not to the particle position, which is assumed will stay unchanged (otherwise, handling the particles may require further communications). Whenever this subroutine returns a check=.false. flag, the particle will be discarded. For all particles that should be kept, it is necessary that check=.true.
Payloads: particle simulations are especially flexible in terms of the quantities that can be dynamically stored in the particle output files. On top of tracking positions and velocities, an arbitrary number of payloads can be assigned to each particle in order to monitor additional physical aspects. As an example, in the 'advect' mode each particle can be assigned to track the local fluid density, which will be then stored in a payload variable and added to the output. A number of default payloads can be calculated and stored for each running mode. Additionally, the user can define custom payloads. The number of default and custom payloads is chosen by the user by setting the parameters ndefpayload and \nusrpayload in the particles_list of the .par file. By default, ndefpayload=1, and payload tracking can be suppressed by setting ndefpayload=0.
The default payloads, depending on the running mode, are:
'advect' mode, the fluid density at the particle location will be tracked and stored in the first payload.'Lorentz' mode, up to four payloads can be updated by default: the particle Lorentz factor (=1 if relativistic=.false.), the particle gyroradius, the magnetic moment, and the local value of \( \textbf{E}\cdot\textbf{B}\).'GCA' mode, there are 14 default payloads:'sample' mode, by default (regardless of the value of ndefpayload in the .par file) there will be a number of payloads n=nw, where nw is the number of variables in the fluid simulation. Each of these payloads samples one of the primitive fluid quantities, and therefore in the .csv output these payloads are named according to the names given to the primitive quantities.A custom payload update routine allows the user to store additional payloads (on top of the default ones). This can be done in the mod_usr.t file via a user-defined routine which must be associated with the usr_update_payload pointer at the beginning of mod_usr.t. The required format for a user-defined payload update routine (e.g. named update_payload) is:
Particle calculations in MPI-AMRVAC can be carried out in two different modes, namely i) concurrently to the (M)HD evolution, or ii) in a fixed (M)HD snapshot. Additionally, the convert operations that can be performed with (M)HD output files also affect the particle outputs. Each of these options is illustrated below.
mod_usr.t file and all parameters for the (M)HD simulation are set in the .par file, the particle module can be activated simply by including the statement hd_particles=.true. or mhd_particles=.true. in the mhd_list of the .par file. This will tell the code to perform the particle calculations according to the parameters specified in the particles_list block of the input file. The (M)HD calculation will behave as usual, and the particle results will be stored in the output according to dedicated parameters (see section "Output and visualisation" below).convert functionality of the code. First, the user must set the parameters convert=.true., autoconvert=.false. in the filelist block of the .par file. Finally, the particle module must be activated via hd_particles=.true. or mhd_particles=.true. in the mhd_list. The particles will be initialised according to the user-defined (or the default) routines. Alternatively, an initial particle snapshot can be used as initial condition for the particles, if provided in the same folder of the fluid snapshot (see the "output and visualisation" section below for more info on particle outputs). In this mode, the particle integration will be carried over until time_max is reached; the code will assume that the initial time is the same as that stored inside the fluid snapshot. This implies that, if the snapshot was saved at time (say) t=9 and the user wishes to integrate particles in this snapshot for a total time of 1, they will have to specify time_max=10 in the .par file.restart_from_file parameter in the .par file for the fluid. If hd_particles=.true. or mhd_particles=.true. in the same .par file, the code will search for a particle snapshot (.dat file) in the same directory of the fluid snapshot used for the restart (and with the same base file name - see next sections for more info on particle output). The particles will then be initialised according to the information stored in the particle snapshot. Note that if the particle snapshot is not found, then the code will restart the fluid calculations from the given fluid snapshot, and it will then initialise the particle module with the user-provided particle initialisation routines (or the default routines if the user-defined ones are not provided). In essence, this allows for either restarting a previously interrupted particle run or for starting a new particle run on top of a restarted fluid run.Below is a description of the various outputs associated to the use of the particle module, together with a brief introduction to the visualisation of the particle results.
Output files and formats: the output of particle calculations (regardless of the chosen mode) consists of two file types, which will be stored in the same folder as the (M)HD results.
Particle snapshots (format .dat) will be produced only when running the particle module along with a time-evolving (M)HD calculation, and will be produced with the same cadence of the (M)HD output. These files contain all particle information in raw binary format, and do not have any direct use beside providing a checkpoint for restarts. The base file name for particle snapshots will be constructed by concatenating the base_filename given in the .par file with an additional _particles and followed by the output number (same as the fluid .dat snapshots). Note that particle snapshots will NOT be produced when running the particle module in a fixed (M)HD snapshot.
Particle individual and ensemble outputs (format .csv) can be used to easily analyse particle results. Whether or not these files are produced is controlled by the parameters write_individual and write_ensemble in the particles_list of the .par file. The output cadence to both individual and ensemble files is controlled via the dtsave_particles parameter in the same list. The output cadence may or may not be the same of the (M)HD output; the two can be specified completely independently.
The standard output quantities stored in the .csv files are, for each particle:
dt;(x1,x2,x3);(u1,u2,u3);pl01, pl02, ..., plN (where N is the number of default payloads specified via ndefpayload);usrpl01, usrpl02, ..., usrplN (where N is the number of custom payloads specified via nusrpayload).Note that when relativistic=.true. in the 'Lorentz' mode, the particle velocity will be replaced by the particle normalised momentum in the output files. In the 'GCA' mode, the quantities stored in (u1,u2,u3) are not the full particle velocity components, but rather the particle parallel velocity (or normalised parallel momentum) in u1 and the magnetic moment (or relativistic magnetic moment) in u2. The Lorentz factor will be stored in u3 if relativistic=.true., otherwise u3=1 will be set by default.
If write_individual=.true., the code will produce one .csv file for each particle that the user flagged with follow=.true. at initialisation (by default, when the user does not provide a custom initial particle setup, all particles are flagged with follow=.false.). The base file name for individual particle .csv's is obtained by concatenating base_filename with _particle_XXXXXX, where XXXXXX is a unique integer index (in %06d format) associated with each particle at initialisation. All information for a single individually followed particle will be stored over time, with cadence equal to dtsave_particles, inside the same .csv file, such that at the end of the run that file will contain the complete history of that single particle. Therefore, individual particle .csv's can be used to visualise the trajectory and time evolution of the quantities associated to specific particles.
If write_ensemble=.true., the code will produce ensemble .csv files which, oppositely to individual .csv files, gather information from all particles at a specific time. A single ensemble .csv file will be produced at each output time specified via dtsave_particles. The base file name for ensemble .csv's is obtained by concatenating base_filename with _ensemble_XXXXXX, where XXXXXX is a number (format %06d) corresponding to the n-th output time, based on the cadence specified via dtsave_particles. Ensemble .csv's are useful to visualise all particles together at a specific moment in time.
As a practical example, suppose the user has chosen dtsave_particles such that 10 particle outputs will be produced during a simulation. Suppose further that both write_individual=.true. and write_ensemble=.true., nparticles=100, and the user has flagged particles with index 36, 47, and 99 with follow=.true. at initialisation. The .csv particle files that will be found in the output folder will then be:
.csv files (one for each dtsave_particles, plus the initial state) containing 100 rows each (one row for each particle);.csv files (one for each of the three followed particles) containing 11 rows each (one row for each output time plus the initial state).An additional .csv file, containing the _destroyed label, may be present in the output folder. In this file, "destroyed" particles (i.e. particles removed from the domain) are stored as the simulation progresses.
Additional parameters in the particles_list of the .par file are available for refining the particle integration process:
relativistic: if .true., the relativistic equations of motion will be solved, instead of the Newtonian ones, when physics_type_particles='Lorentz' or physics_type_particles='GCA'.integrator_type_particles: if physics_type_particles='Lorentz' and relativistic=, the user can set this parameter to the preferred particle integrator, choosing among 'Boris' (for the Boris integrator), 'Vay' (for the Vay integrator), 'HC' (for the Higuera-Cary integrator), or 'LM' (for the Lapenta-Markidis integrator).eta_particles, etah_particles: when physics_type_particles='Lorentz' or physics_type_particles='GCA', the user can replace the resistivity or the hall resistivity from the MHD with a different (constant) resistivity which will be employed for calculating the electric field that enters the particle equations of motion.const_dt_particles: the user can define a constant time step for the particle integration, which will be used instead of calculating the time step from standard procedures.The particle module is still affected by some limitations. Features currently in the works (hence NOT working as of yet) include:
For more specific questions, email Fabio Bacchini, Bart Ripperda, or Jannis Teunissen.
1.9.1