simulation module

Module contains the Simulation class and helper functions for save files.

The Simulation class is used to calculate the electromagnetic fields and potentials generated by moving point charges, which are defined by the Charge base class. The values are determined numerically at each grid point by evaluating the Liénard–Wiechert potentials and corresponding electromagnetic field equations at the retarded time.

The Simulation class also runs simulations using the Dipole class by dynamically calculating the dipole moment at each time step using the run method. Treating the dipole as a Lorentz oscillator, the differential equation of motion is solved at each time step using the RK4 method. The driving electric field acting on the dipole is the component of the external electric field in the direction of polarization generated by the other point charges in the simulation. After the simulation is run, the Simulation object can calculate the electromagnetic fields and potentials generated by the Dipole objects within the timespan of the simulation.

All units are in SI.

`Simulation`

Primary class for electromagnetism calculations using PyCharge.

The Simulation class calculates the electromagnetic fields and potentials generated by point charges that are subclasses of the Charge base class and have predefined trajectories.

The Simulation class also run simulations with Dipole objects over a range of time steps to dynamically determine their dipole moment. The dipoles behave as Lorentz oscillators and their dipole moment and related trajectories are calculated at each time step.

Parameters:

Name	Type	Description	Default
`sources`	`Union[Sequence, Union[Charge, Dipole]]`	Either a single or list of `Dipole` and `Charge` object(s) in the simulation.	required
`tol`	`float`	Tolerance for the scipy.optimize.newton method. Default is `1e-22`.	required

`calculate_A(self, t, x, y, z, exclude_charges=None)`

Calculate the vector potential (A) generated by the point charge(s).

Parameters:

Name	Type	Description	Default
`t`	`float`	Time of simulation.	required
`x`	`ndarray`	Meshgrid of x values.	required
`y`	`ndarray`	Meshgrid of y values.	required
`z`	`ndarray`	Meshgrid of z values.	required
`exclude_charges`	`Sequence`	List of charges to exclude in calculation of the potentials. Defaults to None.	`None`

Returns:

Type	Description
`Tuple[ndarray, ndarray, ndarray]`	List of Ax, Ay, and Az.

Source code in pycharge/simulation.py

def calculate_A(
    self,
    t: float,
    x: ndarray,
    y: ndarray,
    z: ndarray,
    exclude_charges: Sequence = None
) -> Tuple[ndarray, ndarray, ndarray]:
    """Calculate the vector potential (A) generated by the point charge(s).

    Args:
        t: Time of simulation.
        x: Meshgrid of x values.
        y: Meshgrid of y values.
        z: Meshgrid of z values.
        exclude_charges: List of charges to exclude in calculation of the
            potentials. Defaults to None.

    Returns:
        List of Ax, Ay, and Az.
    """
    t_array = np.ones((x.shape))
    t_array[:, :, :] = t
    Ax = np.zeros((x.shape))
    Ay = np.zeros((x.shape))
    Az = np.zeros((x.shape))
    for charge in self.all_charges:
        if exclude_charges is not None and charge in exclude_charges:
            continue
        initial_guess = -1e-12*np.ones((x.shape))
        tr = optimize.newton(charge.solve_time, initial_guess,
                             args=(t_array, x, y, z), tol=self.tol)
        rx = x - charge.xpos(tr)
        ry = y - charge.ypos(tr)
        rz = z - charge.zpos(tr)
        r_mag = (rx**2 + ry**2 + rz**2)**0.5
        vx = charge.xvel(tr)  # retarded velocity - Griffiths Eq. 10.54
        vy = charge.yvel(tr)
        vz = charge.zvel(tr)
        r_dot_v = rx*vx + ry*vy + rz*vz
        # Griffiths Eq. 10.53
        individual_V = charge.q*c/(4*pi*eps*(r_mag*c-r_dot_v))
        # Griffiths Eq. 10.53
        Ax += vx/c**2*individual_V
        Ay += vy/c**2*individual_V
        Az += vz/c**2*individual_V
    return (Ax, Ay, Az)

`calculate_B(self, t, x, y, z, pcharge_field='Total', exclude_charges=None)`

Calculate the magnetic field (B) generated by the point charge(s).

Parameters:

Name	Type	Description	Default
`t`	`float`	Time of simulation.	required
`x`	`ndarray`	Meshgrid of x values.	required
`y`	`ndarray`	Meshgrid of y values.	required
`z`	`ndarray`	Meshgrid of z values.	required
`pcharge_field`	`str`	Determines which field is returned: `Velocity`, `Acceleration`, or `Total`. Defaults to `Total`.	`'Total'`
`exclude_charges`	`Sequence`	List of charges to exclude in calculation of the magnetic field. Defaults to None.	`None`

Returns:

Type	Description
`ndarray`	List of Bx, By, and Bz.

Source code in pycharge/simulation.py

def calculate_B(
    self,
    t: float,
    x: ndarray,
    y: ndarray,
    z: ndarray,
    pcharge_field: str = 'Total',
    exclude_charges: Sequence = None
) -> ndarray:
    """Calculate the magnetic field (B) generated by the point charge(s).

    Args:
        t: Time of simulation.
        x: Meshgrid of x values.
        y: Meshgrid of y values.
        z: Meshgrid of z values.
        pcharge_field: Determines which field is returned:
            `Velocity`, `Acceleration`, or `Total`. Defaults to `Total`.
        exclude_charges: List of charges to exclude in calculation of the
            magnetic field. Defaults to None.

    Returns:
        List of Bx, By, and Bz.
    """
    t_array = np.ones((x.shape))
    t_array[:, :, :] = t
    Bx = np.zeros((x.shape))
    By = np.zeros((x.shape))
    Bz = np.zeros((x.shape))
    for charge in self.all_charges:
        if exclude_charges is not None and charge in exclude_charges:
            continue
        initial_guess = -1e-12*np.ones((x.shape))
        tr = optimize.newton(charge.solve_time, initial_guess,
                             args=(t_array, x, y, z), tol=self.tol)
        E_x, E_y, E_z = self._calculate_individual_E(
            charge, tr, x, y, z, pcharge_field)
        rx = x - charge.xpos(tr)
        ry = y - charge.ypos(tr)
        rz = z - charge.zpos(tr)
        r_mag = (rx**2 + ry**2 + rz**2)**0.5
        Bx += 1/(c*r_mag)*(ry*E_z-rz*E_y)
        By += 1/(c*r_mag)*(rz*E_x-rx*E_z)
        Bz += 1/(c*r_mag)*(rx*E_y-ry*E_x)
    return np.array((Bx, By, Bz))

`calculate_E(self, t, x, y, z, pcharge_field='Total', exclude_charges=None)`

Calculate the electric field (E) generated by the point charge(s).

Meshgrid must use matrix indexing ('ij').

Parameters:

Name	Type	Description	Default
`t`	`float`	Time of simulation.	required
`x`	`ndarray`	Meshgrid of x values.	required
`y`	`ndarray`	Meshgrid of y values.	required
`z`	`ndarray`	Meshgrid of z values.	required
`pcharge_field`	`str`	Determines which field is returned: `Velocity`, `Acceleration`, or `Total`. Defaults to `Total`.	`'Total'`
`exclude_charges`	`Sequence`	List of charges to exclude in calculation of the electric field. Defaults to None.	`None`

Returns:

Type	Description
`ndarray`	List of E_x, E_y, and E_z.

Source code in pycharge/simulation.py

def calculate_E(
    self,
    t: float,
    x: ndarray,
    y: ndarray,
    z: ndarray,
    pcharge_field: str = 'Total',
    exclude_charges: Sequence = None
) -> ndarray:
    """Calculate the electric field (E) generated by the point charge(s).

    Meshgrid must use matrix indexing ('ij').

    Args:
        t: Time of simulation.
        x: Meshgrid of x values.
        y: Meshgrid of y values.
        z: Meshgrid of z values.
        pcharge_field: Determines which field is returned: `Velocity`,
            `Acceleration`, or `Total`. Defaults to `Total`.
        exclude_charges: List of charges to exclude in calculation of the
            electric field. Defaults to None.

    Returns:
        List of E_x, E_y, and E_z.
    """
    t_array = np.ones((x.shape))
    t_array[:, :, :] = t
    E_x = np.zeros((x.shape))
    E_y = np.zeros((x.shape))
    E_z = np.zeros((x.shape))
    for charge in self.all_charges:
        if exclude_charges is not None and charge in exclude_charges:
            continue
        initial_guess = -1e-12*np.ones((x.shape))
        tr = optimize.newton(charge.solve_time, initial_guess,
                             args=(t_array, x, y, z), tol=self.tol)
        E_field = self._calculate_individual_E(
            charge, tr, x, y, z, pcharge_field)
        E_x += E_field[0]
        E_y += E_field[1]
        E_z += E_field[2]
    return np.array((E_x, E_y, E_z))

`calculate_V(self, t, x, y, z, exclude_charges=None)`

Calculate the scalar potential (V) generated by the point charge(s).

Parameters:

Name	Type	Description	Default
`t`	`float`	Time of simulation.	required
`x`	`ndarray`	Meshgrid of x values.	required
`y`	`ndarray`	Meshgrid of y values.	required
`z`	`ndarray`	Meshgrid of z values.	required
`exclude_charges`	`Sequence`	List of charges to exclude in calculation of the potentials. Defaults to None.	`None`

Returns:

Type	Description
`ndarray`	Scalar potential V.

Source code in pycharge/simulation.py

def calculate_V(
    self,
    t: float,
    x: ndarray,
    y: ndarray,
    z: ndarray,
    exclude_charges: Sequence = None
) -> ndarray:
    """Calculate the scalar potential (V) generated by the point charge(s).

    Args:
        t: Time of simulation.
        x: Meshgrid of x values.
        y: Meshgrid of y values.
        z: Meshgrid of z values.
        exclude_charges: List of charges to exclude in calculation of the
            potentials. Defaults to None.

    Returns:
        Scalar potential V.
    """
    t_array = np.ones((x.shape))
    t_array[:, :, :] = t
    V = np.zeros((x.shape))
    for charge in self.all_charges:
        if exclude_charges is not None and charge in exclude_charges:
            continue
        initial_guess = -1e-12*np.ones((x.shape))
        tr = optimize.newton(charge.solve_time, initial_guess,
                             args=(t_array, x, y, z), tol=self.tol)
        rx = x - charge.xpos(tr)
        ry = y - charge.ypos(tr)
        rz = z - charge.zpos(tr)
        r_mag = (rx**2 + ry**2 + rz**2)**0.5
        vx = charge.xvel(tr)  # retarded velocity - Griffiths Eq. 10.54
        vy = charge.yvel(tr)
        vz = charge.zvel(tr)
        r_dot_v = rx*vx + ry*vy + rz*vz
        # Griffiths Eq. 10.53
        V += charge.q*c/(4*pi*eps*(r_mag*c-r_dot_v))
    return V

`run(self, timesteps, dt, file=None, save_E=False, max_vel=2997924.58, progressbar=True)`

Run simulation with Dipole and Charge object(s).

Simulation calculates the dipole moment and corresponding derivatives at each time step by solving the equation of motion using the RK4 method. The position, velocity, and acceleration values are then updated in the Dipole object.

Parameters:

Name	Type	Description	Default
`timesteps`	`int`	Number of time steps run during the simulation.	required
`dt`	`float`	Time step.	required
`file`	`Optional[str]`	File name to load/save `Simulation`. Defaults to `None`.	`None`
`save_E`	`bool`	If True, the driving field acting on each dipole is saved at each time step. However, this comes at a cost to speed and memory. Defaults to `True`.	`False`
`max_vel`	`float`	The maximum possible velocity achieved by the two point charges in `Dipole` object. Defaults to `c/100`.	`2997924.58`
`progressbar`	`bool`	Prints the progressbar if `True`. Defaults to `True`.	`True`

Exceptions:

Type	Description
`ValueError`	Raised if the point charge's velocity in the `Dipole` object becomes larger than `max_vel`.

Source code in pycharge/simulation.py

def run(
    self,
    timesteps: int,
    dt: float,
    file: Optional[str] = None,
    save_E: bool = False,
    max_vel: float = c/100,
    progressbar: bool = True
) -> None:
    """Run simulation with `Dipole` and `Charge` object(s).

    Simulation calculates the dipole moment and corresponding derivatives
    at each time step by solving the equation of motion using the RK4
    method. The position, velocity, and acceleration values are then
    updated in the `Dipole` object.

    Args:
        timesteps (int): Number of time steps run during the simulation.
        dt (float): Time step.
        file: File name to load/save `Simulation`. Defaults to `None`.
        save_E (bool): If True, the driving field acting on each dipole is
            saved at each time step. However, this comes at a cost to speed
            and memory. Defaults to `True`.
        max_vel: The maximum possible velocity achieved by the two point
            charges in `Dipole` object. Defaults to `c/100`.
        progressbar: Prints the progressbar if `True`. Defaults to `True`.

    Raises:
        ValueError: Raised if the point charge's velocity in the `Dipole`
            object becomes larger than `max_vel`.
    """
    self._init_simulation(timesteps, dt, save_E)
    if file is not None and self._load(file):
        return
    for tstep in tqdm(range(self.timesteps-1), total=self.timesteps,
                      initial=1, mininterval=0.5, disable=not progressbar):
        for dipole in self.dipoles:
            moment_disp, moment_vel = self._rk4(tstep, dipole)
            self._update_dipole(tstep, dipole, moment_disp, moment_vel,
                                max_vel)
    if file is not None:
        self._save(file)

`run_mpi(self, timesteps, dt, file=None, save_E=False, max_vel=2997924.58, progressbar=True)`

Run simulation with Dipole and Charge object(s) using MPI.

Method exploits MPI to speed up the simulation. In the ideal case, the number of MPI processors and the number of Dipole objects are equal. Most efficient simulation occurs when save_E is False.

Simulation calculates the dipole moment and corresponding derivatives at each time step by solving the equation of motion using the RK4 method. The position, velocity, and acceleration values are then updated in the Dipole object.

Parameters:

Name	Type	Description	Default
`timesteps`	`int`	Number of time steps run during the simulation.	required
`dt`	`float`	Time step.	required
`file`	`Optional[str]`	File name to load/save `Simulation`. Defaults to `None`.	`None`
`save_E`	`bool`	If True, the driving field acting on each dipole is saved at each time step. However, this comes at a cost to speed and memory. Defaults to `True`.	`False`
`max_vel`	`float`	The maximum possible velocity achieved by the two point charges in `Dipole` object. Defaults to `c/100`.	`2997924.58`
`progressbar`	`bool`	Prints the progressbar if `True`. Defaults to `True`.	`True`

Exceptions:

Type	Description
`NotImplementedError`	Raised if `mpi4py` package is not installed.
`ValueError`	Raised if the point charge's velocity in the `Dipole` object becomes larger than `max_vel`.

Examples:

An MPI program can be run from a script with the command mpiexec: >> $ mpiexec -n 4 python script.py Here, 4 separate processes will be initiated to run script.py.

Source code in pycharge/simulation.py

def run_mpi(
    self,
    timesteps: int,
    dt: float,
    file: Optional[str] = None,
    save_E: bool = False,
    max_vel: float = c/100,
    progressbar: bool = True
) -> None:
    """Run simulation with `Dipole` and `Charge` object(s) using MPI.

    Method exploits MPI to speed up the simulation. In the ideal case, the
    number of MPI processors and the number of `Dipole` objects are equal.
    Most efficient simulation occurs when `save_E` is `False`.

    Simulation calculates the dipole moment and corresponding derivatives
    at each time step by solving the equation of motion using the RK4
    method. The position, velocity, and acceleration values are then
    updated in the `Dipole` object.

    Args:
        timesteps (int): Number of time steps run during the simulation.
        dt (float): Time step.
        file: File name to load/save `Simulation`. Defaults to `None`.
        save_E (bool): If True, the driving field acting on each dipole is
            saved at each time step. However, this comes at a cost to speed
            and memory. Defaults to `True`.
        max_vel: The maximum possible velocity achieved by the two point
            charges in `Dipole` object. Defaults to `c/100`.
        progressbar: Prints the progressbar if `True`. Defaults to `True`.

    Raises:
        NotImplementedError: Raised if `mpi4py` package is not installed.
        ValueError: Raised if the point charge's velocity in the `Dipole`
            object becomes larger than `max_vel`.

    Example:
        An MPI program can be run from a script with the command mpiexec:
            >> $ mpiexec -n 4 python script.py
        Here, 4 separate processes will be initiated to run script.py.
    """
    try:
        from mpi4py import MPI  # pylint: disable=import-outside-toplevel
    except ImportError:
        NotImplementedError('mpi4py package is not installed.')
    self._init_simulation(timesteps, dt, save_E)
    if file is not None and self._load(file):
        return
    mpi_comm = MPI.COMM_WORLD
    mpi_size = mpi_comm.Get_size()  # Number of MPI processes
    mpi_rank = mpi_comm.Get_rank()  # Rank of current MPI process
    process_dipoles = self.dipoles[mpi_rank::mpi_size]  # Dipoles evaluated
    disable_bar = not progressbar or mpi_rank != 0
    for tstep in tqdm(range(self.timesteps-1), total=self.timesteps,
                      initial=1, mininterval=0.5, disable=disable_bar):
        moment_data = []
        for dipole in process_dipoles:  # Only evaluate `process_dipoles`
            moment_disp, moment_vel = self._rk4(tstep, dipole)
            self._update_dipole(tstep, dipole, moment_disp, moment_vel,
                                max_vel)
            moment_data.append((moment_disp, moment_vel))
        # Broadcast the dipole moment disp. and vel. to other MPI processes
        for bcast_rank in range(mpi_size):  # Send moment data
            if mpi_rank == bcast_rank:
                mpi_comm.bcast(moment_data, root=bcast_rank)
            else:  # Receive moment data and update dipoles
                bcast_data = None
                bcast_data = mpi_comm.bcast(bcast_data, root=bcast_rank)
                for dipole, data in zip(self.dipoles[bcast_rank::mpi_size],
                                        bcast_data):
                    self._update_dipole(tstep, dipole, data[0], data[1],
                                        max_vel)
    if mpi_rank == 0 and file is not None:
        self._save(file)

`combine_files(input_files, new_file)`

Combine Simulation objects from many .dat files into a single file.

Parameters:

Name	Type	Description	Default
`input_files`	`Tuple[str, ...]`	List of names of files to be concatenated.	required
`new_file`	`str`	Name of output file, should end in `.dat`.	required

Source code in pycharge/simulation.py

def combine_files(input_files: Tuple[str, ...], new_file: str) -> None:
    """Combine `Simulation` objects from many `.dat` files into a single file.

    Args:
        input_files: List of names of files to be concatenated.
        new_file: Name of output file, should end in `.dat`.
    """
    for file in input_files:
        with open(file, "rb") as f:
            loaded_simulation = pickle.load(f)
        with open(new_file, "ab") as new_f:
            pickle.dump(loaded_simulation, new_f)

`get_file_length(file)`

Return the number of Simulation objects in .dat file.

Parameters:

Name	Type	Description	Default
`file`	`str`	File name.	required

Returns:

Type	Description
`int`	Number of `Simulation` objects in file, `-1` if file not found.

Source code in pycharge/simulation.py

def get_file_length(file: str) -> int:
    """Return the number of `Simulation` objects in `.dat` file.

    Args:
        file: File name.

    Returns:
        Number of `Simulation` objects in file, `-1` if file not found.
    """
    length = 0
    try:
        with open(file, "rb") as f:
            while True:
                pickle.load(f)
                length += 1
    except EOFError:
        return length
    except FileNotFoundError:
        return -1