Modal group index calculation without finite differences #2010
Comments
|
Yeah this would be great to have. By the way one of the reasons for keeping our frontend open (including the local mode solver!) is hoping the community can also contribute with things like these that we may not find time for right now. It sounds like you have quite a lot of expertise on the matter so if you want to give it an initial try, we'd certainly be happy to support as much as possible. |
|
I've poked around at what it would take -- still pretty new to the tidy3D interface so I'm hesitant to write something crappy and submit it 😅 There are some other things I'm going to try my hand at first, in particular dielectric smoothing for generally-anisotropic materials. That'd be a much bigger improvement for the projects I'm working on in the forseeable future (and would make the solution to this problem better!). Might rope @doddgray into it. |
|
I would be glad to help with this. Calculating group index from a single frequency solution should be faster, more accurate and straightforward. It's the ratio of the integrated energy density to the integrated poynting flux dotted into the plane. If you want it to be accurate you need to add the dispersive terms to the energy density. I tried doing this by hand with tidy3D mode data a year or so ago and got confused about whether the field and dielectric tensor elements corresponded to the same spatial grids or whether they were on the Yee lattice. Maybe this is documented now. What is the best/expected way to calculate Poynting flux and (ED + HB) with Tidy3D mode data? That's all there is to it. |
|
Thanks! So the fields in Tidy3D can be either colocated to grid boundaries or returned on the Yee grid, depending on the |
|
Thanks, that sounds similar to what I tried previously so good to know I wasn't too far off. IIRC I ended up with something that seemed nearly correct but I was worried that I was performing unnecessary or slightly incorrect field & permittivity interpolations (mu also 1/I for me). If calculated consistently with the field data grid locations I think this should result in modal group index values that are accurate to smth comparable to the eigenvalue tolerance used for mode solving. I think i was doing some inconsistent interpolations to colocate the fields and this was noticeably degrading that accuracy. I'll dig up the code today and share here. I bet you could fix it in seconds. This reminds me- Do field and tensor elements always/usually live on the Yee lattice in frequency domain solvers (like the mode solver)?. I thought the Yee lattice was needed for stable discretized time domain integration and would have naively thought to use colocalized fields and tensors for the mode solver. |
|
Oh right IIRC out-of-plane field elements are needed for accurate (ED+HB) and i was confused how to do this for a mode solver defined in a 2D plane when the fields are on the Yee grid. |
The Yee grid is actually more about accuracy than about stability, although if a different grid is used in the time domain, a stability analysis would certainly need to be done. But yeah it's just a way to directly do centered derivatives on a grid that feels like twice higher resolution than a collocated one, with the same amount of work and memory cost. For this reason it is ubiquitously used in frequency-domain methods too, like FDFD. (Not sure about specific references, I'm sure there are many, e.g. https://www.youtube.com/watch?v=hv5lIx4u8mY&ab_channel=EMPossible) The mode solver is something like an FDFD solver but with a slightly different matrix and obviously doing an eigen solve rather than a linear solve. You can for sure write a mode solver on a collocated grid too. But apart from the convenience / potentially higher accuracy of the Yee grid, when the mode solver is strongly coupled to FDTD simulations as is typically the case, it's even more crucial for accuracy reasons to solve on the same grid that the FDTD uses. This is both to avoid interpolation errors if you need to interpolate the mode fields to the FDTD grid, and in order to use the same subpixel averaging in both solvers, which again is grid-dependent.
So along the normal direction the fields are actually always collocated to the mode plane. I'm not sure I understand this point fully but if you do need them anywhere out of that plane, you just need to apply the |
|
Thanks for those explanations. Here's what I tried before, which seemed at least slightly wrong at the time. Any corrections or suggested more efficient/accurate ways to do this would be greatly appreciated. # .... create a ModeSolver `mode_solver` ...
# ... solve for modes, generating `mode_data` ....
# ... choose mode `mode_index` to analyze ...
### Modal Energy Density `U` & Integral `W` ###
grid = mode_solver.grid_snapped # or `mode_data.grid_expanded`? or one of these grids' `.centers`?
Ex = mode_data.Ex.sel(mode_index=mode_index)
Ey = mode_data.Ey.sel(mode_index=mode_index)
Ez = mode_data.Ez.sel(mode_index=mode_index)
eps_xx = mode_solver.simulation.epsilon_on_grid(grid=grid, coord_key="Ex").interp_like(Ex)
eps_yy = mode_solver.simulation.epsilon_on_grid(grid=grid, coord_key="Ey").interp_like(Ey)
eps_zz = mode_solver.simulation.epsilon_on_grid(grid=grid, coord_key="Ez").interp_like(Ez)
Dx = eps_xx * Ex
Dy = eps_yy * Ey
Dz = eps_zz * Ez
U = (Ex.conj() * Dx + Ey.conj() * Dy + Ez.conj() * Dz).real # Modal Energy density [J/cm^3]
W = float(U.integrate(coord=["x", "z"]).real) # Modal Energy density per unit length, dU/dy [J/cm]
print(f"Modal Energy Density Per Length `W`:\t{W}")
### Poynting Flux `S` and Integral `P` ###
norm_factor = 1 # not used right now, needed for unitful flux and power quantities
Ex = mode_data.Ex.sel(mode_index=mode_index).squeeze(drop=True) * norm_factor
Ez = mode_data.Ez.interp_like(Ex).sel(mode_index=mode_index).squeeze(drop=True) * norm_factor
Hx = mode_data.Hx.interp_like(Ex).sel(mode_index=mode_index).squeeze(drop=True) * norm_factor
Hz = mode_data.Hz.interp_like(Ex).sel(mode_index=mode_index).squeeze(drop=True) * norm_factor
S = ( (Ez.conj() * Hx - Ex.conj() * Hz) - (Hz.conj() * Ex - Hx.conj() * Ez) ) / 4 # Modal Poynting Flux Density [W/cm^2]
P = float(S.integrate(coord=["x", "z"]).real) # Modal Power [W]
print(f"Modal Power Integral `P`:\t{P}")
### Modal Group Index ###
n_g = W / P # possibly missing factor of td.C_0, if so this is inverse modal group velocity
print(f"Modal Group Index `n_g`:\t{n_g}") |
|
Thanks! So this was giving NaNs at first, I believe something to do with mismatch of various grids on which things are computed. However I fixed it in a way that seems to produce meaningful results like this: The normalization seem off as the values are of the order of 1600-1800 for the three modes I tried for the mode solver in our notebook example. The difference doesn't seem to be (only) a td.C_0 factor. However, the ratio of the values is pretty close to the one computed using the finite difference derivative: Would be great if we can figure the final touches out! That said it's clear this will involve some effort to incorporate natively, to handle and test all edge cases, dispersive materials, etc. But figuring out this simple example should be a good start! |
|
Nice those numbers look reasonably close! Thank you for looking at this. Using the native modal Maybe worth looking at this analytic vs. finite difference group index comparison for a more confined mode ( In terms of final touches to make the analytic group index calculation generally correct and accurate, a few thoughts:
I assume the dielectric function is not re-calculated and smoothed at different frequencies used in the finite difference group index calculation currently used? That would be necessary to check the group index calculation including material dispersion. This may sound insignificant but it is very important for those of us who care about modal GVDs and dispersion engineering for pulses and nonlinear optics. Without material dispersion the modal group index is strictly wrong, typically at the few-percent level in cases I've focused on and the GVD is very wrong. |
The code is here. There are some intricacies to make sure that the flux computed this way matched exactly what the remote solver would also return (e.g. for a perfectly injected mode propagating in the FDTD grid), so the finite step size in the z direction is actually incorporated. This is the only point that I am not sure about whether it would require a re-computation of the flux, or whether it should be incorporated in the
Sorry, to clarify the normalization of the group index computed by the code above seems off, so what I am showing in the end of my message is the ratios of group indexes between the second and the first mode, and between the third and the second mode. The actual values are higher than 2 and the modes are well confined. So there is definitely still the issue of getting the normalization right to do the exact comparison.
Basically, these live at different locations, so any choice would require interpolation. Taking the grid boundaries is what we normally do for integrations. Some time ago we were actually using grid centers but switched to boundaries for some reasons that I don't need to get into, but the integration results would generally change very little. Regarding off-diagonal permittivity, do you mean this is needed when the underlying materials are fully anisotropic? Or are there other cases?
Yeah so like I mentioned above re flux, we may need to look into that, but I would leave it for last.
Yeah probably the averaging only works for lossless modes.
I see, yeah, we'd have to finite difference it, but it will be cheaper. But again, for starters we could just set this up to only work for non-dispersive cases and then expand later.
We run completely independent mode solves so independent structure parsing and permittivity smoothening is applied at every separate freuquency, including the small offset freqs used for the numerical derivatives. |
|
Thanks that all makes sense. I will check/fix the normalization today so the fraction of integrals is actually the group index. I know I figured it out in the past. I'm using tidy3D a lot so definitely excited to get this working. What if anything would be helpful for us to do? Do you want a user submitted PR or should we leave it to the pros? I have some more responses to yours, particularly about fully anisotropic tensor data but will type those up when I get to a computer. |
That would be great!
You can certainly attempt to PR it yourself if you wish, we're always happy to get external contributions. If left up to us, it may take some time to get to it. We are trying to hire a few more people though and this could be a good first project if we get some of them soon. |
|
I think I finally figured this out, or at least this group index calculation based on individual mode solutions seems much closer to correct than my earlier attempt. The main issues with my previous snippet were that I forgot the vacuum permittivity In this version I've fixed those and included both the Beyond that, I think I still need to
import tidy3d as td
import xarray as xr
from tidy3d.constants import ETA_0, EPSILON_0, MU_0, C_0
def calc_ng(
mode_data,
mode_solver,
grid = None,
check_density_balance:bool = False,
squeeze:bool = True,
):
""" Calculate the modal group index as
n_g = td.C_0 * W / P
where `W` is the modal energy density per unit length in [J/μm]
and `P` is the modal power in [W]
`P` = `mode_data.flux` and is typically normalized to 1W."""
if grid is None:
grid = mode_data.grid_expanded # mode_solver.grid_expanded
interp_kwargs = dict(
y=grid.centers.y,
z=grid.centers.z,
kwargs={"fill_value": "extrapolate"},
)
# First compute everything on the Yee grid
Ex = mode_data.Ex.interp(**interp_kwargs)
Ey = mode_data.Ey.interp(**interp_kwargs)
Ez = mode_data.Ez.interp(**interp_kwargs)
Hx = mode_data.Hx.interp(**interp_kwargs)
Hy = mode_data.Hy.interp(**interp_kwargs)
Hz = mode_data.Hz.interp(**interp_kwargs)
# Relative Dielectric Tensor Elements
eps_xx = xr.concat(
[mode_solver.simulation.epsilon_on_grid(grid=grid, coord_key="Ex", freq=freq).expand_dims(dim={'f':[freq,]}) for freq in mode_solver.freqs],
dim='f',
).interp(**interp_kwargs)
eps_yy = xr.concat(
[mode_solver.simulation.epsilon_on_grid(grid=grid, coord_key="Ey", freq=freq).expand_dims(dim={'f':[freq,]}) for freq in mode_solver.freqs],
dim='f',
).interp(**interp_kwargs)
eps_zz = xr.concat(
[mode_solver.simulation.epsilon_on_grid(grid=grid, coord_key="Ez", freq=freq).expand_dims(dim={'f':[freq,]}) for freq in mode_solver.freqs],
dim='f',
).interp(**interp_kwargs)
# Modal E-field Energy density [J/μm³] scaled by C_0 ( C_0 * EPSILON_0 * n² * |E|² = (1/Z_0) * n² * |E|² )
U_E = (
(eps_xx * Ex.abs**2).interp(**interp_kwargs)
+ (eps_yy * Ey.abs**2 ).interp(**interp_kwargs)
+ (eps_zz * Ez.abs**2 ).interp(**interp_kwargs)
) / td.ETA_0
# Modal H-field Energy density [J/μm³] scaled by C_0 ( C_0 * MU_0 * |H|² = Z_0 * |H|² )
U_H = (
( Hx.abs**2).interp(**interp_kwargs)
+ ( Hy.abs**2 ).interp(**interp_kwargs)
+ ( Hz.abs**2 ).interp(**interp_kwargs)
) * td.ETA_0
if check_density_balance:
W_E = U_E.integrate(coord=["y", "z"]).real # Modal E-field Energy density per unit length, dU/dy [J/μm]
print(f"\tW_E:\t{W_E}")
W_H = U_H.integrate(coord=["y", "z"]).real # Modal H-field Energy density per unit length, dU/dy [J/μm]
print(f"\tW_H:\t{W_H}")
W = (W_E + W_H) / 4 # TODO: Where should this factor of 4 come in?
else:
W = ( U_E + U_H ).integrate(coord=["y", "z"]).real / 4 # Modal Energy density per unit length, dU/dy [J/μm]
if squeeze:
W = W.squeeze(drop=True) # Drop singleton spatial dimension coordinate
# n_g = W / P * C_0, but we have already scaled W by C_0 and P = 1 [W] for Tidy3D's normalized mode fields
# so here, n_g = W. So we just return W
return W |
Problem
Modal group index calculations take a while, and cost a lot. That's because they use finite differences.
Aside from the cost, it becomes quite problematic to finite difference when simulating GVD. Not only are you propagating even more numerical error (and are susceptible to choice of step size, though this is a lesser issue), frequently photonics designers are searching for places where GVD is flat or zero -- exactly where numerical noise will be the most visible!
Solution I'd like
It's possible to analytically calculate the modal group index based off the effective index, modal field distribution, and permittivity distribution. This was referenced in a comment on #1169, but seems to have been lost in the ether when that issue closed.
References for equations to implement
The referenced paper here has a rigorous treatment of this (Eqs. 14 and 15). There are additional simplifications which can be made, particularly considering the case of Z-invariant (X-invariant, whatever) propagation, but this is probably easy enough to implement as-is.
Po-Ru and Steven's derivation does not include the effects of anisotropy. It's not a difficult extension -- I've derived it in the past for the case of negligible loss (maintaining a hermetian operator) with fully-anisotropic and dispersive materials. If there's interest in pushing this forward I can dig up those derivations. A super brief version exists in this paper, Section 6 of the supplement, leaving it in a general form and then the solution in the plane-wave basis.
A problem with this would be that material dispersion can strongly affect the EM energy density and must be accounted for. In general you already have dispersive materials so it's juts a question of dielectric smoothing. While dielectric smoothing is underway, it'd be useful to get these calculations going with the restriction that it only works for un-smoothed models (much like fully anisotropic materials are at the moment).
The text was updated successfully, but these errors were encountered: