Numerical instability in two-dimensional simulations
Closed this issue · 4 comments
Hello everyone,
Two-dimensional plasma thermal expansion results show a speckled density distribution. How can I solve this probelm?
I am currently using EPOCH2D version 4.19.2, and the simulation was started with a plasma of homogeneous density and temperature set at a radius of 0.1 mm from the origin. There is no external electromagnetic field, and the plasma is assumed to simply expand thermally.
Here is the spacial distribution of electron and proton.
Can't say I've seen anything like that before....can you post a full input deck? Best I can do right now is check if I can see anything obviously wrong.
Thank you for your reply. Here is my input.deck.
--- CONTROL ---
begin:control
nx = 4000
ny = 4000
t_end = 1e-9
x_min = -2000. * 5e-7
x_max = 2000. * 5e-7
y_min = -2000. * 5e-7
y_max = 2000. * 5e-7
#stdout_frequency = 10
balance_first = T
use_optimal_layout = T
smooth_currents = T
npart = 2e7
physics_table_location = /epoch/epoch2d/src/physics_packages/TABLES
end:control
--- BOUNDARY CONDITION ---
begin:boundaries
bc_x_min = open
bc_x_max = open
bc_y_min = open
bc_y_max = open
end:boundaries
--- CONSTANT ---
begin:constant
rc = 6e-4
xc = -rc
yc = rc
r0 = yc/6.0
T0 = 1000 # initial plasma temprature [eV]
n0 = 1e+24 # initial plasma density [m^-3]
end:constant
--- PROTON ---
begin:species
name = proton
charge = 1.0
mass = 100.0
frac = 0.5
bc_x_min = open
bc_x_max = open
bc_y_min = open
bc_y_max = open
number_density = if(sqrt(xx+yy) lt r0, n0, 0.0)
temp_ev = if(sqrt(xx+yy) lt r0, T0, 0.0)
end:species
--- ELECTRON ---
begin:species
name = electron
charge = -1.0
mass = 1.0
frac = 0.5
bc_x_min = open
bc_x_max = open
bc_y_min = open
bc_y_max = open
number_density = if(sqrt(xx+yy) lt r0, n0, 0.0)
temp_ev = if(sqrt(xx+yy) lt r0, T0, 0.0)
end:species
--- OUTPUT ---
begin:output
name = normal
dt_snapshot = 1e-11
nstep_average = 10
Properties at particle positions
particles = always
#px = always
#py = always
#pz = always
vx = always
vy = always
vz = always
id = always
ex = always + average
ey = always + average
ez = always + average
bx = always + average
by = always + average
bz = always + average
jx = always + average
jy = always + average
jz = always + average
ppc = always + species
particle_weight = always
Properties on grid
ekbar = always + species
grid = always
number_density = always + species + average
temperature = always + species + average
total_energy_sum = always
end:output
begin:output
name = Restart
file_prefix = restart
restartable = T
dt_snapshot = 1e-11
end:output
It seems that this instability comes from the grid size. In this case, the grid size is set to be twice as large as the initial debye length. However, when I decrease the number density of ions and electrons, the result is reasonable.
Happy to hear you resolved the issue. The typical rule for grid resolution is to have cell-sizes which are close to the Debye length. Sometimes you can get away with larger cells if you use higher order particle shapes, or current-smoothing. The discussion on self heating covers some of the ways simulation parameters may be relaxed.
I'll close the issue now, but feel free to continue commenting if you run into problems later.
Cheers,
Stuart