Hybrid Code Crash

[Beforerr]

My init setup is an large-amplitude oblique linearly polarized Alfvén wave using Hybrid PIC model. The wave is propagating in the z-direction.

The code could run for some time, but then crash. The error message is attached below, and the input file is also attached. I am not sure how to debug this issue. Also I observe that magnetic field initiation is not expected in 1D/2D code as there are no variation among z axis. (Simulation also crashed for 1D/2D)

 0: amrex::BLBackTrace::print_backtrace_info(__sFILE*) (in libamrex_3d.dylib) + 64

 1: amrex::BLBackTrace::handler(int) (in libamrex_3d.dylib) + 688

 2: _sigtramp (in libsystem_platform.dylib) + 56

 3: 0x00000001468f37c0 (in warpx_pybind_3d.cpython-311-darwin.so) + 232

 4: 0x00000001468f37c0 (in warpx_pybind_3d.cpython-311-darwin.so) + 232

 5: __kmp_invoke_microtask (in libomp.dylib) + 156

===== TinyProfilers ======
REG::WarpX::Evolve()
WarpX::Evolve()
WarpX::Evolve::step
WarpX::HybridPICEvolveFields()
ablastr::utils::communication::FillBoundary
FabArray::FillBoundaryAndSync()
FillBoundaryAndSync_nowait()
FillBoundary_nowait()

# algo
algo.current_deposition = direct
algo.maxwell_solver = hybrid
algo.particle_shape = 1

# amr
amr.max_grid_size = 32
amr.n_cell = 8 8 128

# boundary
boundary.field_hi = periodic periodic periodic
boundary.field_lo = periodic periodic periodic
boundary.particle_hi = periodic periodic periodic
boundary.particle_lo = periodic periodic periodic

# diag1
diag1.diag_type = Full
diag1.intervals = 8
diag1.write_species = 0

# diag2
diag2.diag_type = Full
diag2.fields_to_plot = none
diag2.intervals = 8
diag2.species = ions
diag2.write_species = 1

# diagnostics
diagnostics.diags_names = diag1 diag2

# geometry
geometry.dims = 3
geometry.prob_hi = 85025.49218531257 85025.49218531257 2720815.7499300023
geometry.prob_lo = -85025.49218531257 -85025.49218531257 0

# hybrid_pic_model
hybrid_pic_model.elec_temp = 61.930612282580206
hybrid_pic_model.n0_ref = 100000000.0
hybrid_pic_model.n_floor = 1562500.0
hybrid_pic_model.plasma_hyper_resistivity = 1e-06
hybrid_pic_model.plasma_resistivity(rho,J) = 1e-06
hybrid_pic_model.substeps = 32

# ions
ions.charge = q_e
ions.density_function(x,y,z) = 100000000.0
ions.initialize_self_fields = 0
ions.injection_style = nrandompercell
ions.mass = 3.6437534806e-28
ions.momentum_distribution_type = gaussian_parse_momentum_function
ions.momentum_function_ux_m(x,y,z) = (0)/299792458.0
ions.momentum_function_ux_th(x,y,z) = (165018.7869319539)/299792458.0
ions.momentum_function_uy_m(x,y,z) = (116685.90326276264 * cos(9.237208079733046e-06 * z))/299792458.0
ions.momentum_function_uy_th(x,y,z) = (165018.7869319539)/299792458.0
ions.momentum_function_uz_m(x,y,z) = (0)/299792458.0
ions.momentum_function_uz_th(x,y,z) = (165018.7869319539)/299792458.0
ions.num_particles_per_cell = 256
ions.profile = parse_density_function

# max_step
max_step = 800

# particles
particles.species_names = ions

# warpx
warpx.B_ext_grid_init_style = parse_b_ext_grid_function
warpx.Bx_external_grid_function(x,y,z) = 8.660254037844386e-08
warpx.By_external_grid_function(x,y,z) = 5.0000000000000004e-08 * cos(9.237208079733046e-06 * z)
warpx.Bz_external_grid_function(x,y,z) = 5.000000000000002e-08
warpx.const_dt = 0.017861933764391047

[roelof-groenewald]

Hi @Beforerr. Thanks for your question and welcome to the WarpX community!

I'd say the first step in debugging this would be to follow-up on the clue of your initial B-field not being what you expect. You say you don't see variation in z for 1d and 2d. But do you get the expected field in 3d?

Could you also please provide some information about how you are running this? Is it on CPU or GPU? How many cores and OMP threads?

[Beforerr]

Hi @roelof-groenewald , thanks for prompt reply.

For 3d, by visual inspection it is fine and the plasma evolved as expected for like 50 ion gyro period then the simulation crashed.

For testing purpose , it is run on CPU with 1 MPI and 12 OMP threads. (If you want to take a look at 1D and 2D cases. The input files could be found at this repo)

[roelof-groenewald]

@Beforerr I'm not able to see the repo you linked - I get 404 Page not found error, likely due to the repo being private.

However, I was able to reproduce the error you saw from your 3d input.

Generally speaking the hybrid model approach has difficulty with Whistler waves that are included in the model physics but because their dispersion runs to the electron cyclotron frequency there is a point where the chosen grid and timestep resolution no longer properly resolves these waves (unless you use a very small timestep in which case the value of the hybrid approach is lost). At this cutoff, the unresolved Whistler waves can cause numerical heating. This effect is counteracted by the inclusion of a finite resistivity to damp these waves. All this is to say, at a given grid resolution the hybrid approach will be unstable if the timestep is too large or the resistivity is too low. I believe this is what you've encountered. Here is a plot of the energy evolution in the system with your input parameters: (note the energy transfer back-and-forth between the ions and magnetic field before the simulation goes unstable).

If I increase the resistivity to 500 I get the following evolution:

The resistivity in this model should not be thought of as a physical resistivity (such as Spitzer resistivity) but rather a numerical resistivity used mainly to damp high frequency Whistler waves. You should play around with the chosen timestep and resistivity value to find a suitable point where the timestep is large enough for your simulation to progress quickly while the resistivity is high enough to damp the unresolved whistler waves without damping desired modes too strongly. You can also increase the number of substeps used in the solver to increase the effective time resolution for the solver (increase the frequency resolution peak). This helps because wave damping at a given resistivity occurs faster for higher frequency waves so you can again use a lower resistivity the more substeps you use. You might also want to increase the spatial resolution as I notice that the ion interial length is just barely resolved.

I hope this helps to clear up your problem. Like I said, I couldn't access your 1d or 2d input scripts to try and reproduce the issue with the initial magnetic field not loading properly. Could you maybe post the input script for one of those cases with the plotting script you used to visualize the fields?

[Beforerr]

Hi @roelof-groenewald , I have made the repo public, sorry for the confusion. You might check the input scripts first and soon I would update the plotting script.

Previously, I also doubted this was caused by the whistler wave. But the CFL I calculated including substeps is already 0.02, much smaller than 1. Do you think that numerical dispersion would still cause the problem?

Also I am just wondering which resistivity you mentioned, is it simple or hyper-resistivity? Previously, I played with the timestep and the substeps, and the problem still exists (just a matter of time). But you may be right, I will try to increase the resistivity and see if it helps.

[roelof-groenewald]

Sorry for the slow response this time around.

I haven't had a chance yet to look at your 1d and 2d results, but to answer your other questions:

Do you think that numerical dispersion would still cause the problem?

Are you referring to CFL = 0.02 using the speed of light or the ion thermal velocity? I assume the latter. But the fact that the simulation stability definitely improved with higher resistivity does suggest to me that the problem is with numerical dispersion.

Also I am just wondering which resistivity you mentioned, is it simple or hyper-resistivity?

I was changing the simple resistivity. The hyper-resistivity is effective at suppressing noise from low density regions (where the density is at or below the floor density value). You can also play with that a bit to see how it affects the algorithm stability though.

[Beforerr]

Hi @roelof-groenewald , thanks for your explanation. I use the following formula to calculate CFL (where $di= c/\omega{pi}$)

$$\Omega_{ci} \Delta t < (\Delta x / d_i)^2 / \pi$$

And here I check the magnetic field and plasma velocity in y direction.

For 1D (here the x label should be z due to yt loading):

For 3D:

You can see the By is not properly initialized in 1D (as expected), but it is in 3D.

Attached is the minimal script to plot the magnetic field data.

# %%
import os
import yt
from functools import partial

dim = 1
directory = "dim_1_beta_0.25_theta_60_eta_200"
field_diag_dir = "diags/diag1"
file_format = "??????"

# %%
os.chdir(directory)
ts_field = yt.load(f"{field_diag_dir}{file_format}")
ds_field = ts_field[0]

# because `yt` will always load the data in 3D, we need to specify the direction corresponding to the simulation direction `z`
if dim == 3:
    direction = "z"
elif dim == 2:
    direction = "y"
else:
    direction = "x"

# %%
plot_field = partial(
    yt.ProfilePlot,
    x_field=direction,
    weight_field=("boxlib", "volume"),
    x_log=False,
    y_log=False,
)

plot_field(ds_field, y_fields="By")

[roelof-groenewald]

I see what you mean by the By field not being initialized properly in 1d. Hopefully I will have some time tomorrow to look into why that is happening.

But in the meantime I wanted to ask, why do you square $\Delta x/di$ in the CFL-like condition calculation? I usually check that $\Omega{ci}d_i \Delta t = v_A \Delta t < \Delta x$ as a CFL-like condition (where $v_A$ is the Alfven speed) and keep $\Delta x < di$. In cases where super-Alfvenic ions exist it is also good to ensure that $v{max}\Delta t < \Delta x$ as a "particle CFL-like" condition.

[RemiLehe]

@Beforerr Thanks for reporting the issue about the external fields in 1D and 2D.

I can reproduce the issue in 1D, but not in 2D. (i.e., in 2D, when using the plotting script that you posted above, I do see oscillations in By).

Here are the simulation scripts that I used inputs_2d.txt inputs_1d.txt

[RemiLehe]

This PR should fix the issue: https://github.com/ECP-WarpX/WarpX/pull/4906

[Beforerr]

I see what you mean by the By field not being initialized properly in 1d. Hopefully I will have some time tomorrow to look into why that is happening.

But in the meantime I wanted to ask, why do you square Δx/di in the CFL-like condition calculation? I usually check that ΩcidiΔt=vAΔt<Δx as a CFL-like condition (where vA is the Alfven speed) and keep Δx<di. In cases where super-Alfvenic ions exist it is also good to ensure that vmaxΔt<Δx as a "particle CFL-like" condition.

Hi @roelof-groenewald , I think you are propably right about CFL condition. It is worth to check the `particle CFL-like` condition. My understanding is that because whistler have dispersion relationship of w = k² , so that spatial step should be squared (I refer this paper https://doi.org/10.1063/5.0146529 ).

[Beforerr]

@RemiLehe . Great thanks for the quick fix. I check and you are right that 2D has no initialization problem (sorry for the confusion).

Just want to make sure that for 1D/2D output, x and Bx are not referring to the same direction and x is actually z?

[RemiLehe]

@Beforerr When using the plotfile format and visualizing the output with yt: yes, I think that you are right x and Bx do not refer to the same direction (which I agree is confusing).

When using the openPMD output (and using e.g. openPMD-viewer): quantities are labeled more consistently. In particular, the z spatial coordinate will be correctly labeled as z.