diff --git a/sarkas/potentials/force_pm.py b/sarkas/potentials/force_pm.py index 7392424d..a9329ddd 100644 --- a/sarkas/potentials/force_pm.py +++ b/sarkas/potentials/force_pm.py @@ -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) diff --git a/sarkas/time_evolution/integrators.py b/sarkas/time_evolution/integrators.py index 65c68d97..18c1b95e 100644 --- a/sarkas/time_evolution/integrators.py +++ b/sarkas/time_evolution/integrators.py @@ -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,7 +1035,6 @@ 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 @@ -1087,10 +1042,10 @@ def pretty_print(self): 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