diff --git a/README_Gmsh_Elmer.md b/README_Gmsh_Elmer.md new file mode 100644 index 000000000..e495bddb8 --- /dev/null +++ b/README_Gmsh_Elmer.md @@ -0,0 +1,85 @@ +# Softwares to Install for Open-source simulation toolchain for Qiskit Metal +## Gmsh +- This is a design rendering and meshing software that we will be using to render our design and then divide it into fine elements (called a `mesh`), in order to solve the necessary differential equations (PDE's of Maxwell's Electromagnetic Laws). This is completely integrated with Qiskit Metal and is semi-automated. +- This will be automatically installed when you make a new `venv` environment and install qiskit metal using the instructions given above (except in Apple Silicon Macs, see the Note below). + +### Test if Gmsh is properly installed +- Run the following in your terminal. A blank Gmsh window should open up. If it does, then hurray, you did it! 😄 +```bash +$ gmsh +``` + +- Open a python REPL in your terminal and type `import gmsh` in the REPL. If `gmsh` gets imported without any errors, then it's installed properly with python. +```python +$ python +... +... +>>> import gmsh # No Error +>>> +``` + +**NOTE:** +- On Apple Silicon Macs (with M1 and M2-series chips), installing Gmsh using `pip install gmsh` might not work for importing the package using python. Therefore, to successfully install Gmsh, and point it to the `gmsh` installed using `pip`, you need to install `gmsh` binary through [Homebrew](https://brew.sh/) as follows: +```bash +$ brew install gmsh +``` +- Restart your terminal, activate the environment, and check if you're able to do `import gmsh` in a python REPL or a scratch jupyter notebook. If the import is successful, then congratulations, you've installed Gmsh successfully! + +- **Note:** If the above steps still gives error in importing gmsh using `import gmsh` in your python environment, refer [this issue](https://gitlab.onelab.info/gmsh/gmsh/-/issues/1705), on the Gmsh GitLab repository, for more information. Also please feel free to contact us through the Qiskit Slack workspace on the `#metal` channel. + +## ElmerFEM +- Congratulations on making it until here! Now we'll see how to install ElmerFEM. +- ElmerFEM is a Finite Element Method (FEM) solver that we'll be using to take the mesh generated by Gmsh, and the solve the Maxwell's equations on top of our meshed design. Right now, Qiskit Metal only supports capacitance-type simulations (Poisson equaiton solver provided by `StatElecSolver`) with ElmerFEM, so only LOM 2.0 analysis can be performed for tuning the qubits using the results obtained from ElmerFEM in the form of a Maxwell Capacitance matrix. +- ElmerFEM doesn't come as a python library and has different installation options on your operating system. +- Please follow the official guide provided by Elmer Foundation CSC, [here](https://github.com/ElmerCSC/elmerfem#elmer-fem) + +**NOTE:** For MacOS, please consider building the software from source rather than using Homebrew on Mac, as the homebrew install isn't very consistent and may cause issues in the future. + +### After you complete your installation +- Add your installation path directory to your system's `PATH` environment variable (you can skip this on Windows and MacOS, as it is automatically done on those OS platforms during the installation steps). +- Restart your terminal. +- Check if ElmerFEM is installed successfully by checking installation for `ElmerGrid` and `ElmerSolver`. +- For testing ElmerGrid, type `ElmerGrid` in your terminal and if the PATH variable is set correctly, it should print a long output text. Check if you see the following at the end of the text: +```bash +$ ElmerGrid +Starting program Elmergrid +****************** Elmergrid ************************ +This program can create simple 2D structured meshes consisting of +linear, quadratic or cubic rectangles or triangles. The meshes may +also be extruded and revolved to create 3D forms. In addition many +mesh formats may be imported into Elmer software. Some options have +not been properly tested. Contact the author if you face problems. + +... +... +... + +Thank you for using Elmergrid! +Send bug reports and feature wishes to elmeradm@csc.fi +``` + +- For testing ElmerSolver, type `ElmerSolver` in your terminal and if the PATH variable is set correctly, it should print a long output text. Check if you see the following: + +```bash +$ ElmerSolver +ELMER SOLVER (v 9.0) STARTED AT: 2023/01/11 09:42:44 +ParCommInit: Initialize #PEs: 1 +MAIN: OMP_NUM_THREADS not set. Using only 1 thread per task. +MAIN: +MAIN: ============================================================= +MAIN: ElmerSolver finite element software, Welcome! +MAIN: This program is free software licensed under (L)GPL +MAIN: Copyright 1st April 1995 - , CSC - IT Center for Science Ltd. +MAIN: Webpage http://www.csc.fi/elmer, Email elmeradm@csc.fi +MAIN: Version: 9.0 (Rev: 02f61d697, Compiled: 2022-11-08) +MAIN: Running one task without MPI parallelization. +MAIN: Running with just one thread per task. +MAIN: ============================================================= +ERROR:: ElmerSolver: Unable to find ELMERSOLVER_STARTINFO, can not execute. +Note: The following floating-point exceptions are signalling: IEEE_OVERFLOW_FLAG +STOP 1 +``` + +If you see the above lines in the text in your terminal, then CONGRATULATIONS!!! You've successfully installed everything needed for using Gmsh and ElmerFEM with Qiskit Metal! 😄 + +**NOTE:** If you face issues at any of the steps above, feel free to contact us on #metal channel on Qiskit slack!! diff --git a/README_developers.md b/README_developers.md index 754360aba..7b9fb8a86 100644 --- a/README_developers.md +++ b/README_developers.md @@ -115,6 +115,12 @@ Here are some things to consider when setting up a development environment: * Library errors when activating conda environments, or initializing jupyter notebook/lab, might indicate a conflict between python libraries in the base and sub environments. Go ahead and manually delete the library from the base environment `site-packages` folder, shown in the error message. You might need to reinstall them in the sub environment, or create a new one. * If Jupyter notebook has trouble finding a dll for a package that works in the new environment outside of Jupyter, then try opening Jupyter notebook from the new environment instead of from `base` +### Installing other dependencies for Open-source Renderers (Gmsh and ElmerFEM) + +If you want to use the recently added open-source renderers for [Gmsh](./qiskit_metal/renderers/renderer_gmsh) and [ElmerFEM](./qiskit_metal/renderers/renderer_elmer) for simulation of your design, please make sure that both of them have been installed successfully in your system before running Qiskit Metal. On Windows, Linux, and MacOS (with x86_64 architecture CPUs), Gmsh will be installed automatically using the `environment.yml` file during the conda installation step above. For more detailed steps to install ElmerFEM, please refer to [this](./README_Gmsh_Elmer.md) document. + +**NOTE:** We would like to give a disclaimer for users on Apple silicon Macs (M1 and M2-series). Currently, Qiskit Metal uses PySide2 which is not natively supported on the ARM architecture. This will lead to error in instantiating the `MetalGUI` as of now. However, if you still want to use Qiskit Metal without the GUI, the process for installing Gmsh software is a bit different and can be found in [this](./README_Gmsh_Elmer.md) document. + # Other Common Issues For other common installation issues, please refer to the [FAQ](https://qiskit.org/documentation/metal/faq.html) @@ -149,4 +155,4 @@ Go to directory with qiskit-metal/.style.yapf file and run the above command to rm /hooks/pre-commit ``` -If you need to uninstall the precommit hook, go to the root of the project and run the above command. \ No newline at end of file +If you need to uninstall the precommit hook, go to the root of the project and run the above command. diff --git a/environment.yml b/environment.yml index af8088ce2..3117af392 100644 --- a/environment.yml +++ b/environment.yml @@ -21,6 +21,7 @@ dependencies: - shapely==2.0.1 - jupyter==1.0.0 - scqubits==3.1.0 + - pyyaml==6.0 - pip==23.0 - pip: - pyaedt==0.6.46 diff --git a/qiskit_metal/config.py b/qiskit_metal/config.py index f4281f24a..1523e8ac6 100644 --- a/qiskit_metal/config.py +++ b/qiskit_metal/config.py @@ -30,6 +30,10 @@ class_name='QQ3DRenderer'), gds=Dict(path_name='qiskit_metal.renderers.renderer_gds.gds_renderer', class_name='QGDSRenderer'), + gmsh=Dict(path_name='qiskit_metal.renderers.renderer_gmsh.gmsh_renderer', + class_name='QGmshRenderer'), + elmer=Dict(path_name='qiskit_metal.renderers.renderer_elmer.elmer_renderer', + class_name='QElmerRenderer'), ) """ Define the renderes to load. Just provide the module names here. diff --git a/qiskit_metal/renderers/__init__.py b/qiskit_metal/renderers/__init__.py index 44c5ae6c0..e28c9dde0 100644 --- a/qiskit_metal/renderers/__init__.py +++ b/qiskit_metal/renderers/__init__.py @@ -107,4 +107,4 @@ from .renderer_ansys.q3d_renderer import QQ3DRenderer from .renderer_gmsh.gmsh_utils import Vec3DArray - from .renderer_gmsh.gmsh_renderer import QGmshRenderer \ No newline at end of file + from .renderer_gmsh.gmsh_renderer import QGmshRenderer diff --git a/qiskit_metal/renderers/renderer_elmer/__init__.py b/qiskit_metal/renderers/renderer_elmer/__init__.py new file mode 100644 index 000000000..a9753a67e --- /dev/null +++ b/qiskit_metal/renderers/renderer_elmer/__init__.py @@ -0,0 +1,14 @@ +# -*- coding: utf-8 -*- + +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2022. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" ElmerFEM QRenderer """ \ No newline at end of file diff --git a/qiskit_metal/renderers/renderer_elmer/elmer_configs.py b/qiskit_metal/renderers/renderer_elmer/elmer_configs.py new file mode 100644 index 000000000..59b36ac32 --- /dev/null +++ b/qiskit_metal/renderers/renderer_elmer/elmer_configs.py @@ -0,0 +1,74 @@ +"""Configuration settings used by ElemerFEM. For details, please refer +to the ElmerFEM documentation http://www.elmerfem.org/blog/documentation/ +""" + +# Dictionary containing a database of materials and their properties to be +# assigned to an ElmerFEM model. +materials = dict( + vacuum={ + "Electric Conductivity": 0.0, + "Relative Permittivity": 1, + }, + silicon={ + "Relative Permittivity": 11.45, + }, +) + +# Dictionary containing valid ElmerFEM simulation types and dimensions with +# their correspoding settings (e.g. steady state 3D, transient 2D, etc.). +simulations = dict( + steady_3D={ + "Max Output Level": 5, + "Coordinate System": "Cartesian", + "Coordinate Mapping(3)": "1 2 3", + "Simulation Type": "Steady state", + "Steady State Max Iterations": 1, + "Output Intervals": 1, + "Timestepping Method": "BDF", + "BDF Order": 1, + "Solver Input File": "case.sif", + "Post File": "case.ep", + "Output File": "case.result", + }) + +# Dictionary containing valid ElmerFEM solvers and their correspoding +# settings (e.g. StatElecSolver, ResultOutputSolve, etc.). +solvers = dict( + capacitance={ + "Equation": "Electrostatics", + "Calculate Electric Field": True, + "Calculate Capacitance Matrix": True, + "Capacitance Matrix Filename": "capacitance.txt", + "Procedure": '"StatElecSolve" "StatElecSolver"', + "Variable": "Potential", + "Calculate Electric Energy": True, + "Exec Solver": "Always", + "Stabilize": True, + "Bubbles": False, + "Lumped Mass Matrix": False, + "Optimize Bandwidth": True, + "Steady State Convergence Tolerance": 1.0e-5, + "Nonlinear System Convergence Tolerance": 1.0e-7, + "Nonlinear System Max Iterations": 20, + "Nonlinear System Newton After Iterations": 3, + "Nonlinear System Newton After Tolerance": 1.0e-3, + "Nonlinear System Relaxation Factor": 1, + "Linear System Solver": "Iterative", + "Linear System Iterative Method": "BiCGStab", + "Linear System Max Iterations": 500, + "Linear System Convergence Tolerance": 1.0e-10, + "BiCGstabl polynomial degree": 2, + "Linear System Preconditioning": "ILU0", + "Linear System ILUT Tolerance": 1.0e-3, + "Linear System Abort Not Converged": False, + "Linear System Residual Output": 10, + "Linear System Precondition Recompute": 1, + }, + postprocessing_gmsh={ + "Exec Solver": "Always", + "Equation": "Result Output", + "Procedure": '"ResultOutputSolve" "ResultOutputSolver"', + "Output File Name": '"case"', + "Output Format": "gmsh", + }, +) diff --git a/qiskit_metal/renderers/renderer_elmer/elmer_renderer.py b/qiskit_metal/renderers/renderer_elmer/elmer_renderer.py new file mode 100644 index 000000000..b35fcbc43 --- /dev/null +++ b/qiskit_metal/renderers/renderer_elmer/elmer_renderer.py @@ -0,0 +1,769 @@ +from typing import Union, Optional +import os +import pandas as pd + +from qiskit_metal import Dict, draw +from qiskit_metal.renderers.renderer_base import QRendererAnalysis +from qiskit_metal.renderers.renderer_gmsh.gmsh_renderer import QGmshRenderer +from qiskit_metal.renderers.renderer_elmer.elmer_runner import ElmerRunner + + +def load_capacitance_matrix_from_file(filename: str) -> pd.DataFrame: + """Loads capacitance matrix from file. + + Args: + filename (str): A '.txt' file containing the capacitance matrix. + + Returns: + pd.DataFrame:: Table containing capacitance matrix. + """ + return pd.read_csv(filename, delimiter=' ', index_col=0) + + +class QElmerRenderer(QRendererAnalysis): + """Extends QRendererAnalysis class and imports meshes from Gmsh using the ElmerFEM python API. + + QElmerRenderer Default Options: + * simulation_type -- Type of ElmerFEM simulation to run. + Default options found in elmer_configs.py under "simulations". + * simulation_dir -- Directory to store simulation-related files. + * mesh_file -- Name of GMSH output mesh file ending in '.msh'. + * simulation_input_file -- Name of ElmerFEM simulation input file ending in '.sif'. + * postprocessing_file -- Name of ElmerFEM postprocessing output file ending in '.msh'. + * output_file -- Name of ElmerFEM results output file ending in '.result'. + + QElmerRenderer Default Setup: + * capacitance -- Default setup for capacitance simulation: + * Calculate_Electric_Field -- Defaults to True. + * Calculate_Electric_Energy -- Defaults to True. + * Calculate_Capacitance_Matrix -- Defaults to True. + * Capacitance_Matrix_Filename -- Defaults to "cap_matrix.txt". + * Linear_System_Solver -- Defaults to "Iterative". + * Steady_State_Convergence_Tolerance -- Defaults to 1.0e-5. + * Nonlinear_System_Convergence_Tolerance -- Defaults to 1.0e-7. + * Nonlinear_System_Max_Iterations -- Defaults to 20. + * Linear_System_Convergence_Tolerance -- Defaults to 1.0e-10. + * Linear_System_Max_Iterations -- Defaults to 500. + * Linear_System_Iterative_Method -- Defaults to "BiCGStab". + * BiCGstabl_polynomial_degree -- Defaults to 2. + + * eigenmode -- Default setup for Eigenmode simulation (to be added in future). + + * materials -- Materials used in the simulation. Defaults to `["vacuum", "silicon"]`. + + * constants -- Default constants to be included in the simulation: + * Permittivity_of_Vacuum -- Defaults to 8.8542e-12 F/m. + * Unit_Charge -- Defaults to 1.602e-19 C. + + QGmshRenderer Default Options: + # Buffer between max/min x and edge of ground plane, in mm + * x_buffer_width_mm -- Buffer between max/min x and edge of ground plane, in mm + * y_buffer_width_mm -- Buffer between max/min y and edge of ground plane, in mm + * mesh -- to define meshing parameters + * max_size -- upper bound for the size of mesh node + * min_size -- lower bound for the size of mesh node + * max_size_jj -- maximum size of mesh nodes at jj + * smoothing -- mesh smoothing value + * nodes_per_2pi_curve -- number of nodes for every 2π radians of curvature + * algorithm_3d -- value to indicate meshing algorithm used by Gmsh + * num_threads -- number of threads for parallel meshing + * mesh_size_fields -- specify mesh size field parameters + * min_distance_from_edges -- min distance for mesh gradient generation + * max_distance_from_edges -- max distance for mesh gradient generation + * distance_delta -- delta change in distance with each consecutive step + * gradient_delta -- delta change in gradient with each consecutive step + * colors -- specify colors for the mesh elements, chips or layers + * metal -- color for metallized entities + * jj -- color for JJs + * dielectric -- color for dielectric entity + """ + + default_options = Dict( + simulation_type="steady_3D", + simulation_dir="./simdata", + mesh_file="out.msh", + simulation_input_file="case.sif", + postprocessing_file="case.msh", + output_file="case.result", + ) + + default_setup = Dict( + capacitance=Dict( + Calculate_Electric_Field=True, + Calculate_Electric_Energy=True, + Calculate_Capacitance_Matrix=True, + Capacitance_Matrix_Filename="cap_matrix.txt", + Linear_System_Solver="Iterative", + Steady_State_Convergence_Tolerance=1.0e-5, + Nonlinear_System_Convergence_Tolerance=1.0e-7, + Nonlinear_System_Max_Iterations=20, + Linear_System_Convergence_Tolerance=1.0e-10, + Linear_System_Max_Iterations=500, + Linear_System_Iterative_Method="BiCGStab", + BiCGstabl_polynomial_degree=2, + ), + + # TODO: Update this later when working on eigenmode setup + eigenmode=Dict(), + materials=["vacuum", "silicon"], + constants=Dict( + Permittivity_of_Vacuum=8.8542e-12, + Unit_Charge=1.602e-19, + ), + ) + + name = "elmer" + """Name""" + + def __init__(self, + design: 'MultiPlanar', + layer_types: Union[dict, None] = None, + initiate=True, + options: Dict = None): + """ + Args: + design ('MultiPlanar'): The design. + layer_types (Union[dict, None]): the type of layer in the format: + dict(metal=[...], dielectric=[...]). + Defaults to None. + initiate (bool): True to initiate the renderer (Default: False). + settings (Dict, optional): Used to override default settings. Defaults to None. + """ + default_layer_types = dict(metal=[1], dielectric=[3]) + self.layer_types = default_layer_types if layer_types is None else layer_types + + super().__init__(design=design, initiate=initiate, options=options) + + @property + def initialized(self): + """Abstract method. Must be implemented by the subclass. + Is renderer ready to be used? + Implementation must return boolean True if successful. False otherwise. + """ + return (hasattr(self, "_elmer_runner") and self.gmsh.initialized) + + def initialize_renderer(self): + """Initializes the Gmsh and Elmer renderer. + NOTE TO USER: this should be used when using + the QElmerRenderer through design.renderers.elmer instance. + + Example: To change one of the options exposed by QGmshRenderer, + one has to first initiate the QElmerRenderer using this method + and then change the property like so, + + ``` + design.renderers.elmer.initialize_renderer() + design.renderers.elmer.gmsh.options.mesh.num_threads = 1 + ... + ``` + """ + self.gmsh = QGmshRenderer(self.design, self.layer_types) + self._elmer_runner = ElmerRunner() + + def _initiate_renderer(self): + """Initializes the Gmsh and Elmer renderer. + NOTE: This is automatically called when the USER + specifically imports QElmerRenderer in a jupyter notebook + and instantiates it.""" + self.gmsh = QGmshRenderer(self.design, self.layer_types) + self._elmer_runner = ElmerRunner() + return True + + def _close_renderer(self): + """Finalizes the Gmsh renderer""" + self.gmsh.close() + + def close(self): + """Public method to close the Gmsh renderer""" + return self._close_renderer() + + def render_design( + self, + selection: Union[list, None] = None, + open_pins: Union[list, None] = None, + box_plus_buffer: bool = True, + draw_sample_holder: bool = True, + skip_junctions: bool = False, + mesh_geoms: bool = True, + ignore_metal_volume: bool = False, + omit_ground_for_layers: Optional[list[int]] = None, + ): + """Render the design in Gmsh and apply changes to modify the geometries + according to the type of simulation. Simulation parameters provided by the user. + + Args: + selection (Union[list, None], optional): List of selected components + to render. Defaults to None. + open_pins (Union[list, None], optional): List of open pins that are open. + Defaults to None. + box_plus_buffer (bool, optional): Set to True for adding buffer to + chip dimensions. Defaults to True. + draw_sample_holder (bool, optional): To draw the sample holder box. Defaults to True. + skip_junctions (bool, optional): Set to True to sip rendering the + junctions. Defaults to False. + mesh_geoms (bool, optional): Set to True for meshing the geometries. + Defaults to True. + ignore_metal_volume (bool, optional): ignore the volume of metals and replace + it with a list of surfaces instead. + Defaults to False. + omit_ground_for_layers (Optional[list[int]]): omit rendering the ground plane for + specified layers. Defaults to None. + """ + + # For handling the case when the user wants to use + # QElmerRenderer from design.renderers.elmer instance. + if not self.initialized: + self._initiate_renderer() + + self.gmsh.render_design(selection=selection, + open_pins=open_pins, + box_plus_buffer=box_plus_buffer, + draw_sample_holder=draw_sample_holder, + skip_junctions=skip_junctions, + mesh_geoms=mesh_geoms, + ignore_metal_volume=ignore_metal_volume, + omit_ground_for_layers=omit_ground_for_layers) + + self.qcomp_geom_table = self.get_qgeometry_table() + self.nets = self.assign_nets(open_pins=open_pins) + + def get_qgeometry_table(self) -> pd.DataFrame: + """Combines the "path" and "poly" qgeometry tables into a single table, and adds + column containing the minimum z coordinate of the layer associated with + each qgeometry. + + Raises: + ValueError: Raised when component in selection of qcomponents is not in QDesign. + + Returns: + pd.DataFrame: Table with elements in both "path" and "poly" qgeometry tables. + """ + + if self.gmsh.case == 0: + qcomp_ids = self.gmsh.qcomp_ids + elif self.gmsh.case == 1: + qcomp_ids = list(set(self.design._components.keys())) + elif self.gmsh.case == 2: + raise ValueError("Selection provided is invalid.") + + metal_layers = self.layer_types["metal"] + + mask = lambda table: table["component"].isin(qcomp_ids) & ~table[ + "subtract"] & table["layer"].isin(metal_layers) + + min_z = lambda layer: min( + sum(self.gmsh.get_thickness_zcoord_for_layer_datatype(layer)), + self.gmsh.get_thickness_zcoord_for_layer_datatype(layer)[1]) + + path_table = self.design.qgeometry.tables["path"] + poly_table = self.design.qgeometry.tables["poly"] + qcomp_paths = path_table[mask(table=path_table)] + qcomp_polys = poly_table[mask(table=poly_table)] + qcomp_geom_table = pd.concat([qcomp_paths, qcomp_polys], + ignore_index=True) + + qcomp_geom_table['min_z'] = qcomp_geom_table['layer'].apply(min_z) + qcomp_geom_table.sort_values(by=['min_z'], + inplace=True, + ignore_index=True) + + return qcomp_geom_table + + def assign_nets( + self, + open_pins: Union[list, + None] = None) -> dict[Union[str, int], list[str]]: + """Assigns a netlist number to each galvanically connected metal region, + and returns a dictionary with each net as a key, and the corresponding list of + geometries associated with that net as values. + + Args: + open_pins (Union[list, None], optional): List of tuples of pins that are open. + Defaults to None. + + Returns: + dict[Union[str, int], list[str]]: dictionary with keys for each net, and list of + values with the corresponding geometries associated with that net as values. + """ + + netlists = dict() + netlist_id = 0 + + qgeom_names = self.qcomp_geom_table["name"] + qcomp_names_for_qgeom = [ + list(self.design.components.keys())[i - 1] + for i in self.qcomp_geom_table["component"] + ] + phys_grps = [ + s1 + '_' + s2 for s1, s2 in zip(qcomp_names_for_qgeom, qgeom_names) + ] + qgeom_idxs = list(range(len(self.qcomp_geom_table))) + id_net_dict = {k: -1 for k in phys_grps} + + while len(qgeom_idxs) != 0: + i = qgeom_idxs.pop(0) + shape_i = self.qcomp_geom_table.iloc[[i]]["geometry"][i] + chip_i = self.qcomp_geom_table.iloc[[i]]["chip"][i] + layer_i = self.qcomp_geom_table.iloc[[i]]["layer"][i] + thick_i, z_coord_i = self.gmsh.get_thickness_zcoord_for_layer_datatype( + layer_i) + id_net_dict[phys_grps[i]] = netlist_id if ( + id_net_dict[phys_grps[i]] == -1) else id_net_dict[phys_grps[i]] + for j in qgeom_idxs: + shape_j = self.qcomp_geom_table.iloc[[j]]["geometry"][j] + chip_j = self.qcomp_geom_table.iloc[[j]]["chip"][j] + layer_j = self.qcomp_geom_table.iloc[[j]]["layer"][j] + thick_j, z_coord_j = self.gmsh.get_thickness_zcoord_for_layer_datatype( + layer_j) + dist = shape_i.distance(shape_j) + + layers_touch = False + if (layer_i == layer_j or z_coord_j == z_coord_i or + z_coord_i + thick_i == z_coord_j or + z_coord_j + thick_j == z_coord_i): + layers_touch = True + + if dist == 0.0 and chip_i == chip_j and layers_touch: + if id_net_dict[phys_grps[j]] == -1: + id_net_dict[phys_grps[j]] = id_net_dict[phys_grps[i]] + elif id_net_dict[phys_grps[j]] != id_net_dict[phys_grps[i]]: + net_id_i = id_net_dict[phys_grps[i]] + for k, v in id_net_dict.items(): + if v == net_id_i: + id_net_dict[k] = id_net_dict[phys_grps[j]] + + if -1 not in id_net_dict.values(): + break + + netlist_id = max(list(id_net_dict.values())) + 1 + + gnd_phys_grps = self.get_gnd_qgeoms(open_pins) + gnd_netlist = list({ + net for (name, net) in id_net_dict.items() + if (name in gnd_phys_grps) + }) + + netlists['gnd'] = list() + for k, v in id_net_dict.items(): + if v in gnd_netlist: + netlists['gnd'].append(k) + else: + if v not in netlists.keys(): + netlists[v] = list() + netlists[v].append(k) + + netlists = {(i - 1 if (k != 'gnd') else k): v + for i, (k, v) in enumerate(netlists.items())} + + return netlists + + def get_gnd_qgeoms(self, open_pins: Union[list, None] = None) -> list[str]: + """ Obtain a list of qgeometry names associated with pins shorted to ground. + + Args: + open_pins (Union[list, None], optional): List of tuples of pins that are open. + Defaults to None. + + Returns: + list[str]: Names of qgeometry components with pins connected to ground plane. + """ + + open_pins = open_pins if open_pins is not None else [] + qcomp_lst = self.design.components.keys() + all_pins = list() + gnd_qgeoms = set() + + for qcomp in qcomp_lst: + qcomp_pins = self.design.components[qcomp].pins.keys() + all_pins += list(zip([qcomp] * len(qcomp_pins), qcomp_pins)) + + nets_table = self.design.net_info + gnd_nets = list(nets_table[nets_table["pin_name"] == "short"]["net_id"]) + gnd_nets_mask = nets_table["net_id"].isin(gnd_nets) + port_nets_table = nets_table[~gnd_nets_mask] + port_pins_names = port_nets_table["pin_name"] + qcomp_names_for_port_pins = [ + list(qcomp_lst)[i - 1] for i in port_nets_table["component_id"] + ] + port_pins = list(zip(qcomp_names_for_port_pins, port_pins_names)) + all_open_pins = list(set(port_pins + open_pins)) + gnd_pins = [pin for pin in all_pins if pin not in all_open_pins] + + for pin in gnd_pins: + qcomp_name, qcomp_pin = pin + pin_qcomp_id = self.design.components[qcomp_name].id + + if pin_qcomp_id in list(self.qcomp_geom_table["component"]): + pin_qcomp_geom_table = self.qcomp_geom_table[ + self.qcomp_geom_table["component"] == pin_qcomp_id] + pin_point = draw.Point( + self.design.components[qcomp_name].pins[qcomp_pin] + ['middle']) + + for _, qgeom_row in pin_qcomp_geom_table.iterrows(): + if pin_point.intersects(qgeom_row["geometry"]): + gnd_qgeoms.add(qcomp_name + '_' + qgeom_row["name"]) + + return list(gnd_qgeoms) + + def run(self, + sim_type: str, + display_cap_matrix: bool = False) -> Union[None, pd.DataFrame]: + """Runs ElmerFEM analysis. + + Args: + sim_type (str): Type of simulation to be executed by ElmerFEM. + display_cap_matrix (bool, optional): Whether or not to return cap matrix. + Defaults to False. + + Returns: + Union[None, pd.DataFrame]: extracted capacitance matrix if + display_cap_matrix = True. + """ + + setup = self.default_setup[sim_type] + sim_dir = self._options["simulation_dir"] + meshfile = self._options["mesh_file"] + sif_name = self._options["simulation_input_file"] + if sim_type == "capacitance": + cap_matrix_file = os.path.join(self._options["simulation_dir"], + setup["Capacitance_Matrix_Filename"]) + + self.logger.info("Running ElmerGrid on input mesh from Gmsh...") + self._elmer_runner.run_elmergrid(sim_dir, meshfile) + self.logger.info(f"Running ElmerSolver for solver type: '{sim_type}'") + self._elmer_runner.run_elmersolver(sim_dir, sif_name) + + self.capacitance_matrix = self._get_capacitance_matrix(cap_matrix_file) + if display_cap_matrix: + return self.capacitance_matrix + + def save_capacitance_matrix(self, path: str): + """Saves capacitance matrix to file. + + Args: + path (str): path where capacitance matrix should be saved. + """ + self.capacitance_matrix.to_csv(path, sep=' ', header=True) + + def _get_capacitance_matrix(self, filename: str) -> pd.DataFrame: + """Reads SPICE capacitance matrix file generated by ElmerFEM, converts + and returns Maxwell capacitance matrix + + Args: + filename (str): Name of capacitance matrix file generated by ElmerFEM + + Returns: + pd.DataFrame: Maxwell capacitance matrix + """ + with open(filename, "r", encoding="utf-8") as f: + lines = f.readlines() + f.close() + + new_lines = [(l.replace(" ", " ").strip() + '\n') for l in lines] + + with open(filename, "w", encoding="utf-8") as f: + f.writelines(new_lines) + f.close() + + df = pd.read_csv(filename, delimiter=" ", header=None) + + # Convert from SPICE cap matrix to Maxwell cap matrix + df2 = df.multiply(-1e15) + row_sum = df2.sum(axis=0) + for i, s in enumerate(row_sum): + df2[i][i] = -s + + name_cap = {k: v[-1] for k, v in self.nets.items() if k != "gnd"} + df2.rename(index=name_cap, columns=name_cap, inplace=True) + + gnd_caps = -1 * df2.sum(axis=0) + df2.loc["ground_plane"] = gnd_caps + df2["ground_plane"] = gnd_caps + df2 = df2.multiply( + self.default_setup["constants"]["Permittivity_of_Vacuum"]) + + # dummy value set to high as it's shorted to ground + df2["ground_plane"]["ground_plane"] = 300 + return df2 + + def add_solution_setup( + self, + sim_name: str = "capacitance", + solver_names: list[str] = ["capacitance", "postprocessing_gmsh"], + equation_name: str = "poisson", + ): + """Initializes ElmerFEM analysis to run. + Currently only supports Electrostatic capacitance extraction + + Args: + sim_name (str, optional): Type of ElmerFEM analysis to initialize. + Defaults to "capacitance". + solver_names (list[str], optional): ElmerFEM solver to use. Defaults + to ["capacitance", "postprocessing_gmsh"]. + equation_name (str, optional): Type of equation for solver. Defaults to "poisson". + """ + setup = self.default_setup + cap_setup = { + k.replace('_', ' '): v + for k, v in self.default_setup[sim_name].items() + } + + # Load simulation + self._elmer_runner.load_simulation(self._options["simulation_type"]) + + # Load constants + constants = self.default_setup["constants"] + for k in constants.keys(): + const_str = k.replace('_', ' ') + self._elmer_runner.add_constant(const_str, 1) + + # Load solvers + solvers = self._elmer_runner.load_solvers(solver_names=solver_names) + self._elmer_runner.solver["capacitance"].update(cap_setup) + + # Load materials + materials = self._elmer_runner.load_materials(setup["materials"]) + + # Add equation + eqn = self._elmer_runner.add_equation(equation_name, solvers=solvers) + + # Add bodies + bodies = self.define_bodies(setup, eqn, materials) + + # Add boundary conditions + boundaries, cap_body = self.define_boundaries() + + if cap_body == 1: + self.logger.warning( + "WARNING: No capacitance bodies added to the model other " + "than those connected to ground. ElmerFEM cannot run a " + "capacitance extraction analysis without any bodies. ") + + # Update capacitance bodies in solver + self._elmer_runner.solver["capacitance"].update( + {"Capacitance Bodies": cap_body - 1}) + + # Write the simulation input file (SIF) + self.write_sif() + + def define_bodies(self, setup: dict, equation: int, + materials: list[int]) -> list[list]: + """Assigns bodies to simulation model. + + Args: + setup (dict): Setup parameters to be used in simulation. + equation (int): Type of equation for solver. + materials (list[int]): Materials used by bodies in the design. + + Returns: + list[list]: Information about the bodies in the design. + """ + bodies = [] + for i, material in enumerate(setup["materials"]): + if material == "vacuum": + bodies += [ + self._elmer_runner.add_body( + "vacuum_box", + [self.gmsh.physical_groups["global"]["vacuum_box"]], + equation, materials[i]) + ] + if material == "silicon": + substrate_geoms = dict() + for layer, ph_geoms in self.gmsh.physical_groups.items(): + if layer in ["global", "chips"]: + continue + + sub = { + k: v + for k, v in ph_geoms.items() + if ("dielectric" in k) and ("_sfs" not in k) + } + substrate_geoms.update(sub) + + bodies += [ + self._elmer_runner.add_body(name, [body], equation, + materials[i]) + for name, body in substrate_geoms.items() + ] + + return bodies + + def define_boundaries(self) -> tuple[list[list], int]: + """Assigns boundaries to simulation model. + + Returns: + tuple[list[list], int]: Information about the bodies in the design, + and total number of capacitance bodies. + """ + + boundaries = list() + phys_grps_dict = dict() + ground_planes = list() + cap_body = 1 + for layer, ph_geoms in self.gmsh.physical_groups.items(): + if layer in ["global", "chips"]: + continue + + for name, geom in ph_geoms.items(): + if ("ground_plane" in name) and ("sfs" in name): + ground_planes.append(geom) + + if "dielectric" in name: + continue + + phys_grps_dict.update(ph_geoms) + + for net, geom_names in self.nets.items(): + geom_id_list = [ + phys_grps_dict[f"{name}_sfs"] for name in geom_names + ] + + if net == "gnd": + ground_planes += geom_id_list + else: + boundaries += [ + self._elmer_runner.add_boundary_condition( + geom_names[-1], + geom_id_list, + options={"Capacitance Body": cap_body}) + ] + cap_body += 1 + + self._elmer_runner.add_boundary_condition( + "ground_plane", ground_planes, options={"Capacitance Body": 0}) + + self._elmer_runner.add_boundary_condition( + "FarField", [self.gmsh.physical_groups["global"]["vacuum_box_sfs"]], + options={"Electric Infinity BC": True}) + + return boundaries, cap_body + + def write_sif(self): + """Write sif (simulation input file) used by ElmerFEM to run analysis""" + filename = self._options["simulation_input_file"] + sim_dir = self._options["simulation_dir"] + self._elmer_runner.write_startinfo_and_sif(filename=filename, + sim_dir=sim_dir) + + def render_chips(self, + chips: Union[str, list[str]] = [], + draw_sample_holder: bool = True, + box_plus_buffer: bool = True): + """Abstract method. Must be implemented by the subclass. + Render all chips of the design. + Calls render_chip for each chip. + """ + pass + + def render_chip(self, chip_name: str): + """Abstract method. Must be implemented by the subclass. + Render the given chip. + + Args: + name (str): chip to render + """ + pass + + def render_layers(self, + draw_sample_holder: bool = True, + layers: Union[list[int], None] = None, + box_plus_buffer: bool = True): + """Abstract method. Must be implemented by the subclass. + Render all layers of the design. + Calls render_layer for each layer. + """ + self.gmsh.render_layers(draw_sample_holder, layers, box_plus_buffer) + + def render_layer(self, layer_number: str, datatype: int = 0): + """Abstract method. Must be implemented by the subclass. + Render the given layer. + + Args: + name (str): layer to render + """ + self.gmsh.render_layer(layer_number, datatype) + + def render_components(self, table_type: str): + """Abstract method. Must be implemented by the subclass. + Render all components of the design. + If selection is none, then render all components. + + Args: + selection (QComponent): Component to render. + """ + self.gmsh.render_components(table_type) + + def render_component(self, component): + """Abstract method. Must be implemented by the subclass. + Render the specified component. + + Args: + component (QComponent): Component to render. + """ + self.gmsh.render_component(component) + + def render_element(self, qgeom: pd.Series, table_type: str): + """Abstract method. Must be implemented by the subclass. + Render the specified element + + Args: + element (Element): Element to render. + """ + self.gmsh.render_element(qgeom, table_type) + + def render_element_path(self, path: pd.Series): + """Abstract method. Must be implemented by the subclass. + Render an element path. + + Args: + path (str): Path to render. + """ + self.gmsh.render_element_path(path) + + def render_element_junction(self, junc: pd.Series): + """Abstract method. Must be implemented by the subclass. + Render an element junction. + + Args: + junc (str): Junction to render. + """ + self.gmsh.render_element_junction(junc) + + def render_element_poly(self, poly: pd.Series): + """Abstract method. Must be implemented by the subclass. + Render an element poly. + + Args: + poly (Poly): Poly to render. + """ + self.gmsh.render_element_poly(poly) + + def save_screenshot(self, path: str = None, show: bool = True): + """Save the screenshot. + + Args: + path (str, optional): Path to save location. Defaults to None. + show (bool, optional): Whether or not to display the screenshot. Defaults to True. + + Returns: + pathlib.WindowsPath: path to png formatted screenshot. + """ + self.gmsh.save_screenshot(path, show) + + def launch_gmsh_gui(self): + """Launch Gmsh GUI for viewing the model. + """ + self.gmsh.launch_gui() + + def export_mesh(self): + """Export Gmsh mesh + """ + self.gmsh.export_mesh(self._options["mesh_file"]) + + def display_post_processing_data(self): + """Import data given by ElmerFEM for Post-Processing in Gmsh + """ + postprocessing_file = os.path.join(self._options["simulation_dir"], + self._options["postprocessing_file"]) + self.gmsh.import_post_processing_data(postprocessing_file) diff --git a/qiskit_metal/renderers/renderer_elmer/elmer_runner.py b/qiskit_metal/renderers/renderer_elmer/elmer_runner.py new file mode 100644 index 000000000..9dcce2b6b --- /dev/null +++ b/qiskit_metal/renderers/renderer_elmer/elmer_runner.py @@ -0,0 +1,641 @@ +from dataclasses import dataclass +from typing import Union, Any +import os +import subprocess +import shutil +from collections import defaultdict +import platform +import yaml + +from qiskit_metal import Dict +from qiskit_metal.renderers.renderer_elmer.elmer_configs import materials, simulations, solvers + + +@dataclass +class Body: + """ Class for storing data for ElmerFEM Body """ + target_bodies: list[int] + equation: int + material: int + body_force: Union[int, None] + initial_condition: Union[int, None] + + +@dataclass +class Boundary: + """ Class for storing data for ElmerFEM boundary """ + target_boundaries: list[int] + other_opts: dict + + +class ElmerRunner: + """Class for generating the Simulation Input File (SIF) for ElmerFEM to run the simulation. + + Attributes: + header: Header information for the SIF + simulation: Simulation information for the SIF + constant: Constants for the SIF + solver: Solver information for the SIF + equation: Equation information for the SIF + material: Material information for the SIF + body_force: Body force information for the SIF + initial_condition: Initial condition information for the SIF + body: Body information for the SIF + boundary_condition: Boundary condition information for the SIF + """ + + def __init__(self) -> None: + """Initialize the ElmerRunner object + + Args: + self: reference to self object + """ + self.header = 'Header\n CHECK KEYWORDS "Warn"\n Mesh DB "." "."\nEnd\n\n' + self.simulation = defaultdict(dict) + self.constant = dict() + self.solver = defaultdict(dict) + self.equation = defaultdict(list) + self.material = defaultdict(dict) + self.body_force = defaultdict(dict) + self.initial_condition = defaultdict(dict) + self.body = defaultdict(Body) + self.boundary_condition = defaultdict(Boundary) + + @property + def solver_info(self) -> list[str]: + """Display the types of solvers used in the simulation. + + Returns: + list[str]: returns the list of solver names + """ + return list(self.solver.keys()) + + @property + def simulation_info(self) -> list[str]: + """Display the types of simulation used in the SIF. + + Returns: + list[str]: returns the list of simulation names + """ + return list(self.simulation.keys()) + + @property + def equation_info(self) -> list[str]: + """Display the types of equations used in the simulation. + + Returns: + list[str]: returns the list of equation names + """ + return list(self.equation.keys()) + + def add_setup_options( + self, + setup_param: Dict, + new_opt: dict[str, Any], + ): + """Function to modify the setup options using the ones + provided by the user in QElmerRenderer. + + Args: + setup_param (Dict): Setup parameters to be written in SIF. + new_opt (dict[str, Any]): user-modified parameters from the + QElmerRenderer setup. + """ + modified_new_opt = {k.replace(' ', '_'): v for k, v in new_opt.items()} + setup_param.update(modified_new_opt) + + def init_simulation(self, name: str, options: dict): + """Function to initialize a simulation. + + Args: + name (str): Simulation name. + options (dict): Simulation options. + """ + if len(self.simulation.keys()) > 1: + raise RuntimeError( + "There can only be one simulation instance at a time.") + + self.simulation.update({name: options}) + + def add_equation(self, name: str, solvers: list) -> int: + """Function to add an equation to the current simulation. + + Args: + name (str): Equation name. + solvers (list): Solvers corresponding to each equation. + + Returns: + int: Equation ID in the simulation. + """ + self.equation.update({name: solvers}) + return len(self.equation) + + def add_constant(self, name: str, value: str): + """Function to add a constant to the simulation. + + Args: + name (str): Constant name. + value (str): Value of the constant in appropriate units. + """ + self.constant.update({name: value}) + + def add_solver(self, name: str, options: dict) -> int: + """Function to add a solver to the simulation. + + Args: + name (str): Solver name. + options (dict): Solver options. + + Returns: + int: Solver ID in the simulation. + """ + self.solver.update({name: options}) + return len(self.solver) + + def add_material(self, name: str, options: dict) -> int: + """Function to add a material in the simulation. + + Args: + name (str): Material name. + options (dict): Material options. + + Returns: + int: Material ID in the simulation. + """ + self.material.update({name: options}) + return len(self.material) + + def add_body(self, + name: str, + target_bodies: list[int], + equation: int, + material: int, + body_force: Union[int, None] = None, + initial_condition: Union[int, None] = None) -> int: + """Function to add a Body in the simulation. + + Args: + name (str): Body name. + target_bodies (list[int]): Target ID of volume in the mesh. + equation (int): Equation associated with the body. + material (int): Material associated with the body. + body_force (Union[int, None]): Body force associated with the body. + initial_condition (Union[int, None]): Initial condition associated with the body. + + Returns: + int: Body ID in the simulation. + """ + sim_body = Body(target_bodies, equation, material, body_force, + initial_condition) + self.body.update({name: sim_body}) + return len(self.body) + + def add_body_force(self, name: str, options: dict) -> int: + """Function to add a body force in the simulation. + + Args: + name (str): Body force name. + options (dict): Body force options. + + Returns: + int: Body force ID in the simulation. + """ + self.body_force.update({name: options}) + return len(self.material) + + def add_boundary_condition(self, name: str, target_boundaries: list[int], + options: dict) -> int: + """Function to add a boundary condition in the simulation. + + Args: + name (str): Boundary condition name. + target_boundaries (list[int]): Target ID of surface/curve in the mesh. + options (dict): Boundary condition options. + + Returns: + int: Boundary condition ID in the simulation. + """ + sim_boundary = Boundary(target_boundaries, options) + self.boundary_condition.update({name: sim_boundary}) + return len(self.boundary_condition) + + def add_initial_condition(self, name: str, options: dict) -> int: + """Function to add an initial condition in the simulation. + + Args: + name (str): Initial condition name. + options (dict): Initial condition options. + + Returns: + int: Initial condition ID in the simulation. + """ + self.initial_condition.update({name: options}) + return len(self.initial_condition) + + def _update_solver_setup(self, yaml_setup_dict: dict, + new_setup_params: dict) -> dict: + """Function to update the solver setup from QElmerRenderer + default setup. + + Args: + yaml_setup_dict (dict): setup dictionary from the YAML file. + new_setup_params (dict): User-modified setup parameters + in QElmerRenderer. + + Returns: + dict: Modifed YAML setup dictionary + """ + yaml_setup_dict.update(new_setup_params) + return yaml_setup_dict + + def _read_yaml_to_dict(self, filename: str) -> dict: + """Function to read a YAML setup file for simulation ad solver details. + + Args: + filename (str): file for the YAML setup. + + Returns: + dict: Returns dictionary populated using the YAML file. + """ + with open(filename, "r", encoding="utf-8") as f: + yaml_setup = yaml.safe_load(f) + f.close() + return yaml_setup + + def load_simulation(self, name: str, config_file: str = None): + """Function to load the simulation. + + Args: + name (str): Simulation name. + config_file (str): Optional YAML configuration file + to load the simulation from. + Defaults to None. + """ + if config_file is not None: + opts = self._read_yaml_to_dict(config_file) + else: + opts = simulations + if name not in opts.keys(): + raise ValueError( + "Simulation not defined in the provided config file. " + "Either add it to your config file, or use `init_simulation` " + "to manually add new simulation to the configuration.") + self.init_simulation(name=name, options=opts[name]) + + def load_solvers(self, + solver_names: list[str], + config_file: str = None) -> list[int]: + """Function to load the solvers (from YAML configutation if specidifed). + + Args: + solver_names (list[str]): Solver names. + config_file (str): Optional YAML configuration file + to load the simulation from. + Defaults to None. + + Returns: + list[int]: Returns a list of solvers that are loaded and active. + """ + if config_file is not None: + opts = self._read_yaml_to_dict(config_file) + else: + opts = solvers + for name in solver_names: + if name not in opts.keys(): + raise ValueError( + "Solver not defined in the provided config file. " + "Either add it to your config file, or use `add_solver` " + "to manually add new solver to the simulation.") + self.add_solver(name=name, options=opts[name]) + + return [i + 1 for i in range(len(self.solver))] + + def load_materials(self, + material_names: list[str], + config_file: str = None) -> list[int]: + """Function to load the materials (from YAML configuration if specified). + + Args: + material_names (list[str]): Material names. + config_file (str): Optional YAML configuration file + to load the simulation from. + Defaults to None. + + Returns: + list[int]: Returns a list of materials that are loaded and active. + """ + if config_file is not None: + opts = self._read_yaml_to_dict(config_file) + else: + opts = materials + for name in material_names: + if name not in opts.keys(): + raise ValueError( + "Material not defined in the provided config file. " + "Either add it to your config file, or use `add_material` " + "to manually add new material to the simulation.") + self.add_material(name=name, options=opts[name]) + + return [i + 1 for i in range(len(self.material))] + + def _write_sif_simulation(self) -> str: + """Write the simulation in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.simulation) > 0: + section_string += f"! {list(self.simulation.keys())[0]}\n" + section_string += "Simulation\n" + for k, v in list(self.simulation.values())[0].items(): + section_string += f" {k} = {v}\n" + section_string += "End\n\n" + return section_string + + def _write_sif_constants(self) -> str: + """Write the constants in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.constant) > 0: + section_string += "Constants\n" + for k, v in self.constant.items(): + section_string += f" {k} = {v}\n" + section_string += "End\n\n" + return section_string + + def _write_sif_equation(self) -> str: + """Write the equations in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.equation) > 0: + for i, k in enumerate(self.equation): + section_string += f"! {k}\n" + section_string += f"Equation {i+1}\n" + solvers = [ + list(self.solver.keys())[i - 1] for i in self.equation[k] + ] + active_solvers = "" + for i, solver in enumerate(solvers): + section_string += f" ! {i+1}: {solver}\n" + active_solvers += f" {i+1}" + + section_string += f" Active Solvers({len(solvers)}) =" + active_solvers + '\n' + section_string += "End\n\n" + section_string += '\n' + return section_string + + def _write_sif_solver(self) -> str: + """Write the solvers in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.solver) > 0: + for i, k in enumerate(self.solver): + section_string += f"! {k}\n" + section_string += f"Solver {i+1}\n" + for key, val in self.solver[k].items(): + section_string += f" {key} = {val}\n" + section_string += "End\n\n" + section_string += '\n' + return section_string + + def _write_sif_material(self) -> str: + """Write the materials in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.material) > 0: + for i, k in enumerate(self.material): + section_string += f"! {k}\n" + section_string += f"Material {i+1}\n" + for key, val in self.material[k].items(): + section_string += f" {key} = {val}\n" + section_string += "End\n\n" + section_string += '\n' + return section_string + + def _write_sif_body(self) -> str: + """Write the bodies in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.body) > 0: + for i, k in enumerate(self.body): + section_string += f"! {k}\n" + section_string += f"Body {i+1}\n" + body_obj = self.body[k] + trgt_bds = "" + for tb in body_obj.target_bodies: + trgt_bds += f" {tb}" + section_string += f" Target Bodies({len(body_obj.target_bodies)}) =" + trgt_bds + '\n' + section_string += f" Equation = {body_obj.equation} ! {list(self.equation.items())[body_obj.equation-1][0]}\n" + section_string += f" Material = {body_obj.material} ! {list(self.material.items())[body_obj.material-1][0]}\n" + if body_obj.body_force is not None: + section_string += ( + f" Body Force = {body_obj.body_force}" + f" ! {list(self.body_force.items())[body_obj.body_force-1][0]}\n" + ) + if body_obj.initial_condition is not None: + section_string += ( + f" Initial Condition = {body_obj.initial_condition}" + f" ! {list(self.initial_condition.items())[body_obj.initial_condition-1][0]}\n" + ) + section_string += "End\n\n" + section_string += '\n' + return section_string + + def _write_sif_boundary_condition(self) -> str: + """Write the boundary conditions in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.boundary_condition) > 0: + for i, k in enumerate(self.boundary_condition): + section_string += f"! {k}\n" + section_string += f"Boundary Condition {i+1}\n" + bdry_obj = self.boundary_condition[k] + trgt_bds = "" + for tb in bdry_obj.target_boundaries: + trgt_bds += f" {tb}" + section_string += f" Target Boundaries({len(bdry_obj.target_boundaries)}) =" + trgt_bds + '\n' + for k, v in bdry_obj.other_opts.items(): + section_string += f" {k} = {v}\n" + section_string += "End\n\n" + section_string += '\n' + return section_string + + def _write_sif_initial_condition(self) -> str: + """Write the initial conditions in SIF using data stored + as a single Python string. + + Returns: + str: Returns the partially written string for SIF. + """ + section_string = "" + if len(self.initial_condition) > 0: + for i, k in enumerate(self.initial_condition): + section_string += f"! {k}\n" + section_string += f"Initial Condition {i+1}\n" + for key, val in self.initial_condition[k].items(): + section_string += f" {key} = {val}\n" + section_string += "End\n\n" + section_string += '\n' + return section_string + + def write_startinfo_and_sif(self, filename: str, sim_dir: str): + """Write the SIF and STARTINFO files for ElmerFEM to + provide the simulation setup. + + Args: + filename (str): filename for the SIF. + sim_dir (str): Directory for storing all the simulation data. + """ + if not os.path.exists(sim_dir): + os.mkdir(sim_dir) + + startinfo_path = os.path.join(sim_dir, "ELMERSOLVER_STARTINFO") + with open(startinfo_path, "w+", encoding="utf-8") as startinfo: + startinfo.write(filename) + startinfo.close() + + file_path = os.path.join(sim_dir, filename) + with open(file_path, "w+", encoding="utf-8") as sif: + sif_data = self.header + sif_data += self._write_sif_simulation() + sif_data += self._write_sif_constants() + sif_data += self._write_sif_equation() + sif_data += self._write_sif_solver() + sif_data += self._write_sif_material() + sif_data += self._write_sif_body() + sif_data += self._write_sif_boundary_condition() + sif_data += self._write_sif_initial_condition() + + sif.write(sif_data) + sif.close() + + def run_elmergrid(self, + sim_dir: str, + meshfile: str, + elmergrid: str = None, + options: list = []): + """Function to run ElmerGrid executable in the system for transforming + Gmsh style mesh to ElmerFEM style mesh. + + Args: + sim_dir (str): Directory for storing simulation data. + meshfile (str): Gmsh mesh file + elmergrid (str): Path to ElmerGrid executable (if not installed + using the default procedure). Defaults to None. + options (list): ElmerGrid additional options. + """ + os_platform = platform.system() + if elmergrid is None: + # On Windows ElmerGrid.exe is not found once gmsh.initialize() was executed. + # Try to use abs-path instead. + if os_platform == 'Windows' and os.path.exists( + "C:/Program Files/Elmer 9.0-Release/bin/ElmerGrid.exe"): + elmergrid = "C:/Program Files/Elmer 9.0-Release/bin/ElmerGrid.exe" + else: + elmergrid = "ElmerGrid" + + args = [elmergrid, "14", "2", os.path.join("..", meshfile)] + options + with open(os.path.join(sim_dir, "elmergrid.log"), + "w+", + encoding="utf-8") as f: + subprocess.run(args, cwd=sim_dir, stdout=f, stderr=f) + + mesh_dir = meshfile.split(".")[-2] + files = os.listdir(mesh_dir) + + if not os.path.exists(sim_dir): + os.mkdir(sim_dir) + + for f in files: + if os.path.exists(os.path.join(sim_dir, f)): + os.remove(os.path.join(sim_dir, f)) + shutil.move(os.path.join(mesh_dir, f), sim_dir) + shutil.rmtree(mesh_dir) + + def run_elmersolver( + self, + sim_dir: str, + sif_file: str, + # out_files: ty.List[str], + elmersolver: str = None, + options: list = []): + """Function to run ElmerSolver executable in the system + running the FEM simulation. + + Args: + sim_dir (str): Directory for storing simulation data. + sif_file (str): Simulation Input File (SIF) path + elmersolver (str): Path to ElmerSolver executable (if not installed + using the default procedure). Defaults to None. + options (list): ElmerSolver additional options. + """ + os_platform = platform.system() + if elmersolver is None: + # On Windows ElmerSolver.exe is not found once gmsh.initialize() was executed. + # Try to use abs-path instead. + if os_platform == 'Windows' and os.path.exists( + "C:/Program Files/Elmer 9.0-Release/bin/ElmerSolver.exe"): + elmersolver = "C:/Program Files/Elmer 9.0-Release/bin/ElmerSolver.exe" + else: + elmersolver = "ElmerSolver" + + args = [elmersolver, sif_file] + options + with open(os.path.join(sim_dir, "elmersolver.log"), + "w+", + encoding="utf=8") as f: + subprocess.run(args, cwd=sim_dir, stdout=f, stderr=f) + + # for f in out_files: + # out_file = os.path.join(sim_dir, f) + # if os.path.exists(os.path.join(".", f)): + # os.remove(os.path.join(".", f)) + # shutil.move(out_file, ".") + + # def run_elmersolver_mpi(self, + # sim_dir: str, + # sif_file: str, + # elmersolver: str = None, + # num_procs: int = 8, + # options: list = []): + + # os_platform = platform.system() + # if elmersolver is None: + # # On Windows ElmerSolver.exe is not found once gmsh.initialize() was executed. + # # Try to use abs-path instead. + # if os_platform == 'Windows' and os.path.exists( + # "C:/Program Files/Elmer 9.0-Release/bin/ElmerSolver_mpi.exe"): + # elmersolver = "C:/Program Files/Elmer 9.0-Release/bin/ElmerSolver_mpi.exe" + # else: + # elmersolver = "ElmerSolver_mpi" + + # # TODO: how to specify number of processes? + # args = [elmersolver, sif_file] + options + # with open(os.path.join(sim_dir, "elmersolver.log", "w+"), + # encoding="utf=8") as f: + # subprocess.run(args, cwd=sim_dir, stdout=f, stderr=f) diff --git a/qiskit_metal/renderers/renderer_gmsh/gmsh_renderer.py b/qiskit_metal/renderers/renderer_gmsh/gmsh_renderer.py index 0f0efc924..cdca676cf 100644 --- a/qiskit_metal/renderers/renderer_gmsh/gmsh_renderer.py +++ b/qiskit_metal/renderers/renderer_gmsh/gmsh_renderer.py @@ -22,7 +22,6 @@ class QGmshRenderer(QRenderer): """Extends QRendererAnalysis class to export designs to Gmsh using the Gmsh python API. Default Options: - # Buffer between max/min x and edge of ground plane, in mm * x_buffer_width_mm -- Buffer between max/min x and edge of ground plane, in mm * y_buffer_width_mm -- Buffer between max/min y and edge of ground plane, in mm * mesh -- to define meshing parameters @@ -96,11 +95,9 @@ def __init__(self, @property def initialized(self): - """Abstract method. Must be implemented by the subclass. - Is renderer ready to be used? - Implementation must return boolean True if successful. False otherwise. - """ - if self.model == self._model_name: + """Returns boolean True if initialized successfully. + False otherwise.""" + if gmsh.isInitialized(): return True return False @@ -214,6 +211,7 @@ def render_design( skip_junctions: bool = False, mesh_geoms: bool = True, ignore_metal_volume: bool = False, + omit_ground_for_layers: Optional[list[int]] = None, ): """Render the design in Gmsh and apply changes to modify the geometries according to the type of simulation. Simulation parameters provided by the user. @@ -233,8 +231,15 @@ def render_design( ignore_metal_volume (bool, optional): ignore the volume of metals and replace it with a list of surfaces instead. Defaults to False. + omit_ground_for_layers (Optional[list[int]]): omit rendering the ground plane for + specified layers. Defaults to None. """ + # For handling the case when the user wants to use + # QGmshRenderer from design.renderers.gmsh instance. + if not self.initialized: + self._initiate_renderer() + # defaultdict: chip -- geom_tag self.layers_dict = defaultdict(list) @@ -253,7 +258,8 @@ def render_design( open_pins=open_pins, box_plus_buffer=box_plus_buffer, draw_sample_holder=draw_sample_holder, - skip_junctions=skip_junctions) + skip_junctions=skip_junctions, + omit_ground_for_layers=omit_ground_for_layers) self.apply_changes_for_simulation( ignore_metal_volume=ignore_metal_volume, @@ -270,7 +276,8 @@ def draw_geometries(self, selection: Union[list, None] = None, open_pins: Union[list, None] = None, box_plus_buffer: bool = True, - skip_junctions: bool = False): + skip_junctions: bool = False, + omit_ground_for_layers: Optional[list[int]] = None): """This function draws the raw geometries in Gmsh as taken from the QGeometry tables and applies thickness depending on the layer-stack. @@ -284,6 +291,8 @@ def draw_geometries(self, draw_sample_holder (bool): To draw the sample holder box. skip_junctions (bool, optional): Set to True to sip rendering the junctions. Defaults to False. + omit_ground_for_layers (Optional[list[int]]): omit rendering the ground plane for + specified layers. Defaults to None. """ self.qcomp_ids, self.case = self.get_unique_component_ids(selection) @@ -296,8 +305,9 @@ def draw_geometries(self, self.render_tables(skip_junction=skip_junctions) self.add_endcaps(open_pins=open_pins) self.render_layers(box_plus_buffer=box_plus_buffer, + omit_layers=omit_ground_for_layers, draw_sample_holder=draw_sample_holder) - self.subtract_from_layers() + self.subtract_from_layers(omit_layers=omit_ground_for_layers) self.gmsh_occ_synchronize() def apply_changes_for_simulation(self, ignore_metal_volume: bool, @@ -602,8 +612,17 @@ def add_endcaps(self, open_pins: Union[list, None] = None): open_pins = open_pins if open_pins is not None else [] for comp, pin in open_pins: + if comp not in self.design.components: + raise ValueError( + f"Component '{comp}' not present in current design.") + qcomp = self.design.components[comp] qc_layer = int(qcomp.options.layer) + + if pin not in qcomp.pins: + raise ValueError( + f"Pin '{pin}' not present in component '{comp}'.") + pin_dict = qcomp.pins[pin] width, gap = self.parse_units_gmsh( [pin_dict["width"], pin_dict["gap"]]) @@ -647,18 +666,22 @@ def add_endcaps(self, open_pins: Union[list, None] = None): def render_layers(self, draw_sample_holder: bool, - layers: Optional[List[int]] = None, + omit_layers: Optional[List[int]] = None, box_plus_buffer: bool = True): """Render all chips of the design. calls `render_chip` to render the actual geometries Args: - layers (Optional[List[int]]): List of layers to render. - Renders all if [] or None is given. Defaults to None. + omit_layers (Optional[List[int]]): List of layers to omit render. + Renders all if [] or None is given. + Defaults to None. draw_sample_holder (bool): To draw the sample holder box. - box_plus_buffer (bool, optional): For adding buffer to chip dimensions. Defaults to True. + box_plus_buffer (bool, optional): For adding buffer to chip dimensions. + Defaults to True. """ - layer_list = list(set(l for l in self.design.ls.ls_df["layer"]) - ) if layers is None else layers + layer_list = list(set(l for l in self.design.ls.ls_df["layer"])) + + if omit_layers is not None: + layer_list = list(l for l in layer_list if l not in omit_layers) for layer in layer_list: # Add the buffer, using options for renderer. @@ -731,9 +754,18 @@ def render_layer(self, layer_number: int, datatype: int = 0): self.layers_dict[layer_number] = [layer_tag] - def subtract_from_layers(self): - """Subtract the QGeometries in tables from the chip ground plane""" + def subtract_from_layers(self, omit_layers: Optional[list[int]] = None): + """Subtract the QGeometries in tables from the chip ground plane + + Args: + omit_layers (Optional[List[int]]): List of layers to omit render. + Renders all if [] or None is given. + Defaults to None. + """ for layer_num, shapes in self.layer_subtract_dict.items(): + if omit_layers is not None and layer_num in omit_layers: + continue + thickness = self.get_thickness_for_layer_datatype( layer_num=layer_num) @@ -742,16 +774,18 @@ def subtract_from_layers(self): shape_dim_tags = [(dim, s) for s in shapes] layer_dim_tag = (dim, self.layers_dict[layer_num][0]) tool_dimtags = [layer_dim_tag] - subtract_layer = gmsh.model.occ.cut(tool_dimtags, shape_dim_tags) + if len(shape_dim_tags) > 0: + subtract_layer = gmsh.model.occ.cut(tool_dimtags, + shape_dim_tags) - updated_layer_geoms = [] - for i in range(len(tool_dimtags)): - if len(subtract_layer[1][i]) > 0: - updated_layer_geoms += [ - tag for _, tag in subtract_layer[1][i] - ] + updated_layer_geoms = [] + for i in range(len(tool_dimtags)): + if len(subtract_layer[1][i]) > 0: + updated_layer_geoms += [ + tag for _, tag in subtract_layer[1][i] + ] - self.layers_dict[layer_num] = updated_layer_geoms + self.layers_dict[layer_num] = updated_layer_geoms def fragment_interfaces(self, draw_sample_holder: bool): """Fragment Gmsh surfaces to ensure consistent tetrahedral meshing @@ -797,6 +831,7 @@ def fragment_interfaces(self, draw_sample_holder: bool): all_geom_dimtags) # Extract the new vacuum_box volume self.vacuum_box = fragmented_geoms[1][0][0][1] + object_dimtag = (3, self.vacuum_box) else: # Get one of the dim=3 objects dim3_dimtag = [ @@ -841,6 +876,33 @@ def fragment_interfaces(self, draw_sample_holder: bool): # junc_dimtags = [(2, junc) for junc in all_juncs] # gmsh.model.occ.fragment([(3, self.vacuum_box)], junc_dimtags) + def get_all_metal_surfaces(self): + metal_geoms = list() + surf_tags = list() + for layer in self.layer_types["metal"]: + thickness = self.get_thickness_for_layer_datatype(layer_num=layer) + + for _, tag in self.juncs_dict[layer].items(): + surf_tags += tag + + if thickness > 0.0: + for _, tag in self.polys_dict[layer].items(): + metal_geoms += tag + for _, tag in self.paths_dict[layer].items(): + metal_geoms += tag + metal_geoms += self.layers_dict[layer] + else: + for _, tag in self.polys_dict[layer].items(): + surf_tags += tag + for _, tag in self.paths_dict[layer].items(): + surf_tags += tag + surf_tags += self.layers_dict[layer] + + for geom in metal_geoms: + surf_tags += list(gmsh.model.occ.getSurfaceLoops(geom)[1][0]) + + return surf_tags + def assign_physical_groups(self, ignore_metal_volume: bool, draw_sample_holder: bool): """Assign physical groups to classify different geometries physically. @@ -854,7 +916,7 @@ def assign_physical_groups(self, ignore_metal_volume: bool, ValueError: if self.layer_types is not a dict ValueError: if layer number is not in self.layer_types """ - layer_numbers = list(self.layers_dict.keys()) + layer_numbers = list(set(l for l in self.design.ls.ls_df["layer"])) for layer in layer_numbers: # TODO: check if thickness == 0, then fragment differently layer_thickness = self.get_thickness_for_layer_datatype( @@ -915,13 +977,22 @@ def assign_physical_groups(self, ignore_metal_volume: bool, layer_name = layer_type + f'_(layer {layer})' layer_tag = self.layers_dict[layer] + all_metal_surfs = self.get_all_metal_surfaces() if len(layer_tag) > 0: if layer_dim == 3: layer_sfs_tags = [] for vol in layer_tag: - layer_sfs_tags += list( + layer_sfs = list( gmsh.model.occ.getSurfaceLoops(vol)[1][0]) + if layer_type == "ground_plane": + layer_sfs_tags += layer_sfs + else: + layer_sfs_tags += [ + sf for sf in layer_sfs + if sf not in all_metal_surfs + ] + if layer_type != "ground_plane" or ( layer_type == "ground_plane" and not ignore_metal_volume): @@ -1279,12 +1350,9 @@ def import_post_processing_data(self, raise ValueError( "Only .msh files supported for post processing views.") - if self.initialized: - self.model = str(self.model) + "_post_processing" - else: - self.model = "post_processing" + self.model = "post_processing" - gmsh.merge(filename) + gmsh.open(filename) if launch_gui: self.launch_gui() diff --git a/qiskit_metal/tests/test_renderers.py b/qiskit_metal/tests/test_renderers.py index 8bfd5f47b..8b3daa4d8 100644 --- a/qiskit_metal/tests/test_renderers.py +++ b/qiskit_metal/tests/test_renderers.py @@ -33,6 +33,7 @@ from qiskit_metal.renderers.renderer_gds.gds_renderer import QGDSRenderer from qiskit_metal.renderers.renderer_mpl.mpl_interaction import MplInteraction from qiskit_metal.renderers.renderer_gmsh.gmsh_renderer import QGmshRenderer +from qiskit_metal.renderers.renderer_elmer.elmer_renderer import QElmerRenderer from qiskit_metal.renderers.renderer_ansys import ansys_renderer @@ -131,6 +132,27 @@ def test_renderer_instantiate_qgmsh_renderer(self): "QGmshRenderer(design, initiate=False, options={}) failed in DesignPLanar" ) + def test_renderer_instantiate_qelmer_renderer(self): + """Test instantiation of QElmerRenderer in elmer_renderer.py.""" + design = designs.MultiPlanar() + try: + QElmerRenderer(design) + except Exception: + self.fail("QElmerRenderer(design) failed") + + try: + QElmerRenderer(design, initiate=False) + except Exception: + self.fail( + "QElmerRenderer(design, initiate=False) failed in MultiPlanar") + + try: + QElmerRenderer(design, initiate=False, options={}) + except Exception: + self.fail( + "QElmerRenderer(design, initiate=False, options={}) failed in MultiPlanar" + ) + def test_renderer_qansys_renderer_options(self): """Test that defaults in QAnsysRenderer were not accidentally changed.""" design = designs.DesignPlanar() @@ -223,6 +245,52 @@ def test_renderer_qansysrenderer_default_setup(self): self.assertEqual(default_setup['q3d']['solution_order'], 'High') self.assertEqual(default_setup['q3d']['solver_type'], 'Iterative') + def test_renderer_qelmer_renderer_default_setup(self): + """Test that default_setup in QElmerRenderer have not been accidentally changed.""" + default_setup = QElmerRenderer.default_setup + + self.assertEqual(len(default_setup), 4) + self.assertEqual(len(default_setup['capacitance']), 12) + self.assertEqual(len(default_setup['eigenmode']), 0) + self.assertEqual(len(default_setup['materials']), 2) + self.assertEqual(len(default_setup['constants']), 2) + + self.assertEqual( + default_setup['capacitance']['Calculate_Electric_Field'], True) + self.assertEqual( + default_setup['capacitance']['Calculate_Electric_Energy'], True) + self.assertEqual( + default_setup['capacitance']['Calculate_Capacitance_Matrix'], True) + self.assertEqual( + default_setup['capacitance']['Capacitance_Matrix_Filename'], + 'cap_matrix.txt') + self.assertEqual(default_setup['capacitance']['Linear_System_Solver'], + 'Iterative') + self.assertEqual( + default_setup['capacitance']['Steady_State_Convergence_Tolerance'], + 1e-05) + self.assertEqual( + default_setup['capacitance'] + ['Nonlinear_System_Convergence_Tolerance'], 1e-07) + self.assertEqual( + default_setup['capacitance']['Nonlinear_System_Max_Iterations'], 20) + self.assertEqual( + default_setup['capacitance']['Linear_System_Convergence_Tolerance'], + 1e-10) + self.assertEqual( + default_setup['capacitance']['Linear_System_Max_Iterations'], 500) + self.assertEqual( + default_setup['capacitance']['Linear_System_Iterative_Method'], + 'BiCGStab') + self.assertEqual( + default_setup['capacitance']['BiCGstabl_polynomial_degree'], 2) + + self.assertEqual(default_setup['materials'], ['vacuum', 'silicon']) + + self.assertEqual(default_setup['constants']['Permittivity_of_Vacuum'], + 8.8542e-12) + self.assertEqual(default_setup['constants']['Unit_Charge'], 1.602e-19) + def test_renderer_qq3d_render_options(self): """Test that defaults in QQ3DRenderer were not accidentally changed.""" design = designs.DesignPlanar() @@ -331,6 +399,21 @@ def test_renderer_qgmsh_renderer_options(self): self.assertEqual(options["colors"]["jj"], (84, 140, 168, 150)) self.assertEqual(options["colors"]["dielectric"], (180, 180, 180, 255)) + def test_renderer_qelmer_renderer_options(self): + """Test that default_options in QElmerRenderer were not accidentally + changed.""" + design = designs.MultiPlanar() + renderer = QElmerRenderer(design) + options = renderer.default_options + + self.assertEqual(len(options), 6) + self.assertEqual(options["simulation_type"], "steady_3D") + self.assertEqual(options["simulation_dir"], "./simdata") + self.assertEqual(options["mesh_file"], "out.msh") + self.assertEqual(options["simulation_input_file"], "case.sif") + self.assertEqual(options["postprocessing_file"], "case.msh") + self.assertEqual(options["output_file"], "case.result") + def test_renderer_ansys_renderer_get_clean_name(self): """Test get_clean_name in ansys_renderer.py""" self.assertEqual(ansys_renderer.get_clean_name('name12'), 'name12') @@ -365,6 +448,12 @@ def test_renderer_qgmsh_renderer_name(self): renderer = QGmshRenderer(design, initiate=False) self.assertEqual(renderer.name, 'gmsh') + def test_renderer_qelmer_renderer_name(self): + """Test name in QElmerRenderer.""" + design = designs.MultiPlanar() + renderer = QElmerRenderer(design, initiate=False) + self.assertEqual(renderer.name, 'elmer') + def test_renderer_gdsrenderer_inclusive_bound(self): """Test functionality of inclusive_bound in gds_renderer.py.""" design = designs.DesignPlanar() diff --git a/requirements.txt b/requirements.txt index 534642cef..92121f6b5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -17,4 +17,5 @@ shapely==2.0.1 scqubits==3.1.0 gmsh==4.11.1 pyaedt==0.6.46 +pyyaml==6.0 # jupyter (if you need a fresh install) or ipykernel (if you prefer to make of this a new kernel to use from an existing jupyter install) diff --git a/tutorials/4 Analysis/B. Advanced - Direct use of the renderers/4.19 Analyze a transmon using ElmerFEM.ipynb b/tutorials/4 Analysis/B. Advanced - Direct use of the renderers/4.19 Analyze a transmon using ElmerFEM.ipynb new file mode 100644 index 000000000..ecdf34c46 --- /dev/null +++ b/tutorials/4 Analysis/B. Advanced - Direct use of the renderers/4.19 Analyze a transmon using ElmerFEM.ipynb @@ -0,0 +1,642 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "075bd937-5d0c-4e29-8351-faa413a5f7d1", + "metadata": {}, + "source": [ + "# Analyse a Transmon Qubit using ElmerFEM\n", + "This notebook demonstrates the new open-source rendering and simulation capabilities of Qiskit Metal using `Gmsh` and `ElmerFEM` for tuning the performance parameters of a transmon qubit. instructions to download and install this branch of Qiskit Metal with Gmsh and ElmerFEM, can he found [here](https://github.com/Qiskit/qiskit-metal/blob/elmer_renderer/tutorials/4%20Analysis/B.%20Advanced%20-%20Direct%20use%20of%20the%20renderers/Gmsh-Elmer-install-instructions.md) The tutorial has the following steps:\n", + "\n", + "## Contents\n", + "### 1. Creating a Transmon Qubit in Qiskit Metal\n", + "### 2. Rendering and meshing your design in the `QGmshRenderer`\n", + "- Rendering the design wireframe\n", + "- Applying a basic mesh\n", + "- Customizing the mesh\n", + " - Setting initial mesh size parameters\n", + " - Using the Intelli-mesh feature in `QGmshRenderer`\n", + "- Saving your mesh to a file\n", + "\n", + "### 3. Rendering your design in `QElmerRenderer` (uses `QGmshRenderer`)\n", + "- Render design and generate mesh\n", + "- Export the mesh\n", + "- Run the capacitance simulation\n", + " - Add solution setup\n", + " - Run the solver\n", + " - Export capacitance matrix\n", + "\n", + "### 4. Perform LOM 2.0 Analysis\n", + "- Import the capacitance matrix\n", + "- Define cells and subsystems\n", + "- Define the composite system\n", + "- Compute the results:\n", + " - Qubit Frequency ($f_Q$)\n", + " - Anharmonicity ($\\alpha$)\n", + " - Readout $\\chi$\n", + " - Coupling $g$" + ] + }, + { + "cell_type": "markdown", + "id": "be3cbb43-2c24-4139-b690-1d56047a03c5", + "metadata": {}, + "source": [ + "# Necessary Imports" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28f8c32e-5341-49d6-b63f-0e16419b162c", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "# Import basic things of Qiskit Metal\n", + "from qiskit_metal import MetalGUI, designs\n", + "from qiskit_metal.qlibrary.qubits.transmon_pocket_6 import TransmonPocket6\n", + "\n", + "# Import the Gmsh renderer\n", + "from qiskit_metal.renderers.renderer_gmsh.gmsh_renderer import QGmshRenderer\n", + "\n", + "# Import the Elmer renderer\n", + "from qiskit_metal.renderers.renderer_elmer.elmer_renderer import QElmerRenderer\n", + "from qiskit_metal.renderers.renderer_elmer.elmer_renderer import load_capacitance_matrix_from_file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a76242c0-3b42-4157-bbcb-a99b3c47937a", + "metadata": {}, + "outputs": [], + "source": [ + "%metal_heading 1. Create a Transmon Qubit in Qiskit Metal" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31a7b793-64a5-44c6-9176-718f9b7c1c66", + "metadata": {}, + "outputs": [], + "source": [ + "# Instantiate a design object\n", + "design = designs.MultiPlanar({}, True)\n", + "\n", + "# Invoke the Qiskit Metal GUI\n", + "gui = MetalGUI(design)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fe517d6-d32f-417c-a8fa-33061c2e54c6", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a dictionary for connection pad options for the transmon\n", + "conn_pads = dict(\n", + " connection_pads = dict(\n", + " readout = dict(loc_W=0, loc_H=-1),\n", + " coupler1 = dict(loc_W=-1, loc_H=1),\n", + " coupler2 = dict(loc_W=1, loc_H=1)\n", + " )\n", + ")\n", + "\n", + "# Create a TransmonPocket6 object\n", + "q1 = TransmonPocket6(design, \"Q1\", options=dict(**conn_pads))\n", + "\n", + "# Rebuild and autoscale the GUI\n", + "gui.rebuild()\n", + "gui.autoscale()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d97e061c-fbcc-4419-a30c-8854b8c63444", + "metadata": {}, + "outputs": [], + "source": [ + "%metal_heading 2. Rendering and meshing your design in the `QGmshRenderer`" + ] + }, + { + "cell_type": "markdown", + "id": "2f118489-71a8-494f-8c1b-b6cc837382c3", + "metadata": {}, + "source": [ + "## Rendering the design wireframe" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "713202e1-dedc-469d-8dc3-9fe99496e42e", + "metadata": {}, + "outputs": [], + "source": [ + "# Instantiate the Gmsh renderer\n", + "gmsh_renderer = QGmshRenderer(design)\n", + "\n", + "# Render the design\n", + "# Tip: `mesh_geoms = False` will not mesh the design,\n", + "# but only draw the wire-frame of the geomtries\n", + "gmsh_renderer.render_design(open_pins=[(\"Q1\", \"coupler1\"), (\"Q1\", \"coupler2\"), (\"Q1\", \"readout\")],\n", + " mesh_geoms=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eef5a352-24f3-4257-968f-be197eeba30d", + "metadata": {}, + "outputs": [], + "source": [ + "# Launch Gmsh GUI to see the rendererd design\n", + "gmsh_renderer.launch_gui()" + ] + }, + { + "cell_type": "markdown", + "id": "cea57dd6-816f-4649-a77c-e9b774dcc153", + "metadata": { + "tags": [] + }, + "source": [ + "## Applying a basic mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "545d393c-8d57-4914-9dbc-92e8a0f32b0e", + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [], + "source": [ + "# Add a basic mesh to the design and luanch the GUI to view\n", + "gmsh_renderer.add_mesh(dim=3, intelli_mesh=False)\n", + "gmsh_renderer.launch_gui()" + ] + }, + { + "cell_type": "markdown", + "id": "6347235c-d319-4c61-a981-7e4005b25d15", + "metadata": { + "tags": [] + }, + "source": [ + "## Customizing the mesh\n", + "### Setting initial mesh size parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b851f19f-e8d0-479e-b0ae-056fa8a47e9f", + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [], + "source": [ + "# Update the mesh options\n", + "gmsh_renderer.options.mesh.min_size = '5um'\n", + "gmsh_renderer.options.mesh.max_size = '20um'\n", + "\n", + "# Render the design\n", + "gmsh_renderer.render_design(open_pins=[(\"Q1\", \"coupler1\"), (\"Q1\", \"coupler2\"), (\"Q1\", \"readout\")],\n", + " mesh_geoms=False)\n", + "\n", + "# Mesh the design (without intelli-mesh)\n", + "gmsh_renderer.add_mesh(intelli_mesh=False)\n", + "gmsh_renderer.launch_gui()" + ] + }, + { + "cell_type": "markdown", + "id": "44bc0221-43d2-4551-b69d-4ab4f3f5b416", + "metadata": { + "tags": [] + }, + "source": [ + "### Using the Intelli-mesh feature in QGmshRenderer" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b02fbf9-2717-45b8-9293-3bd1121862f9", + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [], + "source": [ + "# Update the mesh options\n", + "gmsh_renderer.options.mesh.min_size = '5um'\n", + "gmsh_renderer.options.mesh.max_size = '50um'\n", + "\n", + "# Render and mesh the design (with intelli-mesh)\n", + "gmsh_renderer.render_design(open_pins=[(\"Q1\", \"coupler1\"), (\"Q1\", \"coupler2\"), (\"Q1\", \"readout\")],\n", + " mesh_geoms=True)\n", + "gmsh_renderer.launch_gui()" + ] + }, + { + "cell_type": "markdown", + "id": "f441adcc-510f-4c8a-8516-1caf7dcb9166", + "metadata": {}, + "source": [ + "## Saving your mesh to a file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af6556de-1546-47b5-8f26-9227af19ff86", + "metadata": {}, + "outputs": [], + "source": [ + "# Export the Gmsh generated mesh to a file\n", + "gmsh_renderer.export_mesh(\"test.msh\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3672846f-b035-4eab-80a9-3c199eb4fd7c", + "metadata": {}, + "outputs": [], + "source": [ + "# Close the Gmsh renderer\n", + "gmsh_renderer.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88883f88-342e-449c-a1ce-c253a28c0145", + "metadata": {}, + "outputs": [], + "source": [ + "%metal_heading 3. Rendering your design in `QElmerRenderer` (uses `QGmshRenderer`)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c2f7a6a-f6ed-4892-96c3-a62f866ef37c", + "metadata": {}, + "outputs": [], + "source": [ + "# Instantiate the Elmer renderer\n", + "elmer_renderer = QElmerRenderer(design)\n", + "\n", + "# Elmer renderer uses the Gmsh renderer to\n", + "# generate a mesh for the design\n", + "# Set initial parameters for meshing in Gmsh (IMPORTANT step!!)\n", + "elmer_renderer.gmsh.options.mesh.min_size = \"5um\"\n", + "elmer_renderer.gmsh.options.mesh.max_size = \"50um\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7c297c2-be7e-495b-b98a-4062cab91be0", + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [], + "source": [ + "# Render the design\n", + "elmer_renderer.render_design(open_pins=[(\"Q1\", \"coupler1\"), (\"Q1\", \"coupler2\"), (\"Q1\", \"readout\")],\n", + " skip_junctions=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43145408-a90f-4e75-921d-ca8eafd5ea7e", + "metadata": {}, + "outputs": [], + "source": [ + "# View the generated mesh\n", + "elmer_renderer.launch_gmsh_gui()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab7e7c7a-9207-4943-bf0f-1d183a8d007c", + "metadata": {}, + "outputs": [], + "source": [ + "# Export the generated mesh\n", + "elmer_renderer.export_mesh()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5cc8b22f-f281-45c3-a53a-f2a7c9a29fed", + "metadata": {}, + "outputs": [], + "source": [ + "# Add a solution setup to solve for the capacitance matrix\n", + "elmer_renderer.add_solution_setup('capacitance')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "020dc57e-d668-4d7b-b8fc-62fcceb6468d", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the simulation in ElmerFEM\n", + "elmer_renderer.run('capacitance')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b250b918-1cd7-429d-b015-92b114bf2594", + "metadata": {}, + "outputs": [], + "source": [ + "# Display the capacitnce matrix obtained after the simulation\n", + "elmer_renderer.capacitance_matrix" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "07625769-38c4-4be7-9c64-9e4e7b993506", + "metadata": {}, + "outputs": [], + "source": [ + "# Run this, and in the Gmsh window,\n", + "# Right click on an empty area, click on \"Toggle mesh visibility\"\n", + "# Next, go to: Tools -> Visibility, select \"Physical groups\" in drop down menu\n", + "# Select \"Volume 2\", click on apply\n", + "# On the left pane, click on \"Post-processing\", and select the views that you want to observe\n", + "elmer_renderer.display_post_processing_data()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca5cf807-4323-48d8-b6cd-f5e38f3b24d3", + "metadata": {}, + "outputs": [], + "source": [ + "# Export the capacitance matrix\n", + "elmer_renderer.save_capacitance_matrix(\"cap_matrix.txt\")\n", + "\n", + "# Close the elmer renderer\n", + "elmer_renderer.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b20674ae-918f-4aa7-b4f7-1a337ae714a5", + "metadata": {}, + "outputs": [], + "source": [ + "%metal_heading 4. Perform LOM 2.0 Analysis" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "603f8d95-ac34-40e0-9574-53d0f6cde2b8", + "metadata": {}, + "outputs": [], + "source": [ + "# Import necessary things for LOM 2.0 Analysis\n", + "import numpy as np\n", + "from qiskit_metal.analyses.quantization.lom_core_analysis import CompositeSystem, Cell, Subsystem \n", + "from scipy.constants import speed_of_light as c_light\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf3df357-f1fe-4bd0-9510-02fd7e2a55bf", + "metadata": {}, + "outputs": [], + "source": [ + "# Load the saved capacitance matrix\n", + "cap_matrix = load_capacitance_matrix_from_file(\"cap_matrix.txt\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9705e531-ef68-4b6b-bbc6-6a8b69ba0a12", + "metadata": {}, + "outputs": [], + "source": [ + "# Define cells\n", + "opt1 = dict(\n", + " node_rename = {'Q1_coupler1_connector_pad': 'coupler1', 'Q1_coupler2_connector_pad': 'coupler2', 'Q1_readout_connector_pad': 'readout'}, \n", + " cap_mat = cap_matrix,\n", + " ind_dict = {('Q1_pad_top', 'Q1_pad_bot'): 12}, # junction inductance in nH\n", + " jj_dict = {('Q1_pad_top', 'Q1_pad_bot'): 'j1'},\n", + " cj_dict = {('Q1_pad_top', 'Q1_pad_bot'): 2}, # junction capacitance in fF\n", + ")\n", + "cell_1 = Cell(opt1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d5c8c296-13f9-4ec9-8daf-96f2f882b923", + "metadata": {}, + "outputs": [], + "source": [ + "# Define subsystems\n", + "# subsystem 1: Transmon\n", + "transmon = Subsystem(name='transmon1_sys', sys_type='TRANSMON', nodes=['j1'])\n", + "\n", + "# subsystem 2: Coupler 1\n", + "q_opts = dict(\n", + " f_res = 7.4, # resonator dressed frequency in GHz\n", + " Z0 = 50, # characteristic impedance in Ohm\n", + " vp = 0.404314 * c_light # phase velocity \n", + ")\n", + "coup1 = Subsystem(name='coup1_sys', sys_type='TL_RESONATOR', nodes=['coupler1'], q_opts=q_opts)\n", + "\n", + "# subsystem 3: Coupler 2\n", + "q_opts = dict(\n", + " f_res = 7.2, # resonator dressed frequency in GHz\n", + " Z0 = 50, # characteristic impedance in Ohm\n", + " vp = 0.404314 * c_light # phase velocity \n", + ")\n", + "coup2 = Subsystem(name='coup2_sys', sys_type='TL_RESONATOR', nodes=['coupler2'], q_opts=q_opts)\n", + "\n", + "# subsystem 4: Readout\n", + "q_opts = dict(\n", + " f_res = 7, # resonator dressed frequency in GHz\n", + " Z0 = 50, # characteristic impedance in Ohm\n", + " vp = 0.404314 * c_light # phase velocity \n", + ")\n", + "readout = Subsystem(name='readout_sys', sys_type='TL_RESONATOR', nodes=['readout'], q_opts=q_opts)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f44c894-aa2b-4632-b513-da391f101b93", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the composite system\n", + "composite_sys = CompositeSystem(\n", + " subsystems=[transmon, coup1, coup2, readout], \n", + " cells=[cell_1], \n", + " grd_node='ground_plane',\n", + " nodes_force_keep=['coupler1', 'coupler2', 'readout']\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86f39d01-c88f-4300-987f-0384703daf61", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the circuit graph object\n", + "cg = composite_sys.circuitGraph()\n", + "print(cg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19fee4da-5659-4b3b-ad29-182565ee8933", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a hilbert space with the composite system\n", + "hilbertspace = composite_sys.create_hilbertspace()\n", + "print(hilbertspace)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c18d596f-1778-4bf1-852c-6a3c7c8dc927", + "metadata": {}, + "outputs": [], + "source": [ + "# Add interaction between the subsystems\n", + "hilbertspace = composite_sys.add_interaction()\n", + "\n", + "# Get the total hamiltonian of the composite system\n", + "hilbertspace.hamiltonian()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b4c2cdfb-6ede-46eb-bbaf-c19f75b9a5d5", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the reults from the hamiltonian\n", + "# Qubit frequency\n", + "# Chi matrix (having anharmonicity and readout chi values)\n", + "hamiltonian_results = composite_sys.hamiltonian_results(hilbertspace, evals_count=30)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c01e33b-c8fa-464a-bd02-fa8e87f48b5a", + "metadata": {}, + "outputs": [], + "source": [ + "chi_df = hamiltonian_results['chi_in_MHz'].to_dataframe()\n", + "print(\"Transmon frequency :\", hamiltonian_results['fQ_in_Ghz']['transmon1_sys'], \"GHz\")\n", + "print(\"Transmon Anharmonicity :\", chi_df['transmon1_sys']['transmon1_sys'], \"MHz\")\n", + "print(\"Readou chi :\", chi_df['readout_sys']['transmon1_sys'], \"MHz\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "afdfe557-8d48-4a18-85bf-923aa2b63c07", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the coupling 'g' values between qubit and resonators\n", + "composite_sys.compute_gs().to_dataframe()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "963212c4-b31f-48c0-ac56-0058e477b9e0", + "metadata": {}, + "outputs": [], + "source": [ + "# Hamiltonian parameters for the transmon\n", + "transmon.h_params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44627014-50d7-4f5c-b873-c7c438d6c6da", + "metadata": {}, + "outputs": [], + "source": [ + "# Close the Qiskit Metal GUI\n", + "gui.main_window.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c3428fa-7009-451d-9032-a54ecf888791", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}