Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
Check that the relation div(E) - rho/epsilon_0 is satisfied, with a
relative precision close to the machine precision (directly in spectral space)
Parameters
----------
sim: Simulation object
rho_ions: list of 2d complex arrays (one per mode)
The density of the ions (which are not explicitly present in the `sim`
object, since they are motionless)
"""
# Create a global field object across all subdomains, and copy the fields
global_Nz, _ = sim.comm.get_Nz_and_iz(
local=False, with_damp=False, with_guard=False )
global_zmin, global_zmax = sim.comm.get_zmin_zmax(
local=False, with_damp=False, with_guard=False )
global_fld = Fields( global_Nz, global_zmax,
sim.fld.Nr, sim.fld.rmax, sim.fld.Nm, sim.fld.dt,
zmin=global_zmin, n_order=sim.fld.n_order, use_cuda=False)
# Gather the fields of the interpolation grid
for m in range(sim.fld.Nm):
# Gather E
for field in ['Er', 'Et', 'Ez' ]:
local_array = getattr( sim.fld.interp[m], field )
gathered_array = sim.comm.gather_grid_array( local_array )
setattr( global_fld.interp[m], field, gathered_array )
# Gather rho
global_fld.interp[m].rho = \
sim.comm.gather_grid_array( sim.fld.interp[m].rho + rho_ions[m] )
# Loop over modes and check charge conservation in spectral space
if sim.comm.rank == 0:
global_fld.interp2spect('E')
Returns
-------
A dictionary containing :
- 'E' : 1d array containing the values of the amplitude
- 'w' : 1d array containing the values of waist
- 'fld' : the Fields object at the end of the simulation.
- 'z_centroid' : 1d array containing the values of the centroid position
"""
# Choose the timestep, so that the simulation advances by
# one cell at every timestep
dt = Lz*1./Nz * 1./c
# Initialize the fields object
fld = Fields( Nz, Lz, Nr, Lr, Nm, dt, v_galilean = 0.999*c )
z0 = Lz/2
init_fields( fld, w0, ctau, k0, z0, E0, m )
# Create the arrays to get the waist and amplitude
w = np.zeros(Nt)
E = np.zeros(Nt)
z_center = np.zeros(Nt)
# Get the fields in spectral space
fld.interp2spect('E')
fld.interp2spect('B')
#Create moving window object
mov_win = MovingWindow( fld.interp[0], v= 0. )
# Loop over the iterations
# Project the charge and currents onto the local subdomain
sim.deposit( 'rho', exchange=True, species_list=[ptcl],
update_spectral=False )
sim.deposit( 'J', exchange=True, species_list=[ptcl],
update_spectral=False )
# Create a global field object across all subdomains, and copy the sources
# (Space-charge calculation is a global operation)
# Note: in the single-proc case, this is also useful in order not to
# erase the pre-existing E and B field in sim.fld
global_Nz, _ = sim.comm.get_Nz_and_iz(
local=False, with_damp=True, with_guard=False )
global_zmin, global_zmax = sim.comm.get_zmin_zmax(
local=False, with_damp=True, with_guard=False )
global_fld = Fields( global_Nz, global_zmax,
sim.fld.Nr, sim.fld.rmax, sim.fld.Nm, sim.fld.dt,
n_order=sim.fld.n_order, smoother=sim.fld.smoother,
zmin=global_zmin, use_cuda=False)
# Gather the sources on the interpolation grid of global_fld
for m in range(sim.fld.Nm):
for field in ['Jr', 'Jt', 'Jz', 'rho']:
local_array = getattr( sim.fld.interp[m], field )
gathered_array = sim.comm.gather_grid_array(
local_array, with_damp=True )
setattr( global_fld.interp[m], field, gathered_array )
# Calculate the space-charge fields on the global grid
# (For a multi-proc simulation: only performed by the first proc)
if sim.comm.rank == 0:
# - Convert the sources to spectral space
global_fld.interp2spect('rho_prev')
self.boost = None
# Register time step
self.dt = dt
# Initialize the boundary communicator
cdt_over_dr = c*dt / (rmax/Nr)
self.comm = BoundaryCommunicator( Nz, zmin, zmax, Nr, rmax, Nm, dt,
self.v_comoving, self.use_galilean, boundaries, n_order,
n_guard, n_damp, cdt_over_dr, None, exchange_period,
use_all_mpi_ranks )
# Modify domain region
zmin, zmax, Nz = self.comm.divide_into_domain()
Nr = self.comm.get_Nr( with_damp=True )
rmax = self.comm.get_rmax( with_damp=True )
# Initialize the field structure
self.fld = Fields( Nz, zmax, Nr, rmax, Nm, dt,
n_order=n_order, zmin=zmin,
v_comoving=v_comoving,
use_pml=self.use_pml,
use_galilean=use_galilean,
current_correction=current_correction,
use_cuda=self.use_cuda,
smoother=smoother,
# Only create threading buffers when running on CPU
create_threading_buffers=(self.use_cuda is False) )
# Initialize the electrons and the ions
self.grid_shape = self.fld.interp[0].Ez.shape
self.particle_shape = particle_shape
self.ptcl = []
if n_e is not None:
# - Initialize the electrons
# (This is done in preparation for gathering among procs)
saved_Er = []
saved_Et = []
for m in range(sim.fld.Nm):
saved_Er.append( sim.fld.interp[m].Er.copy() )
sim.fld.interp[m].Er[:,:] = laser_Er[:,:,m]
saved_Et.append( sim.fld.interp[m].Et.copy() )
sim.fld.interp[m].Et[:,:] = laser_Et[:,:,m]
# Create a global field object across all subdomains, and copy the fields
# (Calculating the self-consistent Ez and B is a global operation)
global_Nz, _ = sim.comm.get_Nz_and_iz(
local=False, with_damp=True, with_guard=False )
global_zmin, global_zmax = sim.comm.get_zmin_zmax(
local=False, with_damp=True, with_guard=False )
global_fld = Fields( global_Nz, global_zmax,
sim.fld.Nr, sim.fld.rmax, sim.fld.Nm, sim.fld.dt,
zmin=global_zmin, n_order=sim.fld.n_order, use_cuda=False)
# Gather the fields of the interpolation grid
for m in range(sim.fld.Nm):
for field in ['Er', 'Et']:
local_array = getattr( sim.fld.interp[m], field )
gathered_array = sim.comm.gather_grid_array(
local_array, with_damp=True)
setattr( global_fld.interp[m], field, gathered_array )
# Now that the (gathered) laser fields are stored in global_fld,
# copy the saved field back into the local grid
for m in range(sim.fld.Nm):
sim.fld.interp[m].Er[:,:] = saved_Er[m]
sim.fld.interp[m].Et[:,:] = saved_Et[m]