Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Minor changes to docstrings.
Browse files Browse the repository at this point in the history
lucianogsilvestri committed Mar 29, 2024
1 parent 72ac984 commit d491179
Showing 2 changed files with 15 additions and 75 deletions.
27 changes: 7 additions & 20 deletions sarkas/potentials/force_pm.py
Original file line number Diff line number Diff line change
@@ -32,29 +32,24 @@ def assgnmnt_func(cao, x):
W = zeros(cao)

if cao == 1:

W[0] = 1.0

elif cao == 2:

W[0] = 0.5 * (1.0 - 2.0 * x)
W[1] = 0.5 * (1.0 + 2.0 * x)

elif cao == 3:

W[0] = (1.0 - 4.0 * x + 4.0 * x**2) / 8.0
W[1] = (3.0 - 4.0 * x**2) / 4.0
W[2] = (1.0 + 4.0 * x + 4.0 * x**2) / 8.0

elif cao == 4:

W[0] = (1.0 - 6.0 * x + 12.0 * x**2 - 8.0 * x**3) / 48.0
W[1] = (23.0 - 30.0 * x - 12.0 * x**2 + 24.0 * x**3) / 48.0
W[2] = (23.0 + 30.0 * x - 12.0 * x**2 - 24.0 * x**3) / 48.0
W[3] = (1.0 + 6.0 * x + 12.0 * x**2 + 8.0 * x**3) / 48.0

elif cao == 5:

W[0] = (1.0 - 8.0 * x + 24.0 * x**2 - 32.0 * x**3 + 16.0 * x**4) / 384.0
W[1] = (19.0 - 44.0 * x + 24.0 * x**2 + 16.0 * x**3 - 16.0 * x**4) / 96.0
W[2] = (115.0 - 120.0 * x**2 + 48.0 * x**4) / 192.0
@@ -70,7 +65,6 @@ def assgnmnt_func(cao, x):
W[5] = (1.0 + 10.0 * x + 40.0 * x**2 + 80.0 * x**3 + 80.0 * x**4 + 32.0 * x**5) / 3840.0

elif cao == 7:

W[0] = (
1.0 - 12.0 * x + 60.0 * x * 2 - 160.0 * x**3 + 240.0 * x**4 - 192.0 * x**5 + 64.0 * x**6
) / 46080.0
@@ -166,9 +160,8 @@ def calc_acc_pm(E_x_r, E_y_r, E_z_r, mesh_pos, mesh_points, q_over_m, cao, mesh_
acc = zeros_like(mesh_pos)

for ipart, q_m in enumerate(q_over_m):

ix = mesh_points[ipart, 0]
x = mesh_pos[ipart,0 ] - (ix + mid[0])
x = mesh_pos[ipart, 0] - (ix + mid[0])

iy = mesh_points[ipart, 1]
y = mesh_pos[ipart, 1] - (iy + mid[1])
@@ -196,7 +189,6 @@ def calc_acc_pm(E_x_r, E_y_r, E_z_r, mesh_pos, mesh_points, q_over_m, cao, mesh_
iyn = iy - pshift[1] # min. index along y-axis

for i in range(cao[1]):

# if iyn < 0:
# r_i = iyn + mesh_sz[1]
# elif iyn > (mesh_sz[1] - 1):
@@ -228,9 +220,9 @@ def calc_acc_pm(E_x_r, E_y_r, E_z_r, mesh_pos, mesh_points, q_over_m, cao, mesh_

izn += 1

acc[:,0] = E_x_p
acc[:,1] = E_y_p
acc[:,2] = E_z_p
acc[:, 0] = E_x_p
acc[:, 1] = E_y_p
acc[:, 2] = E_z_p

return acc

@@ -276,7 +268,6 @@ def calc_charge_dens(mesh_pos, mesh_points, charges, cao, mesh_sz, mid, pshift):
# (ix + 0.5)*h_array[0] = midpoint between the two mesh points closest to the particle

for ipart in range(len(charges)):

ix = mesh_points[ipart, 0]
delta_x = mesh_pos[ipart, 0] - (ix + mid[0])

@@ -294,7 +285,6 @@ def calc_charge_dens(mesh_pos, mesh_points, charges, cao, mesh_sz, mid, pshift):
izn = iz - pshift[2] # min. index along z-axis

for g in range(cao[2]):

# if izn < 0:
# r_g = izn + mesh_sz[2]
# elif izn > (mesh_sz[2] - 1):
@@ -306,7 +296,6 @@ def calc_charge_dens(mesh_pos, mesh_points, charges, cao, mesh_sz, mid, pshift):
iyn = iy - pshift[1] # min. index along y-axis

for i in range(cao[1]):

r_i = iyn + mesh_sz[1] * (iyn < 0) - mesh_sz[1] * (iyn > (mesh_sz[1] - 1))

# if iyn < 0:
@@ -472,7 +461,6 @@ def calc_pot_pm(phi_r, mesh_pos, mesh_points, charges, cao, mesh_sz, mid, pshift
pot_p = zeros_like(charges) # potential energy for

for ipart, q in enumerate(charges):

ix = mesh_points[ipart, 0]
x = mesh_pos[ipart, 0] - (ix + mid[0])

@@ -502,7 +490,6 @@ def calc_pot_pm(phi_r, mesh_pos, mesh_points, charges, cao, mesh_sz, mid, pshift
iyn = iy - pshift[1] # min. index along y-axis

for i in range(cao[1]):

# if iyn < 0:
# r_i = iyn + mesh_sz[1]
# elif iyn > (mesh_sz[1] - 1):
@@ -969,9 +956,9 @@ def update(pos, charges, masses, mesh_sizes, mesh_spacings, mesh_volume, box_vol
phi_k = G_k * rho_k

# Charge density
rho_k_real = rho_k.real
rho_k_imag = rho_k.imag
rho_k_sq = rho_k_real * rho_k_real + rho_k_imag * rho_k_imag
# rho_k_real = rho_k.real
# rho_k_imag = rho_k.imag
# rho_k_sq = rho_k_real * rho_k_real + rho_k_imag * rho_k_imag

# Calculate the Electric field's component on the mesh
E_kx, E_ky, E_kz = calc_field(phi_k, kx_v, ky_v, kz_v)
63 changes: 8 additions & 55 deletions sarkas/time_evolution/integrators.py
Original file line number Diff line number Diff line change
@@ -267,7 +267,6 @@ def pot_acc_setup(self, potential):
)

def boundary_condition_setup(self):

self.supported_boundary_conditions = {
"periodic": self.periodic_bc,
"absorbing": self.absorbing_bc,
@@ -332,13 +331,11 @@ def type_setup(self, int_type):
# Assign integrator.update to the correct method

if int_type == "langevin":

self.sigma = sqrt(2.0 * self.langevin_gamma * self.kB * self.thermostat_temperatures / self.species_masses)
self.c1 = 1.0 - 0.5 * self.langevin_gamma * self.dt
self.c2 = 1.0 / (1.0 + 0.5 * self.langevin_gamma * self.dt)

elif int_type == "magnetic_verlet":

# Calculate functions for magnetic integrator
# This could be used when the generalization to Forest-Ruth and MacLachlan algorithms will be implemented
# In a magnetic Velocity-Verlet the coefficient is 1/2, see eq.~(78) in :cite:`Chin2008`
@@ -357,7 +354,6 @@ def type_setup(self, int_type):
int_type = "magnetic_pos_verlet_zdir"

elif int_type == "magnetic_boris":

# In a leapfrog-type algorithm the coefficient is different for the acceleration and magnetic rotation
# see eq.~(79) in :cite:`Chin2008`
self.magnetic_helpers(1.0)
@@ -508,11 +504,6 @@ def magnetic_verlet_zdir(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
Notes
-----
This integrator is faster than `magnetic_verlet` but valid only for a magnetic field in the :math:`z`-direction.
@@ -587,11 +578,6 @@ def magnetic_verlet(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
Notes
-----
:cite:`Chin2008` equations are written for a negative charge. This allows him to write
@@ -655,11 +641,6 @@ def magnetic_boris_zdir(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
"""
# First half step of velocity update: Apply exp(dt * V_F / 2)
ptcls.vel += 0.5 * ptcls.acc * self.dt
@@ -697,11 +678,6 @@ def magnetic_boris(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
"""

# First half step of velocity update: Apply exp(eV_F/2)
@@ -736,11 +712,6 @@ def magnetic_pos_verlet_zdir(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
Notes
-----
This integrator is faster than `magnetic_verlet` but valid only for a magnetic field in the :math:`z`-direction.
@@ -797,11 +768,6 @@ def magnetic_pos_verlet(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
Notes
-----
:cite:`Chin2008` equations are written for a negative charge. This allows him to write
@@ -855,11 +821,6 @@ def cyclotronic_zdir(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
"""
# Drift half step
# Rotate Positions
@@ -914,11 +875,6 @@ def cyclotronic(self, ptcls):
ptcls: :class:`sarkas.particles.Particles`
Particles data.
Returns
-------
potential_energy : float
Total potential energy.
"""
# Drift half step

@@ -1065,9 +1021,9 @@ def pretty_print(self):
integrator_msg += time_msg
if self.potential_type == "qsp":
wp_e = self.species_plasma_frequencies[0]
t_we = (2.0 * pi)/wp_e
t_we = (2.0 * pi) / wp_e
wp_ions = norm(self.species_plasma_frequencies[1:])
t_wi = (2.0 * pi)/wp_ions
t_wi = (2.0 * pi) / wp_ions
qsp_msg = (
f"e plasma frequency = {wp_e:.6e} {self.units_dict['frequency']}\n"
f"w_pe dt = {self.dt * wp_e:.2e} [rad]\n"
@@ -1079,18 +1035,17 @@ def pretty_print(self):
f"ions plasma period (T_i) = {t_wi:.6e} {self.units_dict['time']}\n"
f"dt/T_i = {self.dt/t_wi:.2e}\n"
f"Timesteps per ion plasma cycle ~ {int(t_wi/self.dt)}\n"

)
integrator_msg += qsp_msg

elif self.potential_type == "lj":
integrator_msg += f"The plasma frequency is defined as w_p = sqrt( epsilon / (sigma^2 * mass) )\n"

if self.magnetized:
high_wc = abs(self.species_cyclotron_frequencies).max()
low_wc = abs(self.species_cyclotron_frequencies).min()
t_wc_high = 2.0 * pi/high_wc
t_wc_low = 2.0 * pi/low_wc
high_wc = abs(self.species_cyclotron_frequencies).max()
low_wc = abs(self.species_cyclotron_frequencies).min()
t_wc_high = 2.0 * pi / high_wc
t_wc_low = 2.0 * pi / low_wc
mag_msg = "\nMagnetic Timescales:\n"
if high_wc > low_wc:
mag_msg += (
@@ -1104,7 +1059,7 @@ def pretty_print(self):
f"Cyclotron period (T) = {t_wc_low:.6e} {self.units_dict['time']}\n"
f"dt/T = = {self.dt/t_wc_low:.2e}\n"
f"Timesteps per plasma cycle ~ {int(t_wc_low/self.dt)}\n"
)
)
else:
mag_msg += (
f"Cyclotron frequency (w_c) = {high_wc:.6e} {self.units_dict['frequency']}\n"
@@ -1131,6 +1086,7 @@ def pretty_print(self):
msg += integrator_msg
print(msg)


@jit(void(float64[:, :], float64[:], float64[:], int64[:], float64), nopython=True)
def berendsen(vel, T_desired, T, species_np, tau):
"""
@@ -1192,7 +1148,6 @@ def enforce_pbc(pos, cntr, box_vector):
# Loop over all particles
for p in arange(pos.shape[0]):
for d in arange(pos.shape[1]):

# If particle is outside of box in positive direction, wrap to negative side
# if pos[d,p] > box_vector[d]:
pos[p, d] -= box_vector[d] * (pos[p, d] > box_vector[d])
@@ -1230,7 +1185,6 @@ def enforce_abc(pos, vel, acc, charges, box_vector):
# Loop over all particles
for p in arange(pos.shape[0]):
for d in arange(pos.shape[1]):

# If particle is outside of box in positive direction, remove charge, velocity and acceleration
if pos[p, d] >= box_vector[d]:
pos[p, d] = box_vector[d]
@@ -1269,7 +1223,6 @@ def enforce_rbc(pos, vel, box_vector, dt):
# Loop over all particles
for p in arange(pos.shape[0]):
for d in arange(pos.shape[1]):

# If particle is outside of box in positive direction, wrap to negative side
if pos[p, d] > box_vector[d] or pos[p, d] < 0.0:
# Revert velocity

0 comments on commit d491179

Please sign in to comment.