Nonuniform Mesh API

[tylerflex]

Let's use this to iterate on the nonuniform mesh architecture and API?

General picture

Introducing the notion of a MeshSpec that contains specifications about how the user wants the space to be discretized (minimum grids per wavelength, etc.).

a MeshSpec + Simulation + other parameters ultimately determines a set of coordinates in 1D that define the yee cell boundaries.

If these mesh specs will be contained directly in the Simulation, when Simulation.grid property is called, the coordinates from the mesh specs should be evaluated loaded into a Grid().

If these mesh specs are defined externally to set grid_size, then the code functions as it does now.

Main API Options

1. Add MeshSpec() as another accepted value in one of Simulation.grid_size. When computing grid boundaries, add a special case to handle a MeshSpec().
2. Keep Simulation.grid_size as float or array and make an external plugin to generate a grid_size given MeshSpec(), simulation and other parameters.3.
3. Make MeshSpec() the only way to define discretization, but supply other MeshSpec() types to handle uniform or user specified grid sizes.4.
4. An Interim solution of 2 and then transition to 3 when it is more stable.

Considerations

• a. If no sources, need to define some notion of a wavelength or length scale. Can define this in the MeshSpec directly, or add an optional freq Field to the Simulation, or have it be a kwarg in the plugin init.
• b. What gets written into the json? A MeshSpec would be easier to read than raw coords.
• c. How extendable is the architecture?

Thoughts

• We should avoid code that looks like: if meshSpec1: do this elif meshSpec2 do this. Ideally, the MeshSpec should contain it's own way to generate a 1D grid from a set of parameters, whether that be size & center, structures, simulation, etc.

[momchil-flex]

My opinion is actually do 1 and then transition to 3. Like you say, I believe that we should be writing MeshSpec in the json, not the raw coords. This is not just because it's easier to read. It's also because you're losing potentially valuable information about how the mesh was generated, in case you want to modify that.

What I'm thinking right now is:

Simulation(grid_size, ....) # grid_size can be float, list of floats, or MeshSpec

    @property
    def grid(self):
        mesher = Mesher(self, grid_size)
        return mesher.grid

The Mesher can also be available as a plugin. The point however is that we don't need to use that plugin in order to initialize a simulation. We just need to provide the correct grid_size (which can be a MeshSpec) that we want. This is similar to how you don't need to use the mode solver plugin to define the mode sources and monitors in the simulation. You provide the ModeSpec, not the solved fields. The plugin is just a convenience to examine the modes and their fields before settling on a ModeSpec. Similarly, you provide a MeshSpec for the nonuniform mesh, not the grid or grid boundaries that are computed by the Mesher for a given MeshSpec. But, you can use the Mesher in advance to examine how things look before settling on a MeshSpec.

[tylerflex]

I dont think the Mesher can be a plugin since Simulation depends on it. There's also a circular import issue that sets off some red flags for me:

mesher.py

from .simulation import Simulation
from .types import GridSize
class Mesher:
    sim : Simulation
    grid_size : GridSize

simulation.py

from .mesher import Mesher
class Simulation

    @property
    def grid(self):
        mesher = Mesher(self, grid_size)
        return mesher.grid

It's kind of inherently circular.

One workaround is similar to something shashwat suggested, which is to pass the relevant simulation parameters. In fact a while ago when we considered grid spec, I had something like this, but it was very clunky.

I'm not sure if there's another approach if we go with 1 or 3. We'll need something that takes a Simulation and MeshSpec(s) and spits out a grid. not sure where that code should lie. I'll continue to stew it over

[tylerflex]

I guess the way we're currently doing it is essentially #1 (without mesher) but I was never 100% happy with it.
https://github.com/flexcompute/tidy3d/blob/develop/tidy3d/components/simulation.py#L1302

[momchil-flex]

Yeah you're right. As long as we want to give MeshSpec-s to the simulation (which I do), it seems to me that the mesher is an integral part of the simulation rather than a plugin, because the simulation needs it to make its grid, which is needed in many places.

Maybe I've got it - it just requires a small API change. Instead of Simulation(grid_size=(...), ...), what about Simulation(mesher=td.Mesher(dx=0.1, dy=(0.1, 0.2, 0.3), dz=td.MeshSpec(min_steps_per_wvl=10))? Now we can just use that Mesher object, which is internal to the Simulation, in a very similar way that we thought of using the plugin:

mesh_spec = td.MeshSpec(min_steps_per_wvl=12)
sim = td.Simulation(mesher=td.Mesher(dx=mesh_spec, dy=mesh_spec, dz=mesh_spec))
sim.plot(z=0)
sim.plot_grid(z=0)

# add a structure to a list of refinement structures in Mesher
sim.mesher.refine(box=td.Box(...), medium=td.Medium(...)) 

sim.plot(z=0)
sim.plot_grid(z=0)

Internally, the Mesher just has something like a get_grid(structures, freq0) method that is called in Simulation.grid using self.structures and some sort of simulation central frequency.

It's a separate question if we want to combine the grid plotting in sim.plot.

[tylerflex]

Ok so here the Mesher is kind of like a "global" grid spec and contains the methods needed to construct the grid? Those methods would still take in things like self.size, self.center, self.structures passed from simulation (self)?

I think this is getting better. I think that since the Mesher is ultimately getting turned into specs, it might make more sense for it to just be called MeshSpec and then the individual dx, dy, etc. could be called something like MeshSpec1D, since it's still just a specification with the discretization methods defined internally but not stored.

I'm still leaning towards an interim plugin for generating the grid_size in the same format that we use now, just to be conservative. But maybe it is worth just redoing everything to use mesh spec now, as you suggest. It's just a major change to introduce that will break backwards compatibility. Are there any other risks that you can anticipate?

[momchil-flex]

Ok so here the Mesher is kind of like a "global" grid spec and contains the methods needed to construct the grid? Those methods would still take in things like self.size, self.center, self.structures passed from simulation (self)?

Yeah.

I think this is getting better. I think that since the Mesher is ultimately getting turned into specs, it might make more sense for it to just be called MeshSpec and then the individual dx, dy, etc. could be called something like MeshSpec1D, since it's still just a specification with the discretization methods defined internally but not stored.

Agree on calling it MeshSpec. How would the MeshSpec1D work to allow init via float, list of floats, or nonuniform parameters? Classmethods?

I'm still leaning towards an interim plugin for generating the grid_size in the same format that we use now, just to be conservative. But maybe it is worth just redoing everything to use mesh spec now, as you suggest. It's just a major change to introduce that will break backwards compatibility. Are there any other risks that you can anticipate?

It is a major change so painful yeah, but as usual better earlier than later? Like, the question is if that's what we want to have eventually, and if yes - then we should do it now. I kinda feel like we'll have to change the grid_size part of the API eventually. Technically, we could patch it such that grid_size, if passed, sets the MeshSpec field, and vice versa...

I don't anticipate any other risks apart from just the need to change like every simulation init in all our docs and tests.

[tylerflex]

How would the MeshSpec1D work to allow init via float, list of floats, or nonuniform parameters? Classmethods?

I'd say we just have a few options for MeshSpec1D subclasses, like AutoMeshSpec(min_wvl_mat) or UniformMeshSpec(dl=0.2) or CustomMeshSpec(boundaries=(-1, -.5, 0, .5, 1)). But the base class has a method that takes simulation parameters and the Fields stored in the mesh spec (self) and returns some 1D coords used to construct the grid in the main MeshSpec.

It is a major change so painful yeah, but as usual better earlier than later? Like, the question is if that's what we want to have eventually, and if yes - then we should do it now. I kinda feel like we'll have to change the grid_size part of the API eventually. Technically, we could patch it such that grid_size, if passed, sets the MeshSpec field, and vice versa...

Yea.. just throwing this idea out there (similar to what you said): we add mesh_spec into simulation fields, which can be set however we want to do it eventually. We still support user setting grid_size as it is now, throwing a warning that this will be deprecated, and then the validator constructs the mesh spec corresponding to the grid_size field. For example see the CustomMeshSpec and UniformMeshSpec above. Whether grid_size overrides mesh_spec or the other way around we can discuss further.

[momchil-flex]

AutoMeshSpec(dl0: float, max_scale: float = 1.4, lambda0: float = None)
Error if lambda0 is None and no sources in Simulation and there are disperssive materials

AutoMeshSpec.from_min_steps(min_steps_per_wvl: float, max_scale: float = 1.4, lambda0: float = None)

[weiliangjin2021]

In the branch weiliang/nonuniform, now one can setup the mesh in this way:

mesh_x = td.AutoMeshSpec(wavelength = 1, min_steps_per_wvl = 20) 
mesh_y = td.UniformMeshSpec(dl = 0.1)
mesh_y = td.CustomMeshSpec(dl = [0.1, 0.1, 0.1]) 
mesh_spec = td.MeshSpec(
    mesh_x = mesh_x,
    mesh_y = mesh_y,
    mesh_z = mesh_z,
)

sim = td.Simulation(
    ...
    mesh_spec = mesh_spec,
    ...
)

Some issues and TODOs:

• The wavelength is supplied to AutoMeshSpec along each direction; ideally it should be supplied to MeshSpec. So maybe just a simple change here;
• The validator _warn_grid_size_too_small in simulation.py is yet to be updated.
• MeshSpec.refine is yet to be implemented. Do we need a private variable for this purpose, so that a structure of higher index will be appended to the structure list for generating finer mesh? This structure will only be visible in MeshSpec.
• I prefer to setup AutoMeshSpec from dl0 through classmethod, since setting up from min_steps_per_wvl is more natural, although it relies more on wavelength. The former is to be implemented.

[momchil-flex]

1. I agree, I think we can validate that if an AutoMeshSpec is provided along any of x, y, z, then wavelength is also provided:

mesh_spec = td.MeshSpec(
    wavelength = 1,
    mesh_x = mesh_x,
    mesh_y = mesh_y,
    mesh_z = mesh_z,
)

2. For auto nonuniform, we can just set min_steps_per_wvl to be at least 6 (or 8). For dl0 it will be an easy check. For custom nonuniform, we can try to include it, but it's such an advanced feature that I don't even count it as 100% necessary.

3. We can postpone this for next patch.

4. I agree.

[momchil-flex]

I've created momchil/nonuniform. I've tried to push somewhat atomic commits so that you can follow what I'm doing if you want to later on. So far I have rebased vs. develop, made all tests run apart from the "IO" test which needs some input from Tyler, and finished 1. from the list above, i.e. I put wavelength in MeshSpec. It is not fully satisfying right now because if it's missing and there are no sources, the error will only be raised upon a call to Simulation.grid and not on Simulation init. But I am not sure if it can be safely moved as a validator without anything breaking if the user manually resets sources or mesh_spec.wavelength.

I think we can postpone 3. and 4. from the list above for next patch. What still needs to be done is:

• See if wavelength/sources can be validated upon simulation init
• Update _warn_grid_size_too_small
• Fix the reading from file. Right now I get an error because it cannot parse the correct MeshSpec1D.
• Update all tests on both backend and frontend to use mesh_spec instead of grid_size
• Update all example notebooks to use mesh_spec instead of grid_size. Examine the meshes to make sure they make sense.

[momchil-flex]

Here's another update.

• Move wavelength from AutoMesh to MeshSpec.
• Validating wavelength upon Simulation init.
• Update _warn_grid_size_too_small for the case of uniform mesh.
• Update almost all frontend tests to use mesh_spec instead of grid_size. Left a few to test that grid_size also still works.
• Fix the reading from file. Right now I get an error because it cannot parse the correct MeshSpec1D.
• Update most of the backend tests to use mesh_spec.
• Update all example notebooks to use mesh_spec instead of grid_size. Examine the meshes to make sure they make sense. Make sure that empty simulation normalization, where needed, works well.
• Update the What's New notebook showing how to use the new mesh_spec.

Also, I introduced MeshSpec classmethods .auto and .uniform. They take the same arguments as AutoMesh and UniformMesh, and apply the same mesh strategy on all three directions. Now, there are multiple ways to conveniently define the mesh. For example:

# default min_steps_per_wvl = 10, wavelength defined explicitly or through sources
Simulation(mesh_spec=MeshSpec(wavelength=1.5)) 
Simulation(mesh_spec=MeshSpec(), sources=[source])

# apply a custom min_steps_per_wvl along all three directions
Simulation(mesh_spec=MeshSpec.auto(wavelength=1.5, min_steps_per_wvl=16)) 

# Uniform mesh in all directions if you don't want to define wavelength yet
Simulation(mesh_spec = MeshSpec.uniform(dl=0.1))

After giving it some thought, I am not that enthusiastic about adding a from_dl0 classmethod for the MeshSpec (or the individual AutoMesh). This is because it is not at all clear what the wavelength and min_steps_per_wvl fields should be set to (we can't just do wavelength = 1, min_steps_per_wvl = 1 / dl0, because the wavelength is also used to get the permittivity of dispersive materials). So it seems to me that it will be a bit of an internal mess to have this.

For next patch:

• Class method for auto mesh spec from dl0 (see above)
• MeshSpec.refine or some other way to add a meshing refinement structure.
• Import/export grid boundaries (or steps) so the same grid can be used in different simulations if needed.

[momchil-flex]

Here's a list of all docs notebooks. We need to incorporate the new auto nonuniform grid in all of them. You can do something like

sim.plot(z=0, ax=ax)
sim.plot_grid(z=0, ax=ax)

To have a look at the grid. In most cases the grid will be too dense to show, so if everything looks good, remove the grid plotting. But if you find some cases where it's useful to keep it on, do so!

• Dispersion.ipynb
• Fitting.ipynb
• GDS_import.ipynb
• GratingCoupler.ipynb
• HighQ_Ge.ipynb
• HighQ_Si.ipynb
• L3_cavity.ipynb - Momchil
• Metalens.ipynb - Momchil
• Modal_sources_monitors.ipynb
• Modes_bent_angled.ipynb - Momchil
• ModeSolver.ipynb - Momchil
• Near2Far.ipynb
• Near2FarSphereRCS.ipynb
• ParameterScan.ipynb
• RingResonator.ipynb
• Simulation.ipynb
• SMatrix.ipynb
• StartHere.ipynb
• VizData.ipynb - Momchil
• VizSimulation.ipynb - Momchil
• WebAPI.ipynb - no need
• WhatsNew.ipynb - Momchil

[shashwat-sh]

I have gone through and updated all the notebooks whose check boxes are ticked above. Here are some observations:

• Every time I ran a simulation I got the Could not decode response json error and Could not get estimated cost! warning. The simulations ran fine, but these errors are printed on the console and we may not want them to be visible on the docs, so that may mean having to run all sims once again before releasing the docs, after the errors are fixed.

• For Dispersion.ipynb, I've used a uniform mesh based on the original one, while the nonuniform case gets sorted out.

• For the Lumerical and Comsol benchmarks, I also used a uniform mesh based on the original, because the resolution was picked based on a periodicity criterion that I didn't want to mess with.

• For GratingCoupler.ipynb, which I don't think has any dispersive materials, I had trouble getting the results to match the original ones with a uniform mesh. Here are the cases I compared:

  1. Original (pre-GridSpec): dl set based on 50 cells per wavelength
  2. GridSpec.auto: 40 min_steps_per_wvl (takes much longer to run, I'm now trying 50)
  3. GridSpec.uniform: dl set based on 50 cells per wavelength, to reproduce original

  Note that I had trouble not only with ii, but also with iii, i.e. I was not able to reproduce the original results (I'm guessing the difference is because the new uniform mesh is still different from the old one in terms of snapping to domain size etc.). The deviation is not huge, for example, in the original, the flux in waveguide / flux in computed at the end was 5.15 %, whereas now it's 4.5 % in both ii and iii. The measured angle for the scattered cross section is also away from the design angle by an additional degree (~28 degrees before, ~27 +/- 0.5 degrees now). I'm not sure if these differences are large enough to be of concern, but I still find it a bit surprising. I got similar results when I used ii but with fewer cells per wavelength, it did not look like results were converging as I refined (but it was not a very detailed study so I could be wrong about that)