Automatical mesh refinement in layered structures

[weiliangjin2021]

(associated with https://github.com/flexcompute/tidy3d-core/issues/509)

Background

In many applications, structures are built layer by layer. Materials can be considered uniform along the normal axis.

• When there are metals (including PEC and lossy metal) present in the layer, it's important to snap and refine grid to the metal corners. We also want to refine around the corners/edges.
• Sometimes we want to resolve layer thickness. In most cases, we want lower or upper bound of the layer (or both) to be at grid.

In 2D, we can simply use shapely to find out corners. With corners, we can have boundaries (or edges). With the two, we can automatically generate snapping points and refinement mesh override structures.

Features

• A dedicated file for 2d corner finder;
• A GridRefinement class here for local mesh refinement. It is to make it easier to generate mesh override structures based on the number of grid cells and refined grid size. The refined grid size can be defined by either the refinement factor w.r.t. vacuum grid size, or grid size directly.
• LayerRefinementSpec that defines the layer (it can be bounded on the inplane dimensions so that corners outside the bound are discarded) and various mesh generation specs.
• Added layer_refinement_specs field to GridSpec. The final snapping points and override structures will be internal ones + user supplied ones. The internal ones are generated by LayerRefinementSpec.
• When dl_min is not set by the user, we'll automatically generate a value based on the structures, simulation domain size, and layer_specs.

Examples

If we are fine with default meshing specs, layer_refinement_spec can be easily defined. First, let's defined a layer that is unbounded perpendicular to the normal axis:

layer = td.LayerRefinementSpec.from_unbounded(
    axis=2,
    bounds=(-patch_thickness/2, patch_thickness/2),
)

and grid_spec,
grid_spec = td.GridSpec.auto(wavelength=wavelength, layer_refinement_specs=[layer])

Consider a patch antenna but we drill a big hole on the left, with the above grid_spec, we get

As you can see, metal corners are all identified, and snapping points are applied at those corners. Small mesh refinement structures are placed at corners, too, so that the mesh is denser near those singularities.

Sometimes, you only want refine a portion of the layer, as the other portion is not affecting the simulation as much. E.g., we only want to refine the feeding region,

layer = td.LayerRefinementSpec(
    axis=2,
    center=(0, -10e3, 0),
    size=(10e3,10e3, patch_thickness),
)

Tasks

Major

• Layer refinement specification class to define layer position and refinement parameters: defined here
• Corner finder in 2d with shapely
• Automatic snapping points generation
• Automatically sets dl_min if not set yet:
  • Compute dl_min from layers
  • Compute dl_min from all structures
• Automatic mesh refinement structures generation

Minor

• Allow snapping points to take None as coordinate, so that it's not snapped along that axis. Now snapping point is of type CoordinateOptional. In visualization, snapping points will be snapping line if one of the inplane coordinate is None.
• Shift _filter_structures_plane to geometry/utils.py, as it'll be used in 2d corner finder.

Miscellaneous

• Rebase once QuasiUniformGrid and min_steps_per_sim_size` #2130 is merged. Many conflicts expected.
• Add tests
• Add CHANGELOG
• Better visualization of snapping points that contain None

[weiliangjin2021]

There are still a few items to be done, but they are likely to be minor. You are welcome to start reviewing now, e.g. the general flow, naming, etc.

[dbochkov-flexcompute]

Looks great! Just some minor comments so far

[momchil-flex]

This seems like impressive work, note I mostly reviewed the usage rather than the low level since there seem to be many eyes on it.

[momchil-flex]

I assume this is not always enforced as it could lead to very small step size for thin layers?

[weiliangjin2021]

It will not be enforced if dl_min is set. However, if dl_min=0/None, we'll internally set a dl_min based on material response and those refinement factors, meaning that it will be enforced for thin layers. In a way, that is expected because in many classes of simulation you'll want to resolve the thin layer.

[momchil-flex]

Ok, I can see how that's probably what's needed in RF. However one thing I'm a little wary about is if people put some arbitrary small thickness for layers that they want to say are "very thin". Like ideally they should put 0 thickness and use a 2D medium but they may try instead to put like a 10nm plate with a layer refinement because they don't know better, and then get a giant number of grid cells and time steps.

[lucas-flexcompute]

It's not easy to differentiate between the 2 scenarios (intended vs. non-intended refinement). Maybe we could issue a warning if the simulation is in optical frequencies, but I'm not super convinced it would be useful, since we already have validations in place for the number of time steps (and cost)…

[weiliangjin2021]

like ideally they should put 0 thickness and use a 2D medium but they may try instead to put like a 10nm plate with a layer refinement because they don't know better

Probably low chance, as now by default, thickness will not be refined, and we only snap one bound to the grid.

[lucas-flexcompute]

This looks good. From what I can tell, we can use the interface easily enough to build automatic meshes for RF. It would be good to have a more clear understanding of the bounds_refinement and corner_refinement fields, i.e., how exactly do they work (maybe with an illustration?)

Am I understanding it correctly that the bounds_refinement will create GridRefinement.num_cells at both bounds? Are the cells inside the layer or outside the layer?

And for corner_refinement, where are the cells created? Inside or outside the polygon?

[dmarek-flex]

When a metal layer is very thin and we do not want to mesh both sides of the metal, so as to enable a larger delta t, would it be possible to automatically choose the best interface (top or bottom) for placing a grid boundary. Consider the case attached. Is it correct that placing a grid boundary along the black line would be more economical/accurate than at the red line, since it would reduce the number of times a control volume contains >2 materials? (Near the metal corners is where I am thinking)

[weiliangjin2021]

Am I understanding it correctly that the bounds_refinement will create GridRefinement.num_cells at both bounds? Are the cells inside the layer or outside the layer?

Yes, at both bounds. The refinement is applied in the region that centered at the bound, so it's both inside and outside. Do you see a need only outside/inside?

And for corner_refinement, where are the cells created? Inside or outside the polygon?

They are centered at the vertices, e.g. those magenta dashed boxes below

[lucas-flexcompute]

Am I understanding it correctly that the bounds_refinement will create GridRefinement.num_cells at both bounds? Are the cells inside the layer or outside the layer?

Yes, at both bounds. The refinement is applied in the region that centered at the bound, so it's both inside and outside. Do you see a need only outside/inside?

No, that's perfect!

[weiliangjin2021]

When a metal layer is very thin and we do not want to mesh both sides of the metal, so as to enable a larger delta t, would it be possible to automatically choose the best interface (top or bottom) for placing a grid boundary. Consider the case attached. Is it correct that placing a grid boundary along the black line would be more economical/accurate than at the red line, since it would reduce the number of times a control volume contains >2 materials? (Near the metal corners is where I am thinking)

When the thickness of metal layer is much thinner than the Yee cell size, it looks chances are equal with black/red line? Could you share a figure on the control volume on this (let's take this discussion off-line)?

Separately, it is interesting to study the behavior of conformal meshing in a simple microstrip when the grid line is placed at black/red. I think they will lead to slightly different averaged permittivity/permeability, as the red line go through 2 materials, while black 3 materials.

[dmarek-flex]

This PR is pretty great! I am still playing around with all of the options to understand what they are doing. I found perhaps one non-ideal scenario. When I used a default bounds_refinement=td.GridRefinement() I get an overrefined grid. The calculated dl for the two created MeshOverrideStructure is ~0.6 mm, but the resulting grid is much more refined.

Could it be because the two MeshOverrideStructure associated with the bounds of the layer are offset by a small amount?

Second pic is when bounds_refinement=None.

[dmarek-flex]

As I said, this is great! Still trying it on examples, but I think I found one bug.

[weiliangjin2021]

I found perhaps one non-ideal scenario. When I used a default bounds_refinement=td.GridRefinement() I get an overrefined grid. The calculated dl for the two created MeshOverrideStructure is ~0.6 mm, but the resulting grid is much more refined.

Could you share the code here so that I can reproduce?

[weiliangjin2021]

In plot_grid, to distinguish between lines from override structures and snapping lines, now I set the former to be dotted, and latter to be dashed. Here is an example on the visualization effects (the dash-dotted is coming from the overlap of the two types):