Convenience for setting simulation run time

[momchil-flex]

Is your feature request related to a problem? Please describe.
It is not always clear how to set a good simulation run_time.

We could do something like this: allow run_time to accept a float or a new object RunTimeSpec, with something like

class RunTimeSpec(Tidy3dBaseModel):
    """Stores specifications on how to automatically compute simulation run time.
    """
    source_time: SourceTime = pd.Field(
        ..., description="Source pulse to use. We could default it to the first source in the list of sources,"
        "but we have to make sure empty list doesn't error. "
    )

    propagation_length: float = pd.Field(
        None, description="Longest optical path in the simulation; defaults to longest simulation axis."
    )
    
    highest_q: float = pd.Field(
        None, description="Q-factor of the longest-lived resonance; defaults to 1, i.e. non-resonant."
    )

For e.g. PIC simulations, the user can do something like run_time=RunTimeSpec() and the default values should just work for most simulations. For other things (e.g. metalens), it still gives them an option to get a good run time estimate with much less effort than currently.

[tylerflex]

we’d need to come up with an algorithm for setting the run time based on various properties of the RunTimeSpec / Simulation. Maybe we can just work that detail out before breaking ground on any feature?

as for the API i think it would make sense to change

Simulation.run_time : Union[float, RunTimeSpec]

and then add a @property

Simulation._run_time

that computes run time based on the Simulation.run_time + other information (such as Simulation.size, Simulation.medium, Simulation.sources).
then it should be a matter of replacing the backend to use _run_time instead of run_time?

Alternatively, instead of either float or RunTimeSpec, we could add Simulation.run_time_spec and do some kind of handling for if both run_time and run_time_spec are defined.

[momchil-flex]

Yep, I think the first API option you describe sounds good.

Regarding the formula, it should be something related to what we discussed in the adjoint issue. Something like

source_decay_time = C / source_time.fwidth
propagation_time = propagation_length / C_0 / highest_material_index_at_source_time_freq0
propagation_time *= highest_q
run_time = source_decay_time + propagation_time

I think @lucas-flexcompute mentioned he's doing something similar, could you comment any feedback or tips?

[lucas-flexcompute]

Yes, I have something similar, but I default the highest Q to 100 and leave the refractive index at 1. The propagation length I get from the largest bounding box dimension.

fwidth = max(fmax - fmin, 0.1 * frequencies.mean())
source_time = 10 / fwidth
q_factor = 100
propagation_lenght = (bounds[1] - bounds[0]).max()
propagation_time = propagation_lenght / C_0 * q_factor
run_time = source_time + propagation_time

[tylerflex]

Thanks @lucas-flexcompute and this works pretty well?
I'm wondering about the idea behind fwidth = max(fmax - fmin, 0.1 * frequencies.mean())
, are the fmax and fmin from the monitors or the sources? what about frequencies?

Working towards some kind of spec that we can capture the free parameters of this formula with (when added to a Simulation). Right now it seems the spec would mainly contain

class RunTimeSpec(Tidy3DBaseModel):
    q_factor : pd.PositiveFloat = 100
    fwidth_factor1(?) : 0.1 (?)
    fwidth_factor2(?) : 10 (?)

and then in Simulation.run_time, we'd do something like

@property
def _run_time(self):
    run_time_spec = self.run_time_spec

    # either this
    fwidth = max(fmax - fmin, run_time_spec.fwidth_factor * frequencies.mean())
    source_time = run_time_spec.fwidth_factor2 / fwidth

    # or this
    source_time = self.run_time_of_longest_source(...)

    propagation_lenght = (self.bounds[1] - bounds[0]).max()
    propagation_time = propagation_lenght / C_0 * run_time_spec.q_factor

    return source_time + propagation_time

Do you have any thoughts on whether we should use your formula for source time or try to grab it from the Simulation sources? (note, I guess we can't reliably use the source_time.fwidth, we might also need to take offset into account.)

[lucas-flexcompute]

Those values and formulas for fwidth were taken from the smatrix plugin, to match the Gaussian pulse used there. Those seem to be reasonable defaults for common integrated devices. In the general case, I think grabbing from the source would be a better idea. In fact, each time source could provide it's own implementation of source.end_time() to be used here.

The estimation usually works with a good margin (shutoff around 30% of the total run time) for passive non-resonant components. It does not work for high Q resonators, but that's expected.

[tylerflex]

I see, so perhaps with that in mind, this could be a modification of the _run_time property?

@property
def _run_time(self):
    """The run time of the simulation (if `run_time` is not supplied)`."""
    run_time_spec = self.run_time_spec

    # contribution from the time of of the source pulses
    source_times = [src.source_time.end_time() for src in self.sources]
    if not source_times:
        raise ...
    source_time_max = np.max(source_times)
    source_time = run_time_spec.source_factor * source_time_max

    # contribution from field decay out of the simulation
    propagation_length = np.max((np.array(self.bounds[1]) - np.array(self.bounds[0])))
    maximum_n = self.maximum_refractive_index
    propagation_time = run_time_spec.q_factor * maximum_n * propagation_length / C_0

    return source_time + propagation_time

Note:

1. I added a free parameter run_time_spec.source_factor that roughly controls how much we want to wait for the source to decay.
2. I included maximum refractive index in the propagation time (maybe overkill? or can be already captured somewhat by the q_factor)?

[lucas-flexcompute]

It's looking good to me. In my opinion, if it's easy to get the maximum index, it's better to use it. I don't because I don't want photonforge to depend on any internals of Tidy3D that are not strictly required.

[tylerflex]

Cool. Regarding maximum index, I'm just slightly worried that including that AND Q factor might lead to some overestimation. I assume the two things are kind of linked so not sure. Also would a default Q of 100 be overkill?

[lucas-flexcompute]

I don't think 100 is overkill. That's equivalent to a very bad resonator, which I chose because it helps if there are unintended resonances in the structure (think of an MMI with multiple back-reflections or small cavities that appear in adjoint optimizations). I wouldn't consider anything below 500 a real resonator in optics (RF/MW are another story).