Hi Steven, Hello meep users, I have been interested in the details of the timing of electric and magnetic sources in meep, since if one puts an electric and magnetic source at the same location, it is possible to make a source that emits only to one side. This can be useful sometimes.
However, it seems to me that the magnetic sources are evaluated one time step ahead of the electric ones: they are out of sync. If you give them the same src-time object, then the fields for the electric source come out later. Usually this small delay does not cause major problems, of course, since most simulations do not depend of the precise timing of the source, and one time step is supposed to be a small delay anyway. But if one tries to combine electric and magnetic sources, then the difference is important, and in any case, an extra delay is just not proper. To demonstrate the issue, I have tried to come up with a very simple example: looking directly at the electric and magnetic fields emitted by a pulsed source in empty space, in one dimension. In this case the solution to Maxwell’s equations is of course well known: for a current source j(t) in the origin, the electric field is E(z,t)= - (1/2) j(t-|z|), and the magnetic field is H(z,t)=E(z,t) sign(z). This simple analytical result facilitates a direct comparison with meep’s calculated fields. For definiteness, I have attached below a .ctl file that I used to obtain the meep fields, as well as a small matlab script to plot the results. Looking at the fields, it immediately becomes apparent that the sources are not currents j(t) (as the notation in the meep manual suggests), but rather polarizations p(t) (as Steven has pointed out in the discussion list sometime), with j(t)=d/dt p(t). The electric and magnetic fields from meep are exactly (anti)symmetric in z, for both types of source (E and H), with the signs as expected, and the magnitudes are in excellent agreement as well (but the timing is more delicate). The (anti)symmetry in z confirms that my understanding of the z coordinates is ok (the script obtains the z positions from the hdf file, and the time of the field snapshot from the file name). Looking at the timing of the fields, it appears that the H field lags the E field by half a time step (perhaps not so surprising for an FDTD calculation). A closer look at the source code in “dft.cpp” confirms the time ordering: the update routine “fields_chunk::update_dfts(double timeE, double timeH)” gets called in “fields::update_dfts()” as “update_dfts(time(), time() - 0.5 * dt);”, suggesting that the magnetic field is indeed lagging half a time step. For an H source, it appears that meep’s E field nicely matches the analytical expression, based on the time in the file name, and the H field nicely matches the analytical expression if we take the time half a time step earlier, as is normal for the Yee algorithm. However, for an E source, it appears that both fields are lagging one extra time step. The matlab script plots the E and H fields, with extra signs automatically added so all the amplitudes have the same sign, regardless of the type of source or the sign of z, so that all the curves overlap. The result is shown in “ehexample1.png”. A zoom in on the peak , in “ehexample2.png”, clearly shows the half time step delay for the H field. A comparison of the peak for an E vs. an H source (“ehexample2.png” vs. “ehexample3.png”, same axis scaling) clearly shows the extra time step delay for the E source (compare the peak positions, or the top axes of the two graphs). (The extra delay is taken into account in the analytical result, so the solid and dashed curves keep overlapping in the plot.) Taking a closer look at the source code, I think I am able to point out where the discrepancy arises. Here is an excerpt from ./src/step.cpp: calc_sources(time() - 0.5 * dt); // for H sources step_h(); if (!disable_sources) step_h_source(); calc_sources(time()); // for E sources step_d(); update_e_from_d(); update_from_e(); t += 1; update_dfts(); In step_h_source(), the field is updated using current*dt, where current(t)=((dipole(t+dt) - dipole(t)) / dt) (according to meep.hpp). Since t is time()-0.5*dt from calc_sources() above, the update current is based on dipole(time()+0.5*dt)-dipole(time()-0.5*dt), which is appropriate to update H from H(time()-0.5*dt) to H(time()+0.5*dt), and the H sources are fine. In contrast, in update_e_from_d.cpp, d_minus_p is updated using dipole(t), where t is still time(), from the second call to calc_sources() above. The update happens *after* step_h() and step_d(), i.e., the fields are already propagated to time()+dt. So presumably the update_e_from_d() should use sources evaluated at time()+dt, instead of at time(). I think this is the reason that the E and H fields for an E source are delayed one extra time step. What would be the best way to fix this? Is just replacing calc_sources(time()); // for E sources by calc_sources(time()+dt); // for E sources ok? Best regards, Mischa (See attached file: ehexample.ctl)(See attached file: ehexample.m)(See attached file: ehexample1.png)(See attached file: ehexample2.png)(See attached file: ehexample3.png)
ehexample.ctl
Description: Binary data
ehexample.m
Description: Binary data
<<attachment: ehexample1.png>>
<<attachment: ehexample2.png>>
<<attachment: ehexample3.png>>
_______________________________________________ meep-discuss mailing list [email protected] http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
