How Does MEEP Represent CustomSource?

[kombatEldridge]

Hello, I am working on creating a dipole emitter in MEEP, but would just like to point to some lines of code in the project that correlate the CustomSource function in python (and its src_func) to an electric field in the simulation. In my project, I have an array of how a dipole moment changes over time, and would like to produce an electric field from its changes using CustomSource. Can anyone help me? Thank you!

[stevengj]

With a source, you are (by default) specifying the current density J in Maxwell's equations. There is no simple relationship (in general) between the current source and the resulting electric field, because it depends on the surrounding materials.

This is a Meep FAQ.

If you have a dipole moment (polarization) density P, the current density J=dP/dt is the time derivative of this. However, if you set the flag is_integrated=True in CustomSource, then your custom source is precisely the dipole moment density P.

[kombatEldridge]

Thank you for your help!

[kombatEldridge]

Hello, circling back to this Q&A @stevengj. Can you help me understand the physics of flipping is_integrated to True? At a basic level, my code places the dipole moment into an dict where a CustomSource can look at for its amplitude. Each moment is the value to a key which is named by the time step (in Meep Units). When the CustomSource is called, it looks for the value in the dict corresponding to the current time.

mp.Source(
    mp.CustomSource(src_func=self.chirpx, is_integrated=True),
    center=self.positionMolecule,
    component=mp.Ex
)
...

def chirpx(self, t):
    return self.dipoleResponse['x'].get(str(round(t, self.decimalPlaces)), 0)

When printing the time for chirpx, here is an example:

DEBUG: CustomSource being called at 0.000375.
DEBUG: CustomSource being called at 0.000625.
DEBUG: CustomSource being called at 0.000375.
DEBUG: CustomSource being called at 0.000500.
DEBUG: CustomSource being called at 0.000750.
DEBUG: CustomSource being called at 0.000500.
INFO: Getting Electric Field at the molecule at time 0.000500.
DEBUG: CustomSource being called at 0.000625.
DEBUG: CustomSource being called at 0.000875.
DEBUG: CustomSource being called at 0.000625.
DEBUG: CustomSource being called at 0.000750.
DEBUG: CustomSource being called at 0.001000.
DEBUG: CustomSource being called at 0.000750.
INFO: Getting Electric Field at the molecule at time 0.000750.

However, if is_integrated is set to False, the output looks like I would expect (where CustomSource is called every dt/2):

DEBUG: CustomSource being called at 0.000375.
DEBUG: CustomSource being called at 0.000375.
DEBUG: CustomSource being called at 0.000500.
DEBUG: CustomSource being called at 0.000500.
INFO: Getting Electric Field at the molecule at time 0.000500.
DEBUG: CustomSource being called at 0.000625.
DEBUG: CustomSource being called at 0.000625.
DEBUG: CustomSource being called at 0.000750.
DEBUG: CustomSource being called at 0.000750.
INFO: Getting Electric Field at the molecule at time 0.000750.

Can you explain why the is_integrated=True CustomSource is being called at these odd times?

(You can see more about my project here in case you need more context.)

[stevengj]

The first thing to realize is that, at every half-timestep, the update() method of the current source is called, which for the current time $t$ computes both the current $J$ and (possibly) the dipole moment $P$. The two are related by $J(t) = \partial P/\partial t$, and caches the values in the source object. These cached values are then re-used whenever the current or the dipole moment is needed.

When you specify a Python function as the time dependence, it gets wrapped in a custom_src_time object in C++. There are two cases:

• If is_integrated=False, your custom function is interpreted as $J(t)$. (It is also used as $P(t)$ in the update function, but this value ends up not being used. That's why your function is called twice per half-timestep.) During the subsequent field update, this current $J$ is added into the displacement field $D$ update from Ampere's law $\partial D/\partial t= \nabla \times H - J$.

• If is_integrated=True, your custom function is interpreted as the dipole moment $P(t)$ — this is one evaluation per half-timestep. (The update function also computes $J = dP/dt$ a finite-difference approximation $J ≈ (P(t+\Delta t)-P(t))/\Delta t$, but it ends up not using this value. That's why your custom source is called at two extra times you don't expect, one of which is offset by Δt! It's harmless — the value isn't used.) In the subsequent field timestep, the dipole moment P directly when computing $E = D - P$ (your source term is included in the P).

In short, the times at which your custom function are evaluated include some extra evaluations that Meep ends up discarding. We don't worry about this because the time dependence is only evaluated once per (half) timestep, so its cost is normally negligible compared to the cost of timestepping the fields. But I realize that it could be confusing if you are inspecting the timeseries as you are doing above.