search

SmileiPIC / Smilei

Particle-in-cell code for plasma simulation
https://smileipic.github.io/Smilei
335 stars 119 forks source link

running faults related to momentum_initialization using numpy array in namelist.py #628

Closed zk-fdu closed 1 year ago

zk-fdu commented 1 year ago

Description

Hi,developers, I encountered a runniing error related to momentum_initialization using numpy array in namelist.py.When i setup momentum_initialization as the product of a number ranging from 0 to 1 and a number representing species proton with mass 1836, the result would not be output normally. The following is momentum_initialization of species proton using array in namelist.py . ``

momentum_initializtion defined as array_momentum

array_momentum = np.zeros((3,npart)) # momenta x,y,z mass_proton = 1836.15267343 ## proton mass array_momentum[0,:] = np.multiply(0.0,np.ones(npart))mass_proton array_momentum[1,:] = np.multiply(0.5,np.ones(npart))mass_proton array_momentum[2,:] = np.multiply(0.1,np.ones(npart))*mass_proton``

i.e ,assuming simulation_time is defined as 5, the running time maybe 1.2 and abruptly.

I don't know where error is. When i use happi block via Python3, the hint and picture plot of px vs. time is given .


happi/_Diagnostics/TrackParticles.py:631: RuntimeWarning: invalid value encountered in greater
  self._rawData['brokenLine'] += self._np.abs(dudt).max(axis=0) > 1.
... done

## Parameters

- Smilei version: 4.7
- Input file : as enclosure
-
[namelist.py.txt.txt](https://github.com/SmileiPIC/Smilei/files/11595415/namelist.py.txt.txt)

![px](https://github.com/SmileiPIC/Smilei/assets/83404979/0a4221ba-aa27-43ce-9866-836acd468343)

- Computing resources
  - C++ compiler (to get some information `g++ --version`): [GCC 9.4.0]
  - Python version (`python --version`): Python 3.8.10
zk-fdu commented 1 year ago

When i use happi block via Python3, the hint and picture plot of px vs. time is given .

hint

happi/_Diagnostics/TrackParticles.py:631: RuntimeWarning: invalid value encountered in greater self._rawData['brokenLine'] += self._np.abs(dudt).max(axis=0) > 1. ... done

namelist.py.txt.txt px

zk-fdu commented 1 year ago

I setup momentum_initialziation in units of according to Unit document. I checked the hint /happi/_Diagnostics/TrackParticles.py:631, which says the reason is caused by velocity> c；however, i was not against the rule. i don't know where the error is. The pointed code will be attached as enclosure. array_momentum = np.zeros((3,npart)) # momenta x,y,z mass_proton = 1836.15267343 ## proton mass array_momentum[0,:] = np.multiply(0.0,np.ones(npart))mass_proton array_momentum[1,:] = np.multiply(0.5,np.ones(npart))mass_proton array_momentum[2,:] = np.multiply(0.1,np.ones(npart))*mass_proton TrackParticles-631

mccoys commented 1 year ago

I don't understand where is the bug. What happens exactly? Does your simulation run? Is it a problem with post processing only?

zk-fdu commented 1 year ago

Hi,Mccoys！ I encountered a running error related to momentum_initializtion.Let me introduce some specific details.Firstly,my simulation can run normally during simulation time;however, the output is wrong:the output is as usual when code run in former part of simulation while the latter is wrong,during which output is NaN. Then,i show my setup about momentum_initializtion.In the namelist.py.txt, I setup momentum_initialization as the product of a number ranging from 0 to 1 and a number representing species proton with mass 1836,in units of $m_e c$ according to Unit document. The following is momentum_initialization of species proton using array in namelist.py . ``

momentum_initializtion defined as array_momentum

array_momentum = np.zeros((3,npart)) # momenta x,y,z mass_proton = 1836.15267343 ## proton mass array_momentum[0,:] = np.multiply(0.0,np.ones(npart))mass_proton array_momentum[1,:] = np.multiply(0.5,np.ones(npart))mass_proton array_momentum[2,:] = np.multiply(0.1,np.ones(npart))*mass_proton``

Finally,when running code was over and h5 file is dealt with, some strange things appeared.Let's take the momentum in the x direction for example（px.pngpx vs.time).During post processing via happi block of Python , the hint given is "Smilei-4.7/happi/_Diagnostics/TrackParticles.py:631: RuntimeWarning: invalid value encountered in greater self._rawData['brokenLine'] += self._np.abs(dudt).max(axis=0) > 1." For the picture of px vs time,within about 1.2 T_r, the data is normal,while the rest is blank abruptly.[assuming simulation_time is defined as 5].When I checked the hint /happi/_Diagnostics/TrackParticles.py:631, which says the reason is caused by velocity> c；however, i was not against this rule. i don't know where the error is. The pointed code will be attached as enclosure(TrackParticle-631.png).I would appreciate it very much if you should help me! Best wishes, Zhao

Parameters Smilei 4.7 Input file: namelist.py.txt Computing resources Linux: Ubuntu 20.04.5 LTS (GNU/Linux 5.14.0-1059-oem x86_64) C++ compiler g++ [GCC 9.4.0] on linux Python version: Python 3.8.10

-----原始邮件----- 发件人:mccoys @.> 发送时间:2023-05-30 14:43:13 (星期二) 收件人: SmileiPIC/Smilei @.> 抄送: zk-fdu @.>, "State change" @.> 主题: Re: [SmileiPIC/Smilei] running faults related to momentum_initialization using numpy array in namelist.py (Issue #628)

I don't understand where is the bug. What happens exactly? Does your simulation run? Is it a problem with post processing only?

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you modified the open/close state.Message ID: @.***>

----------------------------------------------------------------------------------------

SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI

----------------------------------------------------------------------------------------

Physical Inputs in normalized units, converted from SI units

import math import numpy as np

import scipy

from scipy.constants import c, epsilon_0, e, m_e # constants in SI units wavelength_SI= 0.8e-6 # laser wavelength, m

Normalization quantities[same as Smilei Units]

L_r = wavelength_SI # reference length is the wavelength wr = c/L_r # reference angular frequency T_r =1/wr #reference time

Laser wavelength and period in normalized units

l0 = wavelength_SI / L_r t0 = wavelength_SI/c / T_r

Grid size and total simulated time in normalized units

Lx =6. wavelength_SI / L_r Ly =4.wavelength_SI / L_r Lz = 4.wavelength_SI / L_r Tsim =5wavelength_SI/c / T_r # duration of the simulation

#################### Simulated domain and time interval #######################

Spatial and temporal resolution

resx = 50 # nb of cells in one laser wavelength

rest = 15. # nb of timesteps in one optical cycle

Mesh and integration timestep

dx = l0/resx dy = l0/resx dz = l0/resx dt = 0.95/math.sqrt(1./(dxdx) + 1./(dydy)+1./(dz*dz)) cell_length = [dx,dy,dz] grid_length = [Lx,Ly,Lz]

Namelists

Main( geometry = "3Dcartesian",

interpolation_order = 2 ,

cell_length = cell_length,
grid_length = grid_length,

number_of_patches = [4,4,4],

#time_fields_frozen = simulation_time,
maxwell_solver="Yee",
timestep = dt,
simulation_time = Tsim,

EM_boundary_conditions = [
    ["silver-muller", "silver-muller"],
    ["silver-muller", "silver-muller"],
    ["silver-muller", "silver-muller"],
],
solve_relativistic_poisson="True",
#solve_poisson="False",   
print_every = 1,

reference_angular_frequency_SI = wr,
#patch_arrangement = "linearized_XYZ"

)

T=0.1

npart=1 lambda0 = 0.8e-6 # laser wavelength, m c=299792458 # lightspeed, m/s omega0 = 2math.pic/lambda0 # laser frequency eps0=8.854187e-12 # Vacuum permittivity, F/m e=1.60217646e-19 # Elementary charge, C me=9.10938215e-31 # Electron mass, kg ncrit = eps0*omega02*me/e*2 # Plasma critical number density [m-3] normalized_species_charge = -1 # For electrons Q_bunch = -1e-19 Q_part = Q_bunch/npart cell_volume = 0.2(c/omega0)3 weight = Q_part/(cell_volumencrite*normalized_species_charge)

array_position = np.zeros((4,npart)) # positions x,y,z, weight array_momentum = np.zeros((3,npart)) # momenta x,y,z array_position[0,:] = np.multiply(0.01,np.ones(npart)) array_position[1,:] = np.multiply(0.01,np.ones(npart)) array_position[2,:] = np.multiply(0.01,np.ones(npart)) array_position[3,:] = np.multiply(np.ones(npart),weight) mass_proton = 1836.15267343 ## proton mass array_momentum[0,:] = np.multiply(0.0,np.ones(npart))mass_proton array_momentum[1,:] = np.multiply(0.5,np.ones(npart))mass_proton array_momentum[2,:] = np.multiply(0.1,np.ones(npart))*mass_proton

__

LoadBalancing( every = 20, cell_load = 1., frozen_particle_load = 0.1 ) Species( name = "proton", position_initialization =array_position, momentum_initialization =array_momentum, mass = 1836.15267343, charge = 1.0, pusher = "boris", boundary_conditions = [ ["remove", "remove"], ["remove", "remove"], ["remove", "remove"], ],

)

constant B field along x-y plane

B_field_amplitude =5000 B_field_direction = np.array([0,0,1]) B_field_direction = B_field_direction / np.linalg.norm(B_field_direction) B_field_vector = B_field_amplitude* B_field_direction

ExternalField( field = "Bx", profile = constant(B_field_vector[0]) )

ExternalField( field = "By", profile = constant(B_field_vector[1]) )

ExternalField( field = "Bz", profile = constant(B_field_vector[2]) )

Vectorization(

mode = "off",

)

DiagFields( every = 10, fields = ['Ex','Ey','Ez','Bx','By','Bz'] )

DiagTrackParticles( species = "proton", attributes = ["x","y","z","px","py", "pz","w"], every = 1, flush_every = 1

)

mccoys commented 1 year ago

I believe the particle has simply exited the box so it cannot appear anymore. It has nothing to do with the Warning

zk-fdu commented 1 year ago

What you mean is that if I add simulation grid length，this phenomenon will be improved?

mccoys commented 1 year ago

No. The particle goes back to negative values of x. There is no box for negative x. So you should place your particle further inside the box (larger x) at the beginning

zk-fdu commented 1 year ago

Oh,I see. However, I setup species proton with positive charge +1 and positive components of velocity, and does the particle move along the positive axis ?how the particle goes back to negative values of x? This point confuses me terribly,!

zk-fdu commented 1 year ago

No. The particle goes back to negative values of x. There is no box for negative x. So you should place your particle further inside the box (larger x) at the beginning

what you mean is there is no distinction between positive axis and negative axis for the box ,and the distinction between positive and negative direction is based on species initialzition? If I make the particle possess the velocity +0.1c along x axis , the increasing index of grid is positive direction ?

mccoys commented 1 year ago

No

The box starts at x = 0, up to the value you indicate in grid_length. The particles cannot go outside the box. If they cross the boundary, they are either reflected or destroyed depending on what you chose in boundary_conditions

If you set a particle velocity equal to 0.1 c, it will move towards positive x. But in your graph, it acquires a negative value of velocity so that it goes back towards small values of x until probably exiting the box

zk-fdu commented 1 year ago

Oh， I see what you mean. I have a try. Thanks for your help!