Warwick-Plasma / epoch

Particle-in-cell code for plasma physics simulations
https://epochpic.github.io
GNU General Public License v3.0
170 stars 54 forks source link

Numerical instability in two-dimensional simulations #674

Open naokikatsura opened 1 week ago

naokikatsura commented 1 week ago

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.

ne0020 ni0020

TomGoffrey commented 1 week ago

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.

naokikatsura commented 1 week ago

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