 |
Sapphire++
v1.3.0-121-g2682132 |
|
Sapphire++
v1.3.0-121-g2682132 |
In this example we simulate the time-dependent transport and acceleration of particles at a parallel shock. The setup is similar to the quick start, but we use an advanced and more realistic configuration. You can find a detailed description of the setup in [15].
The astrophysical shock is assumed to move along the \(x\)-axis. In the rest frame of the shock wave, the background plasma's velocity is represented as \(\mathbf{U}_1 = U_{\rm sh} \hat{\mathbf{e}}_x\) upstream ( \(x<0\)) and \(\mathbf{U}_1 = U_{\rm sh}/r \hat{\mathbf{e}}_x\) downstream ( \(x>0\)). Here, \(r\) is the compression ratio of the shock and \(U_{\rm sh}\) is its velocity. As this represents a simple plain shock wave, we only require a single dimension in configuration space for the simulation. We also presuppose a constant magnetic field that is parallel to the shock normal, i.e. \(\mathbf{B} = B_0 \hat{\mathbf{e}}_x\). This kind of setup is referred to as a parallel shock.
In astrophysical shocks, particles can be accelerated to high energies through the process of diffusive shock acceleration [6] [12]. This involves injecting particles at the shock with a momentum of \(p_{\rm inj}\) at a rate \(Q\). These particles originate from the thermal background plasma and propagate around the shock. These particles may be scattered by Alvén waves travelling through the background plasma, a process we model with the help of scattering frequency \(\nu\). In [15] a constant scattering frequency is assumed, i.e. it is independent of \(p\) and \(x\):
\[ \nu = \nu_0 \, . \]
For the case of a constant scattering frequency, an analytic solution for the steady-state distribution function in mixed coordinates is presented in [6] :
\[ f(x, p, \theta) = \frac{Q}{p_{\rm inj}U_{1}} \frac{3r}{r - 1} \left(\frac{p}{p_{\rm inj}}\right)^{-3r/(r - 1)} \times \begin{cases} e^{3U_{1} \nu / V^{2} x} \left[1 - 3 U_{1}/V \cos\theta\right] &\quad \text{ for } x < 0 \\ 1 &\quad \text{ for } x \geq 0 \end{cases} \]
where \(V = p/\gamma m\) is the particle velocity in mixed coordinates. Decomposing this solution into spherical harmonics, we find that the isotropic part of the distribution function is given by
\[ f_{000}(x, p) = \sqrt{4\pi} \frac{Q}{p_{\rm inj}U_{1}} \frac{3r}{r - 1} \left(\frac{p}{p_{\rm inj}}\right)^{-3r/(r - 1)} \times \begin{cases} e^{3U_{1} \nu / V^{2} x} &\quad \text{ for } x < 0 \\ 1 &\quad \text{ for } x \geq 0 \end{cases} \]
Moreover, in the upstream we have a non-trivial anisotropy, namely
\[ f_{100}(x, p) = -3 U_{1}/V \sqrt{\frac{4\pi}{3}} \frac{Q}{p_{\rm inj}U_{1}} \frac{3r}{r - 1} \left(\frac{p}{p_{\rm inj}}\right)^{-3r/(r - 1)} e^{3U_{1} \nu / V^{2} x} \quad \text{ for } x < 0 \,, \]
while the downstream is isotropic, i.e.
\[ f_{100} = 0 \quad \text{ for } x \geq 0 \,. \]
An approximate analytic expression for the temporal evolution of the isotropic part of the distribution function \(f\) has been given in [7] and [8] . At the shock wave
\[ F(t, x = 0, p) = f_{000}(x = 0, p) \phi(t) \]
where \(F\) is the time-dependent isotropic part of the distribution function and
\[ \phi(t) = \frac{1}{2}\left[\exp\left(\frac{2 c^{2}_{1}}{c_{2}}\right) \mathrm{erfc}\left(\sqrt{\frac{c^{3}_{1}}{2 t c_{2}}} + \sqrt{\frac{c^{\phantom{3}}_{1} t}{2 c_{2}}}\right) + \mathrm{erfc}\left(\sqrt{\frac{c^{3}_{1}}{2 t c_{2}}} - \sqrt{\frac{c^{\phantom{3}}_{1}t}{2 c_{2}}}\right)\right] \,. \]
The constants \(c_1\) and \(c_2\) are the mean and the variance of the acceleration time. They are
\[ \begin{align} c_1 &= \frac{r}{2 U^{2}_{1} \nu} \frac{r + 1}{r - 1} \ln\left(\frac{1 + p^{2}}{1 + p^{2}_{\rm inj}}\right) \, \\ c_2 &= \frac{1}{3\nu^{2}} \frac{r}{U^{4}_{1}} \frac{r^{3} + 1}{r - 1} \left[ \frac{1}{1 + p^{2}} - \frac{1}{1 + p^{2}_{\rm inj}} + \ln\left(\frac{1 + p^{2}}{1 + p^{2}_{\rm inj}}\right)\right] \,. \end{align} \]
More details can be found in [15] .
A feature of Sapphire++, not covered in the quick start guide, is the use of parameter files. Its implementation uses the ParameterHandler class of deal.II. Parameter files allow you to define simulation parameters that are set at runtime, referred to hereafter as runtime parameters. This is in contrast to compile-time parameters.
Runtime parameters are read from a parameter file, such as parameter.prm, at the start of the simulation. You can pass a parameter file to Sapphire++ as the first command-line argument, e.g.
This enables you to modify simulation parameters without recompiling the code, which is particularly useful for parameter studies. Sapphire++ includes several predefined runtime parameters, such as the grid dimensions and the expansion order of the spherical harmonic expansion of \(f\).
Parameters in Sapphire++ are organised into categories and subcategories, primarily focusing on Output and VFP equation parameters. Key parameters are:
| Parameter Name | Category | Description |
|---|---|---|
| Results directory | Output | Specifies the directory for output files. |
| Format | Output | Defines the output format (e.g., vtu/pvtu/hdf5). |
| Output frequency | Output | Determines the frequency of output. |
| Simulation identifier | Output | Name of the simulation run, i.e. subfolder for the simulation. |
| Checkpoint frequency | Output | Time step frequency at which checkpoints are created. |
| Expansion order | VFP:Expansion | Sets the expansion order in spherical harmonics. |
| Polynomial degree | VFP:Finite Elements | Indicates the polynomial degree of the finite element method. |
| Grid type | VFP:Mesh | Specifies the grid type (e.g., Shock grid/Hypercube/File). |
| Final time | VFP:Time stepping | Sets the final time of the simulation. |
| Time step size | VFP:Time stepping | Determines the time step size of the simulation. |
| Method | VFP:Time stepping | Defines the time stepping method (e.g., CN/ERK4/...). |
A comprehensive and descriptive sample parameter file, template-parameter.prm, is located in the sapphirepp folder. This file is automatically generated when running Sapphire++ without any arguments. Here's an example of the syntax used in .prm files:
For more details on the predefined parameters in Sapphire++ we also refer to the documentation of the sapphirepp::Utils::OutputParameters and the sapphirepp::VFP::VFPParameters class.
Users can also define their own parameters to describe their physical setup. An example of this can be found in the Implementation section.
The implementation of this example can be found in the example directory sapphirepp/examples/vfp/parallel-shock. The following sections explain it in detail. As shown in the quick start, there are only a few lines in the config.h file that need to be adjusted.
As can be seen in the analytic solution presented above, the distribution function \(f\) depends only on \((x, p, \theta)\). This means that we need a single dimension in configuration space. Together with the \(p\) dependence this results in a two-dimensional reduced phase space. The \(\theta\) dependence is handled by the spherical harmonics with \(l>0\). Because we don't have any dependence on \(\varphi\), spherical harmonics with \(m>0\) will be zero.
Next, we need to specify which terms of the VFP equation are necessary to model the transport and acceleration of particles at a parallel shock. We actually have to solve the full VFP equation,
\[ \frac{\partial f}{\partial t} + (\mathbf{U} + \mathbf{v}) \cdot \nabla_{x} f - \gamma m \frac{\mathrm{D} \mathbf{U}}{\mathrm{D} t} \cdot \nabla_{p}f - \mathbf{p} \cdot\nabla_{x} \mathbf{U}\cdot \nabla_{p} f + q \mathbf{v} \cdot \left( \mathbf{B} \times \nabla_{p} f \right) = \frac{\nu}{2} \Delta_{\theta, \varphi} f + S \,. \]
We therefore need to activate all flags in the vfp_flags variable. In addition, we know that the velocity field and magnetic field are time independent. This also holds true for the source term.
This section outlines the process of implementing user-defined runtime parameters that describe their physical setup. The example parameters include:
| Name | Symbol | Code | Description |
|---|---|---|---|
| Shock velocity | \(U_{\rm sh}\) | u_sh | Velocity of the shock. |
| Compression ratio | \(r\) | compression_ratio | Compression ratio of the shock. |
| Injection momentum | \(p_{\rm inj}\) | p_inj | Momentum of the injected particles. |
| Injection rate | \(Q\) | Q | Rate of injected particles. |
| Scattering frequency | \(\nu_0\) | nu0 | Scattering frequency of the particles. |
| Magnetic field | \(B_0\) | B0 | Magnetic field strength. |
Additionally, there are several numerical parameters that are not part of the physical setup:
| Name | Symbol | Code | Description |
|---|---|---|---|
| Shock width | \(d_{\rm sh}\) | shock_width | Width of the shock. |
| Injection width in \(p\) | \(\sigma_p\) | sig_p | Width of the injection in momentum. |
| Injection width in \(x\) | \(\sigma_x\) | sig_x | Width of the injection in configuration space. |
| Injection position | \(x_{\rm inj}\) | x_inj | Position of the injection. |
Implementing these parameters in Sapphire++ involves a two-step process:
Define the parameters within the PhysicalParameters class in the config.h file. This step informs the compiler about their existence and sets default values to be used if the parameter is not specified in the parameter file.
Add the parameters to the parameter file to ensure that the parameter parser expects them and automatically sets their values. The deal.II class ParameterHandler prm and its add_parameter method are used for this purpose:
The declare_parameters() method is edited to include the following code:
We now implement the functions modelling the physical scenario. We start with the scattering frequency. As stated above, we assume a constant scattering frequency, i.e. \(\nu = \nu_0\). This translates to
The function calculates the scattering frequency at multiple points in the \((x,\ln p)\)–domain. We use the index q_index to refer to an individual point. Notice that we access the runtime parameter \(\nu_0\) via the (previously modified) PhysicalParameters class prm.
Furthermore, we would like to inject mono-energetic particles at the shock. To model this we use a Gaussian profile centred at \((x_{\rm inj}, p_{\rm inj})\) with width \(\sigma_x\) and \(\sigma_p\),
\[ S(x, p) = \frac{Q}{2\pi\sigma_x\sigma_p} \exp\left[ -\frac{(x - x_{\rm inj})^2}{2\sigma_x^2} - \frac{(p - p_{\rm inj})^2}{2\sigma_p^2} \right] \, . \]
As already explained in the quick start, the source term needs to be decomposed into spherical harmonics. As we inject an isotropic distribution of particles, we can set all \(l>0\) components to zero. Notice that we use the index i to refer to the components of the spherical harmonic decomposition, \(i(l,m,s)\), at a single point in the \((x,\ln p)\)–domain.
We emphasise that the second component of point, i.e. point[1], corresponds to the momentum variable \(\ln p\).
As described above, we use a constant magnetic field parallel to the shock, \(\mathbf{B} = B_0 \hat{\mathbf{e}}_x\). It does not influence the steady-state solution, but nevertheless we include it for completeness.
The velocity field is the same as in the quick start example. We use a \(\tanh\) profile to model the transition at the shock wave. The only modification is, that we now use the user-defined runtime parameters \(U_{\rm sh}\), \(r\) and \(d_{\rm sh}\). Otherwise, \(\mathbf{U}(\mathbf{x})\), its divergence \(\nabla \cdot \mathbf{U}(\mathbf{x})\), the material derivative \(\mathrm{D} \mathbf{U}/\mathrm{D} t\) and the Jacobian \(\partial U_{x}/\partial x\) stay the same.
Background velocity field value \(\mathbf{U}(\mathbf{x})\):
Background velocity divergence \(\nabla \cdot \mathbf{U}(\mathbf{x})\):
Background velocity material derivative \(\mathrm{D}\mathbf{U}/\mathrm{D} t\):
Background velocity Jacobian \(\partial U_{x}/ \partial x\):
To run the simulation, implement the above functions in the config.h file and recompile Sapphire++:
Alternatively, you can use the implementation in the examples/vfp/parallel-shock folder. This executable will be named parallel-shock instead of sapphirepp. Note that for this to work, Sapphire++ must be configured with the -DEXAMPLES=ON option.
We recommend running the simulation with the parameters given in examples/vfp/parallel-shock/parameter.prm:
Run the simulation with:
We provide a sapphireppplot script to generate the plots:
Sometimes simulations fail. This can be related to a bug - in that case you probably need to rerun the simulation from start after fixing it. But it can also be due to hardware failure or because of time limits for large simulations on HPC clusters. In this case it is imperative to be able to restart simulation from an earlier time step. For this special output files are needed, called checkpoints. (The normal output is processed so it can be easily visualized, but it looses information on the underlying discontinuous Galerkin (dG) degrees of freedom (DoFs).)
To enable checkpointing, you need to set the Checkpoint frequency in the parameter file. In this example we use a Checkpoint frequency = 50, which means that a checkpoint will be created every 50 simulation time steps. This is quite excessive and chosen for educational purposes. A good checkpoint frequency is a compromise of how long it takes to run a time step, the expected failure rate/cluster time limit and the time it needs to write the checkpoints. See the deal.II checkpoint/restart tutorial for more details.
After running the simulation, you can find these checkpoints in the results folder:
yields
These checkpoint files save basic metadata, the grid and solution DoFs. Only the most recent checkpoint is kept to limit storage usage.
Restarting a simulation is done with the --resume (or -r) command-line flag. Since the checkpoint does not save the parameters used for the simulation, you also have to provide a parameter file to Sapphire++. This allows you to change some parameters between restarts. But to keep parameters consistent, the log.prm file can be used:
Note, that the number of MPI processes can differ for the restart. In case of a successful restart, you should see the following message:
Before we analyse the results, we check that the solution has converged to a steady state. To demonstrate this, we show the time series of the 2D data in a logarithmic scale.
From the first figure we can conclude that the solution reached steady state even for large momenta. The second figure shows a cut-out of the solution in the shock region. This serves an illustrative purpose demonstrating the stretched Shock grid used in the simulation.
The analytic solution predicts that the energy spectrum of the particles at shock follows a \(p^{-3r/(r - 1)}\) power law for \(p > p_{\rm inj}\). Given a compression ratio \(r=4\), we anticipate a \(p^{-4}\) power law. We plot \(f_{000}(x = 0, \ln p)\) and compare it with our expectation.
In the upstream region ( \(x<0\)), we expect an exponential cut-off, which is dependent on the scattering frequency \(\nu\). As discussed in the introduction, we also observe an upstream anisotropy, \(f_{100}\). Note that the anisotropic part disappears in the downstream region. When comparing our simulation results to the analytic solution, we find they are in good agreement. The figure below shows \(f_{000}(x, p = 10)\) and \(f_{100}(x, p = 10)\):
In the last figure we present the temporal evolution of the energy spectrum at the shock for a specific \(p\), namely \(p=10\). The simulation results are compared to the approximate analytical solution mentioned in the introduction.
We note that there exists an exact analytical solution for a specific choice of the scattering frequency in the up- and downstream, see, for example, eq. 16 and 21 in [7] . In [4] it is demonstrated that Sapphire++ reproduces it accurately.
1.16.1