H and B fields mismatch

[mochen4]

The H and B field returned by get_dft_array disagree even without magnetic material. Here is a simple test case without any material:

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

mp.verbosity(0)
resolution=20
cell_size = mp.Vector3(4,4)
pml_layers = [mp.PML(1.0)]


fcen = 1/1.55
width = 0.3
fwidth = width * fcen
source_size    = mp.Vector3(0,0,0)
src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)
source = [mp.Source(src,
                    component=mp.Ez,
                    center=mp.Vector3(),
                    size = source_size
                    )]

sim = mp.Simulation(cell_size=cell_size,
                    boundary_layers=pml_layers,
                    geometry=[],
                    sources=source,
                    resolution=resolution)

bhde = sim.add_dft_fields([mp.Bx, mp.Hx, mp.Dz, mp.Ez], [1/1.55], center=mp.Vector3(0.5,0), size=mp.Vector3(0, 0.2), yee_grid=True)

sim.run(until=500)

ez=sim.get_dft_array(bhde, mp.Ez, 0)
dz=sim.get_dft_array(bhde, mp.Dz, 0)
hx=sim.get_dft_array(bhde, mp.Hx, 0)
bx=sim.get_dft_array(bhde, mp.Bx, 0)
print("Ez:", ez)
print("Dz:", dz)
print("Hx:", hx)
print("Bx:", bx)
hnorm = np.sum(np.conj(hx)*bx) 
print("<Hx, Bx>", hnorm)
print("Re/Im", np.real(hnorm)/np.imag(hnorm), "resolution:",resolution)

And it prints

Ez: [-0.87498033+1.1439195j -0.8447642 +1.1783391j -0.8342971 +1.1897874j -0.84476435+1.178339j  -0.874980+1.1439193j]
Dz: [-0.87498033+1.1439195j -0.8447642 +1.1783391j -0.8342971 +1.1897874j -0.84476435+1.178339j  -0.874980+1.1439193j]
Hx: [ 0.2837132 -0.23033781j  0.16987896-0.14916535j  0.05652314-0.05164168 -0.05652316+0.05164165j -0.16987893+0.14916532j -0.28371346+0.23033786j]
Bx: [ 0.2950155 -0.21567228j  0.17721598-0.14036983j  0.05906626-0.04871255j -0.05906628+0.04871247 -0.17721586+0.14036942j -0.29501557+0.21567239j]
<Hx, Bx> (0.38055018+0.019299306j)
Re/Im 19.718334 resolution: 20

The E and D fields agree, as expected; but H and B fields don't. In fact, Re<H, B> / Im<H, B> is approximately the resolution.

[stevengj]

Maybe add some print statements to update_dfts to see if they are DFT-ing the same grid points etcetera?

meep/src/dft.cpp

Lines 259 to 268 in abea7ed

	void dft_chunk::update_dft(double time) {
	if (!fc->f[c][0]) return;
	const int Nomega = omega.size();
	for (int i = 0; i < Nomega; ++i)
	dft_phase[i] = polar(1.0, omega[i] * time) * scale;
	int numcmp = fc->f[c][1] ? 2 : 1;
	PLOOP_OVER_IVECS(fc->gv, is, ie, idx) {

[mochen4]

It seems that H fields and B fields are not updated at the same time. is_magnetic is False for B fields and True for H fields

meep/src/dft.cpp

Line 254 in abea7ed

cur->update_dft(is_magnetic(cur->c) ? timeH : timeE);

so H fields are updated at timeH but H fields are updated at timeE.

[stevengj]

Good catch! Looks like an easy fix.