x parity constraints for eigenmodes

[stevengj]

This would be nice, as discussed in NanoComp/mpb#72. However, implementing it in MPB is rather subtle as discussed in NanoComp/mpb#74 because of the interaction with MPB's transverse basis.

Seems like it might be more straightforward to just do a cyclic rotation of the coordinates in Meep's src/mpb.cpp so that the x direction is always the propagation direction.

[stevengj]

i.e. if the waveguide is in the y direction, rotate the coordinates to (y,z,x), and if the waveguide is in the z direction, rotate them to (z,x,y)

[stevengj]

(The tricky part is just to make sure that we do this consistently. You have to transform all the inputs, lattice vectors etcetera, and also all of the outputs like fields and group velocity.)

[stevengj]

Examples of place that would need to be change:

• meep/src/mpb.cpp

  Lines 52 to 63 in 5f2d84f

  	vec p(eps_data->dim == D3 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1], o[2] + r[2] * s[2])
  	: (eps_data->dim == D2 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1]) :
  	/* D1 */ vec(o[2] + r[2] * s[2])));
  	const fields *f = eps_data->f;
  	eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
  	eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));
  	eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega));
  	ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
  	ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Z, p, omega)), 0);
  	ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, omega)), 0);

• meep/src/mpb.cpp

  Lines 52 to 63 in 5f2d84f

  	vec p(eps_data->dim == D3 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1], o[2] + r[2] * s[2])
  	: (eps_data->dim == D2 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1]) :
  	/* D1 */ vec(o[2] + r[2] * s[2])));
  	const fields *f = eps_data->f;
  	eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
  	eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));
  	eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega));
  	ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
  	ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Z, p, omega)), 0);
  	ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, omega)), 0);

• meep/src/mpb.cpp

  Lines 299 to 328 in 5f2d84f

  	switch (gv.dim) {
  	case D3:
  	o[0] = eig_vol.in_direction_min(X);
  	o[1] = eig_vol.in_direction_min(Y);
  	o[2] = eig_vol.in_direction_min(Z);
  	s[0] = eig_vol.in_direction(X);
  	s[1] = eig_vol.in_direction(Y);
  	s[2] = eig_vol.in_direction(Z);
  	kcart[0] = kpoint.in_direction(X);
  	kcart[1] = kpoint.in_direction(Y);
  	kcart[2] = kpoint.in_direction(Z);
  	break;
  	case D2:
  	o[0] = eig_vol.in_direction_min(X);
  	o[1] = eig_vol.in_direction_min(Y);
  	s[0] = eig_vol.in_direction(X);
  	s[1] = eig_vol.in_direction(Y);
  	kcart[0] = kpoint.in_direction(X);
  	kcart[1] = kpoint.in_direction(Y);
  	kcart[2] = beta; // special_kz feature
  	empty_dim[2] = true;
  	break;
  	case D1:
  	o[2] = eig_vol.in_direction_min(Z);
  	s[2] = eig_vol.in_direction(Z);
  	kcart[2] = kpoint.in_direction(Z);
  	empty_dim[0] = empty_dim[1] = true;
  	break;
  	default: abort("unsupported dimensionality in add_eigenmode_source");
  	}

• meep/src/mpb.cpp

  Line 395 in 5f2d84f

  kmatch = G[d - X][d - X] * k[d - X]; // k[d] in cartesian

  • In particular, the d direction is in Meep coordinates, but kdir and G will be in MPB coordinates, so you will have to use a "rotated" version of d as an index into these arrays.

• meep/src/mpb.cpp

  Lines 410 to 413 in 5f2d84f

  	k[d - X] = kmatch * R[d - X][d - X]; // convert to reciprocal basis
  	if (eig_vol.in_direction(d) > 0 &&
  	fabs(k[d - X]) > 0.4) // ensure k is well inside the Brillouin zone
  	k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;

• meep/src/mpb.cpp

  Line 504 in 5f2d84f

  k[d - X] = kmatch * R[d - X][d - X];

• meep/src/mpb.cpp

  Line 645 in 5f2d84f

  return vec(edata->Gk[0], edata->Gk[1], edata->Gk[2]);

• meep/src/mpb.cpp

  Line 418 in 5f2d84f

  set_maxwell_data_parity(mdata, parity);

  (you'll have to translate the parity to the rotated coordinates)

Note that you have to cycle both the (x,y,z) coordinates and the field components (Ex,Ey,Ez). There will be some coordcycle which says to cycle the coordinates by 0, 1, or 2. Going from Meep to MPB you would cycle by +coordcycle and going from MPB to Meep you would cycle in the opposite direction (by -coordcycle)

[stevengj]

For example, d - X becomes (d - X + coordcycle) % 3.

[stevengj]

Note that the parity flags are defined here. In particular, you should think of parity as a bit field of flags: each bit of the integer is set by one of the *_PARITY flags.

We define our own bits for EVEN_X_PARITY = 16 (4th bit) and ODD_X_PARITY = 32 (the 5th bit). The caller can pass any combination of these. Your code will then have to shift these flags (e.g. if X shifts to Y then you have to shift the 4th and 5th bits to the 2nd and 3rd bits).

[oskooi]

Two questions related to implementation:

1. How should the direction of the "waveguide" (which can only be one of the three coordinate directions x/y/z) be determined from the user data? One approach is to set it to the direction normal to eig_vol (assuming it is a line in 1d or a plane in 2d). However, the situation is more complicated if eig_vol has all of its dimensions with non-zero size (e.g., a photonic-crystal waveguide). Assume the user does not set direction=mp.NO_DIRECTION and specifies the waveguide direction using kpoint_func.
2. Should coordcycle be a function (with a single input parameter eig_vol) or a global variable? How can the existing interface of meep_mpb_eps and fields::get_eigenmode (the two functions in src/mpb.cpp where coordcycle is to be used) be preserved?

[stevengj]

1. see how it is already determining the k-point direction.
2. coordcycle will have to be stored in the eigenmode_data