-
Notifications
You must be signed in to change notification settings - Fork 191
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Hybrid Code Crash #4883
Comments
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? |
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) |
@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: 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 I hope this helps to clear up your problem. |
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. |
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:
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.
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. |
Hi @roelof-groenewald , thanks for your explanation. I use the following formula to calculate And here I check the magnetic field and plasma velocity in For 1D (here the x label should be You can see the 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") |
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 |
@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 Here are the simulation scripts that I used |
This PR should fix the issue: |
Hi @roelof-groenewald , I think you are propably right about CFL condition. It is worth to check the `particle CFL-like` condition. |
@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, |
@Beforerr When using the When using the |
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)
The text was updated successfully, but these errors were encountered: