Simulation of Electromagnetic Waves in Python with Meep

In this blog, we will walk through the process of simulating electromagnetic waves in a 3D cube using Meep, a popular finite-difference time-domain (FDTD) simulation software. The goal of this simulation is to visualize electromagnetic field distributions, flux through different faces of the cube, and the field intensity at various time steps. We will also introduce the concept of adding sources, dielectric materials, and flux monitors to track the flux through the cube.

1. Setting Up the Environment and Installing Required Packages

Before diving into the simulation, we need to ensure that all the required libraries and environments are set up. We first check if conda is installed. If it’s not, we install condacolab and set up a conda environment. We then check for the required packages, including pymeep, nlopt, and others, to ensure the simulation runs smoothly.

import os
import subprocess

# Function to check if a package is installed
def is_package_installed(package_name):
    try:
        subprocess.check_call([package_name, '--version'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
        return True
    except (subprocess.CalledProcessError, FileNotFoundError):
        return False

# Function to check if a conda environment exists
def is_conda_env_installed(env_name):
    try:
        output = subprocess.check_output(['conda', 'env', 'list'], stderr=subprocess.STDOUT).decode()
        return env_name in output
    except subprocess.CalledProcessError:
        return False

# Install condacolab if conda is not available
if not is_package_installed('conda'):
    subprocess.check_call(['pip', 'install', '-q', 'condacolab'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
    import condacolab
    condacolab.install()

# Check again if conda is installed after condacolab installation
if is_package_installed('conda'):
    # Install pymeep using mamba if not already installed
    try:
        import meep as mp
    except ImportError:
        subprocess.check_call(['mamba', 'install', '-c', 'conda-forge', 'pymeep', '-y'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)

    # Create the parallel Meep environment if it doesn't exist
    if not is_conda_env_installed('pmp'):
        subprocess.check_call(['mamba', 'create', '-n', 'pmp', '-c', 'conda-forge', 'pymeep=*=mpi_mpich_*', '-y'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)

    # Check if nlopt is installed, install if not
    try:
        import nlopt
    except ImportError:
        subprocess.check_call(['pip', 'install', '-q', 'nlopt'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)

    # Check if everything is installed correctly
    subprocess.check_call(['conda', 'list'], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
else:
    raise EnvironmentError("Conda installation failed")

2. Defining the Simulation Parameters and Geometry

In the next step, we define the parameters for the simulation, including the size of the cube, the refractive index of the medium, and the frequency range of the sources. We also set up the geometry of the cube, where we place sources at the corners of the cube and define a spherical object in the center. The cube’s sides are defined by LL, and the sources are distributed across the corners.

import meep as mp
import numpy as np
import matplotlib.pyplot as plt

# Simulation parameters
resolution = 10
n_r = 1.45  # refractive index of the medium
lambda_min = 0.4
lambda_max = 0.8
f_min = 1 / lambda_max
f_max = 1 / lambda_min
fcen = (f_min + f_max) / 2  # center frequency
df = f_max - f_min  # frequency width
nfreq = 40  # number of frequency points
timesteps = 100  # Number of time steps to collect

# Define the material for the object inside the cube
material = mp.Medium(index=n_r)

# Define the cube cell
LL = 10  # Length of the cube sides
L = LL / 2
cell = mp.Vector3(LL, LL, LL)  # Simulation domain size

# Define sources at each corner of the cube with varying frequencies
source_frequencies = [fcen, fcen + df/4, fcen - df/4, fcen + df/8, fcen - df/8, fcen, fcen + df/2, fcen - df/2]
source_locs = [
    mp.Vector3(-L / 2, -L / 2, -L / 2), mp.Vector3(L / 2, -L / 2, -L / 2),
    mp.Vector3(-L / 2, L / 2, -L / 2), mp.Vector3(L / 2, L / 2, -L / 2),
    mp.Vector3(-L / 2, -L / 2, L / 2), mp.Vector3(L / 2, -L / 2, L / 2),
    mp.Vector3(-L / 2, L / 2, L / 2), mp.Vector3(L / 2, L / 2, L / 2)
]

sources = [mp.Source(
    src=mp.ContinuousSource(frequency=freq),
    component=mp.Ez,
    center=loc
) for freq, loc in zip(source_frequencies, source_locs)]

# Define the geometry (e.g., add a sphere in the center of the cube)
geometry = [mp.Sphere(radius=L / 2, center=mp.Vector3(0, 0, 0), material=material)]

# Create the simulation
sim = mp.Simulation(
    cell_size=cell,
    geometry=geometry,
    sources=sources,
    resolution=resolution
)

3. Flux Regions and Monitoring

Flux regions are defined to monitor the flux through the faces of the cube. We create six flux regions, each corresponding to a face of the cube. The flux is recorded as the simulation runs, allowing us to visualize the total flux passing through each face over time.

# Define FluxRegions on the faces

flux_planes = [
    mp.FluxRegion(center=mp.Vector3(-L / 2, 0, 0), size=mp.Vector3(0, LL, LL)),  # Left face
    mp.FluxRegion(center=mp.Vector3(L / 2, 0, 0), size=mp.Vector3(0, LL, LL)),   # Right face
    mp.FluxRegion(center=mp.Vector3(0, -L / 2, 0), size=mp.Vector3(LL, 0, LL)),  # Bottom face
    mp.FluxRegion(center=mp.Vector3(0, L / 2, 0), size=mp.Vector3(LL, 0, LL)),   # Top face
    mp.FluxRegion(center=mp.Vector3(0, 0, -L / 2), size=mp.Vector3(LL, LL, 0)),  # Back face
    mp.FluxRegion(center=mp.Vector3(0, 0, L / 2), size=mp.Vector3(LL, LL, 0))    # Front face
]

# Add flux monitors for the selected planes
flux_monitors = [sim.add_flux(fcen, df, nfreq, plane) for plane in flux_planes]

# Initialize storage for total flux data (only storing summed flux)
flux_storage = {i: [] for i in range(len(flux_planes))}  # Dictionary to store flux data per face

# Collect field data less frequently (every 5th time step)
field_data = []
field_intensity = []
# Function to collect total flux data at each time step (sum over frequencies)

def collect_flux(sim):
    center_field = sim.get_field_point(mp.Ez, mp.Vector3(0, 0, 0))
    field_intensity.append(center_field**2)  # Store |E|^2

    field = sim.get_array(center=mp.Vector3(0, 0, 0), size=cell, component=mp.Ez)
    field_data.append(field)
    for i, monitor in enumerate(flux_monitors):
        flux_data = sim.get_flux_data(monitor)  # Get flux data for each monitor
        total_flux = np.sum(flux_data)  # Sum flux over all frequencies at this time step
        flux_storage[i].append(total_flux)  # Store summed flux data for each face

# Run the simulation and collect flux data every 5th time step
sim.run(mp.at_every(5, collect_flux), until=timesteps)

4. Visualization of Source Locations, Flux, and Field Intensity

We now proceed to visualize the simulation results. First, we plot the source locations within the cube. Then, we visualize the flux through each face of the cube and the field intensity at the center of the sphere. Additionally, we can inspect the field data at different time steps.

Source Locations

# Extract source coordinates
source_coords = np.array([[loc.x, loc.y, loc.z] for loc in source_locs])

# Plot source locations in 3D
fig = plt.figure(dpi=100)
ax = fig.add_subplot(111, projection='3d')

# Plot cube boundaries
ax.set_xlim(-L / 2, L / 2)
ax.set_ylim(-L / 2, L / 2)
ax.set_zlim(-L / 2, L / 2)
ax.set_box_aspect([1, 1, 1])  # Equal aspect ratio

# Add sources as red markers
ax.scatter(source_coords[:, 0], source_coords[:, 1], source_coords[:, 2],
           color='red', label='Sources', s=100)

# Add labels and legend
ax.set_title('Source Locations in Cube')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
ax.legend()

plt.show()

Flux Over Time

# Visualize the total flux over time for each face
for i, flux_data in flux_storage.items():
    total_flux_over_time = np.array(flux_data)  # Convert to numpy array for easier manipulation
    actual_time_steps = np.arange(len(total_flux_over_time))  # Time steps array
    plt.plot(actual_time_steps, total_flux_over_time, label=f'Face {i + 1}')

plt.xlabel("Time Step")
plt.ylabel("Total Flux")
plt.title("Total Flux Through Cube Faces Over Time")
plt.legend()
plt.show()

Field Intensity

# Plot the field intensity over time
plt.figure(dpi=100)
plt.plot(actual_time_steps, field_intensity, label="Field Intensity at Sphere Center")
plt.xlabel("Time Step")
plt.ylabel("Field Intensity (|E|^2)")
plt.title("Field Intensity Reaching Sphere Center Over Time")
plt.legend()
plt.grid(True)
plt.show()

Field Data at Specific Time Steps

# Convert field data into a numpy array for easier handling
field_data = np.array(field_data)

# Visualize field data at specific time steps
time_indices = [0, 5, 10]

for t in time_indices:
    ez_data = field_data[t]
    z_slice = ez_data.shape[2] // 2  # Take a slice along the z-axis

    plt.figure(dpi=100)
    plt.imshow(ez_data[:, :, z_slice], interpolation='spline36', cmap='inferno')
    plt.colorbar(label='Ez field intensity')
    plt.title(f'Ez Field Distribution at Time Step {t}')
    plt.axis('off')
    plt.show()

Dielectric Constant Visualization

# Get the dielectric constant array across the entire simulation domain
eps_data = sim.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric)

# Select a specific slice along the z-axis (the middle slice of the cube)
z_slice = eps_data.shape[2] // 2  # Selecting the middle slice along the z-axis

# Plot the selected slice
plt.figure(dpi=100)
plt.imshow(eps_data[:, :, z_slice].transpose(), interpolation='spline36', cmap='binary')
plt.colorbar(label="Dielectric Constant")
plt.title(f"Dielectric Constant Slice at z = {z_slice - eps_data.shape[2]//2}")
plt.axis('off')
plt.show()

Conclusion

In this post, we demonstrated how to set up a 3D electromagnetic wave simulation using Meep, with sources placed at the cube’s corners and a spherical object placed at the center. We monitored the flux through the faces of the cube and visualized the results, including field intensity at various time steps and dielectric constant distribution. This type of simulation can be expanded to study more complex scenarios involving various materials and geometries.