Simulations using COMSOL
========================

Once the design has been completed using *Qiskit-Metal*, the design object can be used to run simulations in COMSOL. Currently, *SQDMetal* supports:

- Capacitance matrix simulations
- RF s-parameter simulations

The appropriate simulation parameters are generated by *SQDMetal*. The resulting simulation can either be run by opening COMSOL via the generated mph file or directly from the Python environment.

Preparing the COMSOL simulation
-------------------------------

The API utilises the ``mph`` Python package. In order to run the commands below, **COMSOL must be installed**.

.. code-block:: python

    from SQDMetal.COMSOL.Model import COMSOL_Model
    from SQDMetal.COMSOL.SimCapacitance import COMSOL_Simulation_CapMats
    from SQDMetal.COMSOL.SimRFsParameter import COMSOL_Simulation_RFsParameters

    #Initialise the COMSOL engine (needs to only be run once)
    COMSOL_Model.init_engine()

    #Instantiate a COMSOL model
    cmsl = COMSOL_Model('leModel')

    #Create simulations to setup - in this case capacitance matrix and RF s-parameter
    sim_capMats = COMSOL_Simulation_CapMats(cmsl)
    sim_sParams = COMSOL_Simulation_RFsParameters(cmsl, adaptive='Multiple', modal_min_freq_num_eigs=(5.2e9,25))

    #(A) - Initialise model from Qiskit-Metal design object: design
    cmsl.initialize_model(design, [sim_capMats, sim_sParams], bottom_grounded=True)

    #(B) - Add metallic layers (use 0 if selecting ground plane)
    cmsl.add_metallic(1, threshold=1e-10, fuse_threshold=1e-10)
    cmsl.add_metallic(2, evap_mode='separate_delete_intersections', group_by_evaporations=True)
    #Alternative syntax to combine multiple layers (via a list) in one go. In doing so, smooth_radius can be used to round off the corners
    cmsl.add_metallic([3,4], threshold=2e-9, fuse_threshold=2e-9, smooth_radius=0.5e-6, multilayer_fuse=True)

    #(C) - OPTIONAL: add ground plane
    cmsl.add_ground_plane()

    #(D) - Combine all contiguous metals into a single metallic electrodes
    cmsl.fuse_all_metals()

    #(E) - OPTIONAL: Fix the indices of the different electrodes via components on those electrodes
    cmsl.reorder_conds_by_comps(['block0', 'launch0', 'launch1', 'launch2'])

    #(F) - OPTIONAL: Set everything within a rectangle to be finely meshed
    cmsl.fine_mesh_in_rectangle(-3e-3, -140e-6, -2.9e-3, 167e-6, 2e-7, 5e-6)
    #(G) - OPTIONAL: Set conductors within a rectangle to be finely meshed
    cmsl.fine_mesh_conductors_in_rectangle(-2.86e-3, -145e-6, -2.53e-3, 167e-6, 2e-7, 5e-6)

    #(H) - Build model
    cmsl.build_geom_mater_elec_mesh(skip_meshing=True, mesh_structure='Fine')

    #(I) - Plot model
    cmsl.plot()

    #(J) - Save model to file
    cmsl.save('Test')

The simulation treats the electrodes as 2D planar objects that live on top of a 3D block representing the dielectric subtrate chip. The chip is enclosed by a rectangular prism that represents the exterior boundaries of the simulation. The simulation automatically accounts for changes in the device geometry due to [shadow evaporations](PVD.md).

Note the following features:
- (A) - To initialise the model, one must provide:
    - *Qiskit-Metal* design object
    - List of COMSOL simulation objects (they all implement `COMSOL_Simulation_Base`).
    - Optional `pad_x`, `pad_y` and `pad_z` (defaults to 0.5mm) that provide the padding from the chip to the exterior boundaries along the x, y and z axes respectively.
    - Optional `bottom_grounded` parameter sets the z-padding below the chip to be zero. That is, the chip will sit flush on the exterior boundary on the bottom; which, for all intents and purposes grounds the bottom of the chip.
- (B) - The metallic polygons are added to the top of the chip as planar objects. They are added per layer in the *Qiskit-Metal* design object. Note that contiguous metallic objects are fused as single electrodes. One must provide:
    - Layer index corresponding to the *Qiskit-Metal* design object
    - Optional `threshold` (defaults to `-1` where it is disabled) given as a positive length in metres. Idea is that:
        - After fusing all contiguous metallic objects, the resulting polygons may have adjacent points that are too close; this causes unnecessarily fine meshing in COMSOL.
        - **If adjacent points in a given polygon electrode are within the given `threshold`, they are combined on the spot.**
        - Use this to clean up a polygon by setting a size that is smaller than the smallest actual feature size.
    - Optional `fuse_threshold` (defaults to `1e-12`) given as a positive length in metres. Idea is that:
        - The `shapely` package does not properly combine polygons in certain cases due to floating-point rounding errors.
        - **If adjacent polygon electrodes are within the given `fuse_threshold` of each other, they are combined into a single polygon.**
    - Optional `evap_mode` (defaults `'separate_delete_below'`) and `evap_trim` (defaults `20e-9`) on how different electrodes across different are to be combined or separated. The options are the same as discussed in the usage of `plot_layer` in the [`PVD_Shadows` class](PVD.md).
    - Optional `group_by_evaporations` (defaults `False`) is relevant when using `'separate_delete_intersections'` or `'separate_delete_below'` for `evap_mode`. If `True`, all angled evaporation steps that are separated due to some deletion, are grouped together as a single metallic electrode - as discussed later in capacitance matrix simulations.
- (C) - The ground plane is a metallic layer corresponding to the ground cut-outs/subtractions in *Qiskit-Metal*. Note that it currently does not account for any shadow evaporations for it is assumed that the fills are evaporated using single vertical evaporations.
- (D) - Metallic electrods are added layer-by-layer. In cases of overlap between layers, it is important to treat them as single electrodes contiguous electrodes. Note that:
    - The function `fuse_all_metals` does not modify the underlying geometry - i.e. overlapping polygons across multiple layers are NOT combined into single polygons.
    - Instead `fuse_all_metals` combines the polygonal selections used to define the electrodes in the COMSOL simulations.
- (E) - In cases where the simulation is swept, it is useful to fix the indices of different electrodes. This can be done deterministically via the `reorder_conds_by_comps` function. Note the following:
    - This function just requires a list of the names of Qiskit-Metal components that are attached to different electrodes. The electrodes are indexed according to the order given in this list. If the number of electrodes exceeds the number of supplied names, the remaining electrodes are automatically indexed.
    - The given component must be completely contiguous and **must exist entirely within a given electrode**. For example, a capacitor is not a good choice as it will split across two COMSOL electrodes (unless the two electrodes are electrically connected via other metallic polygons).
- (F) - Everything on top of the chip within a given rectangular region is meshed via custom fine-meshing parameters. The arguments are:
    - First four parameters are `x1`, `y1`, `x2` and `y2` where `(x1,y1)` is the bottom-left corner while `(x2,y2)` is the top-right corner of the rectangular region.
    - Optional `minElementSize` (defaults `1e-7`) sets the minimum mesh-element size of the triangular mesh in metres.
    - Optional `maxElementSize` (defaults `5e-6`) sets the maximum mesh-element size of the triangular mesh in metres.
- (G) - Instead of meshing everything finely, `fine_mesh_conductors_in_rectangle` is a counterpart to `fine_mesh_in_rectangle` (with the same arguments) in which only conductive regions falling within the rectangle are finely meshed with the given meshing parameters.
- (H) - Build the geometry, setup the materials (i.e. the substrate dielectric and the metals), setup the individual simulation studies and build the resulting mesh in COMSOL. One may provide:
    - Optional `skip_meshing` (defaults `False`). If `True`, the mesh is setup but not explicitly built.
    - Optional `mesh_structure` (defults `'Normal'`) to define the general size of the meshing in COMSOL. Can be: `'Extremely fine'`, `'Extra fine'`, `'Finer'`, `'Fine'`, `'Normal'`, `'Coarse'`, `'Coarser'`, `'Extra coarse'` or `'Extremely coarse'`.
    - Optional `substrate_permittivity` (defaults `11.7` for silicon) to set for the dielectric substrate chip.
- (I) - The 2D model plot shows the different electrodes (and their selection indices) that have been setup for simulation. It should be used to ensure that the electrodes are properly mapped (e.g. when tuning the different threshold parameters in `add_metallic`). The plot function may take the following arguments:
    - Optional `plot_chip_boundaries` (defaults `True`) to include the chip boundaries if `True`. Useful to check if the different metallic elements fit within the chip in the simulation.
    - Optional `plot_fine_mesh_bounds` (defaults `True`) to include the fine mesh boundaries if `True`.
- (J) - This command saves the current COMSOL model into an mph file; in this example, the file will be called `'Test.mph'`. Can include a valid folder path; do not include the mph extension as it will be automatically added.

Note that the order of operations given above must be generally satisfied. For example:
- `fuse_all_metals` should be called after adding all metallic regions
- `reorder_conds_by_comps` should be called after `fuse_all_metals`
- Almost everything (including the fine-meshing) should be called before `build_geom_mater_elec_mesh`.
- However, `save` may be called any time and the model in its current state will be saved.

The following sections will outline the different simulations, where one may directly run a given simulation from the Python environment.

Capacitance matrix simulations
------------------------------

Given a COMSOL model, the DC capacitances between different electrodes can be calculated via the `COMSOL_Simulation_CapMats` class. The indices of the different electrodes are assigned automatically or manually via the `reorder_conds_by_comps` command as outlined in the earlier section. The capacitance matrix can be extracted directly from the `COMSOL_Simulation_CapMats` as follows:

.. code-block:: python

    #Must call build_geom_mater_elec_mesh first and have sim_capMats registered when calling initialize_model
    capMat = sim_capMats.run()


Note that the exterior boundaries are considered to be ground and held as a zero potential reference. Thus, no electrode should intersect with the exterior boundaries as setting the unit potential will fail if the conductor touches a region of zero potential. 

RF s-parameter simulations
--------------------------

Given a COMSOL model, the s-parameters can be calculated between different ports. The `COMSOL_Simulation_RFsParameters` class is instantiated via the COMSOL model and some optional parameters:

- `adaptive` (defaults `'None'`) specifies the kind of simulation to setup:
    - `'None'` - the simulation is a frequency domain sweep where the solutions for every frequency sequentially. This gives the most accurate global phase, but it is the slowest simulation type.
    - `'Single'` - the simulation finds the resonant frequency and width and interpolates the frequency scan given this resonant peak. This is fast and accurate if it is guaranteed that there is only a single resonant peak of interest in the s-parameters.
    - `'Multiple'` - the simulation performs an eigenfrequency scan to find the resonant peaks and then interpolates the frequency scan given these peaks. This is fast and accurate if **all** resonant peaks within the given frequency range of interest are correctly found.
- `relative_tolerance` (defaults `0.01`) is only relevant when using `adaptive='Single'` and gives the relative simulation tolerance in the converged values regarding the resonant peak.
- `modal_min_freq_num_eigs` (defaults `(1e9,7)`) is only relevant when using `adaptive='Multiple'`:
    - Given as a tuple of the the minimum frequency (in Hertz) and number of expected eigenfrequency peaks.
    - Note that the minimum frequency should be less than or equal to the frequency range that is set for the simulation. All searched eigenfrequencies are assumed to be larger than this frequency.
    - The number of eigenfrequencies should be at least larger than the expected number of resonant peaks.


**Before running `build_geom_mater_elec_mesh`**, one must define the ports via any of the `create_port_` commands:


.. code-block:: python

    #(A) - Create ports on LaunchpadWirebond objects
    sim_sParams.create_port_CPW_on_Launcher('launch1')
    sim_sParams.create_port_CPW_on_Launcher('launch2')


The supported `create_port_` commands work as follows:
- (A) - `create_port_CPW_on_Launcher`:
    - Creates a CPW-fed port on a *Qiskit-Metal* `LaunchpadWirebond` component at the start of the launcher. The CPW feed extends onto the ground plane.
    - The function requires the name of said `LaunchpadWirebond` component in the *Qiskit-Metal* design object.
    - Optionally provide `len_launch` (defaults `20e-6`) to define the length of the CPW port inlet along the length of the launcher.

The order of the s-parameter ports are given in the order of calling the `create_port_` commands.

Once the ports are defined and the simulation has been built via `build_geom_mater_elec_mesh`, one may set the frequency range upon which to calculate the s-parameters:


.. code-block:: python

    #Set frequency scan via an linspace-like construct (1000 points between 2GHz and 5GHz in this case):
    sim_sParams.set_freq_range(2e9, 5e9, 1000)

    #Set frequency scan via absolute frequency values (given as a list or numpy array):
    sim_sParams.set_freq_values([2e9,2.5e9,3e9,5e9])


To run the simulation from the Python environment:


.. code-block:: python

    freq_s_params = sim_sParams.run()


The result is a numpy array where the rows are: frequencies, S11 values, S21 values, S31 values, ..., SN1 values (for N defined ports). That is the excitation currently is set to the first port.

Finally, note that on running the simulation (as shown below), all field solutions are stored for all frequencies except for the individual frequency solutions (i.e. Multi-Modal) when using `adaptive='Multiple'`, where only the fields on the ports (i.e. required for s-parameter calculations) are stored.

Magnetic field simulations
--------------------------

Given a COMSOL model, a magnetic field simulation can be done via the `COMSOL_Simulation_MagneticField` class. This simulation is specifically designed to send in a DC current, observe the resulting magnetic field and subsequently calculate the net flux through some unit area (e.g. a SQUID loop).

**Before running `build_geom_mater_elec_mesh`**, one must define the following:

.. code-block:: python

    sim_BField.set_current_feed_on_CPW_on_Route('flux_line')
    sim_BField.set_Bfield_integration_area(coords)

The first command attaches a CPW-friendly current feeder (i.e. feeds the centre line while grounding both adjacent ground-planes simultaneously). The second command defines an area to which the $B_z$ component of the field is integrated. Note that a handy function in the `PVD_Shadows` class can be used to find the coordinates of the largest interior (i.e. hole) in a SQUID loop after shadow evaporation calculations:

.. code-block:: python

    pvd = PVD_Shadows(design)
    # pvd.plot_layer(2)
    coords, poly = pvd.get_shadow_largest_interior_for_component('jj_qubit')

To get the flux, simply run:

.. code-block:: python

    net_flux = sim_BField.run()

The units are in **Webers** (i.e. flux is the surface integral of $\mathbf{B}$).


TODO: Update this upon documenting COMSOL classes