Initialize E field succeeds before sim.run but drops to 0 right after the first time step

[m-h-hsieh]

Hello everyone,

I'm a newbie to MEEP and am excited to apply MEEP to my research!

I'm doing a simple 1D simulation for an optical cavity consisting of two PECs with a thickness of one. I would like to initialize the E field and H field to values I pre-calculated. For now, it's just the fundamental mode for the cavity, i.e. a sine and a cosine for E and H, respectively.

I could initialize E with sim.fields.initialize_field() and got the desired E from sim.get_array(component=mp.Ex). However, once I started to run sim.run(), at the first time step, the E field becomes zero, while the H field is propagated as expected. The resulting field at the center of the cavity is plotted here. Note the sudden drop of Ex (blue line) right after the simulation starts.

The code is as follows:

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

l=12
res=10
Nz = res*l
dz=1/res
courant_factor = 0.1
dt = courant_factor / res

cell = mp.Vector3(0, 0, l)

# Define the left mirror (PEC)
left_mirror = mp.Block(
    size=mp.Vector3(mp.inf, mp.inf, 1),
    center=mp.Vector3(0,0,-(l-2)/2-1/2),
    material=mp.perfect_electric_conductor
)

# Define the right mirror (PEC)
right_mirror = mp.Block(
    size=mp.Vector3(mp.inf, mp.inf,1),
    center=mp.Vector3(0,0, (l-2)/2+1/2),
    material=mp.perfect_electric_conductor
)

sim = mp.Simulation(
    cell_size=cell,
    geometry=[left_mirror, right_mirror],
    resolution=res,
    Courant=courant_factor,
    dimensions=1
)

sim.init_sim()

L=Nz - 2*res
E = np.zeros(Nz)
E[res : -res] = np.sin(np.pi/L * np.arange(0, L))

def custom_E_field(p):
    shifted_x = p.z + l / 2
    index = int(shifted_x * res)
    index = max(0, min(Nz - 1, index))
    return E[index]

H = np.zeros(Nz)
H[res : -res] = np.cos(np.pi/L * np.arange(0.5, L + 0.5))

def custom_H_field(p):
    shifted_x = p.z + l / 2
    index = int(shifted_x * res)
    index = max(0, min(Nz - 1, index))
    return H[index]

sim.fields.initialize_field(mp.Ex, custom_E_field)
sim.fields.initialize_field(mp.Hy, custom_H_field)

# Function to store Ez field at each timestep
def store_field(sim):
    ez_snapshot = sim.get_array(component=mp.Ex)
    hy_snapshot = sim.get_array(component=mp.Hy)
    ez_data_list.append(ez_snapshot) 
    hy_data_list.append(hy_snapshot)

# List to store field snapshots at each timestep
ez_data_list = []
ez_data_list.append(sim.get_array(component=mp.Ex))
hy_data_list = []
hy_data_list.append(sim.get_array(component=mp.Hy))

sim.run(mp.at_every(dt, store_field), until=10)

ez_data_save = np.array(ez_data_list)
hy_data_save = np.array(hy_data_list)

plt.plot(ez_data_save[:,50], label='Ex')
plt.plot(hy_data_save[:,50], label='Hy')
plt.ylabel("amplitude at x=0")
plt.xlabel("time step")
plt.legend()

I'm not sure what is causing this sudden drop. Is there any renormalization process in MEEP when starting the simulation? I'd appreciate any feedback and suggestions. Thank you for your help and time!

[stevengj]

During timestepping, E is computed from D = εE, so you really need to initialize D, I think.

Probably we should compute D from E (or vice versa) automatically when you manually initialize the fields…

[m-h-hsieh]

Hi Professor Johnson,

Initializing D field works! Thank you so much for your help. I appreciate it.