Improve single transmon example (#36)

* docs: Improve single transmon example

* Address review comments and make small edits to single transmon example text

Author: gpeairs

Diff for: docs/src/assets/single_transmon_mesh.png

@@ -1,7 +1,8 @@
 [deps]
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DeviceLayout = "ebf59a4a-04ec-49d7-8cd4-c9382ceb8e85"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JSONSchema = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
-DeviceLayout = "ebf59a4a-04ec-49d7-8cd4-c9382ceb8e85"
+PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"

Diff for: docs/src/examples/singletransmon.md

@@ -8,6 +8,9 @@ import .SchematicDrivenLayout.ExamplePDK
 import .SchematicDrivenLayout.ExamplePDK: LayerVocabulary, L1_TARGET, add_bridges!
 using .ExamplePDK.Transmons, .ExamplePDK.ReadoutResonators
 import .ExamplePDK.SimpleJunctions: ExampleSimpleJunction
+import DeviceLayout: uconvert
+
+using PRIMA
 
 """
     single_transmon(
@@ -34,34 +37,47 @@ function single_transmon(;
     cap_width=24μm,
     cap_length=620μm,
     cap_gap=30μm,
+    total_length=5000μm,
     n_meander_turns=5,
     hanger_length=500μm,
     bend_radius=50μm,
     save_mesh::Bool=false,
-    save_gds::Bool=false
+    save_gds::Bool=false,
+    mesh_order=2
 )
     #### Reset name counter for consistency within a Julia session
     reset_uniquename!()
-    # Compute implicit parameters
-    w_grasp = cap_width + 2 * cap_gap
-    total_height = 444μm + hanger_length + (3 + n_meander_turns * 2) * bend_radius
 
-    #### Create abstract components
+    #### Assemble schematic graph
+    ### Compute additional/implicit parameters
     cpw_width = 10μm
     cpw_gap = 6μm
     PATH_STYLE = Paths.SimpleCPW(cpw_width, cpw_gap)
     BRIDGE_STYLE = ExamplePDK.bridge_geometry(PATH_STYLE)
-
+    coupling_gap = 5μm
+    w_grasp = cap_width + 2 * cap_gap
+    arm_length = 428μm # straight length from meander exit to claw
+    total_height =
+        arm_length +
+        coupling_gap +
+        Paths.extent(PATH_STYLE) +
+        hanger_length +
+        (3 + n_meander_turns * 2) * bend_radius
+    ### Create abstract components
+    ## Transmon
     qubit = ExampleRectangleTransmon(;
         jj_template=ExampleSimpleJunction(),
         name="qubit",
         cap_length,
         cap_gap,
         cap_width
     )
+    ## Resonator
     rres = ExampleClawedMeanderReadout(;
         name="rres",
         coupling_length=400μm,
+        coupling_gap,
+        total_length,
         w_shield,
         w_claw,
         l_claw,
@@ -73,14 +89,6 @@ function single_transmon(;
         bend_radius,
         bridge=BRIDGE_STYLE
     )
-
-    ##### Build schematic graph
-    g = SchematicGraph("single-transmon")
-    qubit_node = add_node!(g, qubit)
-    rres_node = fuse!(g, qubit_node, rres)
-    # Equivalent to `fuse!(g, qubit_node=>:readout, rres=>:qubit)`
-    # because `matching_hooks` was implemented for that component pair
-
     ## Readout path
     readout_length = 2700μm
     p_readout = Path(
@@ -92,22 +100,31 @@ function single_transmon(;
     straight!(p_readout, readout_length / 2, PATH_STYLE)
     straight!(p_readout, readout_length / 2, PATH_STYLE)
 
-    ## Readout lumped ports - one at each end
+    # Readout lumped ports - squares on CPW trace, one at each end
     csport = CoordinateSystem(uniquename("port"), nm)
     render!(
         csport,
         only_simulated(centered(Rectangle(cpw_width, cpw_width))),
         LayerVocabulary.PORT
     )
+    # Attach with port center `cpw_width` from the end (instead of `cpw_width/2`) to avoid corner effects
     attach!(p_readout, sref(csport), cpw_width, i=1) # @ start
     attach!(p_readout, sref(csport), readout_length / 2 - cpw_width, i=2) # @ end
+
+    #### Build schematic graph
+    g = SchematicGraph("single-transmon")
+    qubit_node = add_node!(g, qubit)
+    rres_node = fuse!(g, qubit_node, rres)
+    # Equivalent to `fuse!(g, qubit_node=>:readout, rres=>:qubit)`
+    # because `matching_hooks` was implemented for that component pair
     p_readout_node = add_node!(g, p_readout)
+
     ## Attach resonator to feedline
     # Instead of `fuse!` we use a schematic-based `attach!` method to place it along the path
     # Syntax is a mix of `fuse!` and how we attached the ports above
     attach!(g, p_readout_node, rres_node => :feedline, 0mm, location=1)
 
-    ## Position the components.
+    #### Create the schematic (position the components)
     floorplan = plan(g)
     add_bridges!(floorplan, BRIDGE_STYLE, spacing=300μm) # Add bridges to paths
 
@@ -120,8 +137,11 @@ function single_transmon(;
     chip = centered(Rectangle(substrate_x, substrate_y), on_pt=center_xyz)
     sim_area = centered(Rectangle(substrate_x, substrate_y), on_pt=center_xyz)
 
+    # Define bounds for bounding simulation box
     render!(floorplan.coordinate_system, sim_area, LayerVocabulary.SIMULATED_AREA)
+    # postrendering operations in solidmodel target define metal = (WRITEABLE_AREA - METAL_NEGATIVE) + METAL_POSITIVE
     render!(floorplan.coordinate_system, sim_area, LayerVocabulary.WRITEABLE_AREA)
+    # Define rectangle that gets extruded to generate substrate volume
     render!(floorplan.coordinate_system, chip, LayerVocabulary.CHIP_AREA)
 
     check!(floorplan)
@@ -134,37 +154,43 @@ function single_transmon(;
     # resolution near edges of the geometry.
     meshing_parameters = SolidModels.MeshingParameters(
         mesh_scale=1.0,
-        mesh_order=2,
+        α_default=0.9,
+        mesh_order=mesh_order,
         options=Dict("General.Verbosity" => 1) # General Gmsh option input
     )
     render!(sm, floorplan, tech, meshing_parameters=meshing_parameters)
 
     if save_mesh
+        # SolidModels.gmsh.option.set_number("General.NumThreads", 1) # Force single-threaded (deterministic) meshing
         SolidModels.gmsh.model.mesh.generate(3) # runs without error
         save(joinpath(@__DIR__, "single_transmon.msh2"), sm)
     end
 
     if save_gds
-        # Render to GDS as well.
+        # Render to GDS as well, may be useful to debug SolidModel generation
         c = Cell("single_transmon", nm)
-        render!(c, floorplan, L1_TARGET, strict=:no)
+        # Use simulation=true to render simulation-only geometry, `strict=:no` to continue from errors
+        render!(c, floorplan, L1_TARGET, strict=:no, simulation=true)
         flatten!(c)
         save(joinpath(@__DIR__, "single_transmon.gds"), c)
     end
     return sm
 end
 
 """
-    configfile(sm::SolidModel; palace_build=nothing)
+    configfile(sm::SolidModel; palace_build=nothing, solver_order=2, amr=0)
 
 Given a `SolidModel`, assemble a dictionary defining a configuration file for use within
 Palace.
 
   - `sm`: The `SolidModel`from which to construct the configuration file
   - `palace_build = nothing`: Path to a Palace build directory, used to perform validation of
     the configuration file. If not present, no validation is performed.
+  - `solver_order = 2`: Finite element order (degree) for the solver. Palace supports arbitrary
+    high-order spaces.
+  - `amr = 0`: Maximum number of adaptive mesh refinement (AMR) iterations.
 """
-function configfile(sm::SolidModel; palace_build=nothing)
+function configfile(sm::SolidModel; palace_build=nothing, solver_order=2, amr=0)
     attributes = SolidModels.attributes(sm)
 
     config = Dict(
@@ -175,9 +201,9 @@ function configfile(sm::SolidModel; palace_build=nothing)
         ),
         "Model" => Dict(
             "Mesh" => joinpath(@__DIR__, "single_transmon.msh2"),
-            "L0" => DeviceLayout.ustrip(m, 1SolidModels.STP_UNIT), # um is Palace default; record it anyway
+            "L0" => 1e-6, # um is Palace default; record it anyway
             "Refinement" => Dict(
-                "MaxIts" => 0 # Increase to enable AMR
+                "MaxIts" => amr # Nonzero to enable AMR
             )
         ),
         "Domains" => Dict(
@@ -231,8 +257,8 @@ function configfile(sm::SolidModel; palace_build=nothing)
             ]
         ),
         "Solver" => Dict(
-            "Order" => 1,
-            "Eigenmode" => Dict("N" => 2, "Tol" => 1.0e-6, "Target" => 2, "Save" => 2),
+            "Order" => solver_order,
+            "Eigenmode" => Dict("N" => 2, "Tol" => 1.0e-6, "Target" => 1, "Save" => 2),
             "Linear" => Dict("Type" => "Default", "Tol" => 1.0e-7, "MaxIts" => 500)
         )
     )
@@ -257,12 +283,12 @@ Given a configuration dictionary, write and optionally run the Palace simulation
 
 Writes `config.json` to `@__DIR__` in order to pass the configuration into Palace.
 
-  - `config` - A configuration file defining the required fields for a Palace configuration file
-  - `palace_build` - Path to a Palace build.
-  - `np = 0` - Number of MPI processes to use in the call to Palace. If greater than 0 attempts
+  - `config`: A configuration file defining the required fields for a Palace configuration file
+  - `palace_build`: Path to a Palace build.
+  - `np = 0`: Number of MPI processes to use in the call to Palace. If greater than 0 attempts
     to call palace from within the Julia shell. Requires correct specification of `ENV[PATH]`.
-  - `nt = 1` - Number of OpenMp threads to use in the call to Palace (requires Palace built with
-    OpenMp)
+  - `nt = 1`: Number of OpenMP threads to use in the call to Palace (requires Palace built with
+    OpenMP)
 """
 function palace_job(config::Dict; palace_build, np=0, nt=1)
     # Write the configuration file to json, ready for Palace ingestion
@@ -293,24 +319,136 @@ function palace_job(config::Dict; palace_build, np=0, nt=1)
         freq = CSV.File(joinpath(postprodir, "eig.csv"); header=1) |> DataFrame
 
         println("Eigenmode Frequencies (GHz): ", freq[:, 2])
+        return freq[:, 2]
     end
     return nothing
 end
 
 """
-    main(palace_build, np=0)
+    compute_eigenfrequencies(palace_build, np=1; solver_order=2, mesh_order=2, cap_length=620μm, total_length=5000μm))
 
 Given a build of Palace found at `palace_build`, assemble a `SolidModel` of the single
 transmon example, mesh the geometry, define a configuration file for the geometry, and if
-`np > 0`, launch Palace from the provided `palace_build`.`
+`np > 0`, launch Palace from the provided `palace_build`. (If `0`, use `palace_build`
+only to validate config.)
+
+# Keyword arguments
+
+  - `solver_order = 2`: Finite element order (degree) for the solver. Palace supports arbitrary
+    high-order spaces.
+  - `mesh_order = 2`: Polynomial order used to represent the element geometries in the mesh.
+  - `cap_length = 620μm`: Length of transmon island capacitor.
+  - `total_length = 5000μm`: Total length of readout resonator.
 """
-function main(palace_build, np=0)
+function compute_eigenfrequencies(
+    palace_build,
+    np=1;
+    solver_order=2,
+    mesh_order=2,
+    cap_length=620μm,
+    total_length=5000μm
+)
     # Construct the SolidModel
-    @time "SolidModel + Meshing" sm = single_transmon(save_mesh=true)
+    @time "SolidModel + Meshing" sm = single_transmon(
+        save_mesh=true;
+        cap_length=cap_length,
+        total_length=total_length,
+        mesh_order=mesh_order
+    )
     # Assemble the configuration
-    @time "Configuration" config = configfile(sm; palace_build)
+    @time "Configuration" config = configfile(sm; palace_build, solver_order=solver_order)
     # Call Palace
-    @time "Palace" palace_job(config; palace_build, np)
+    @time "Palace" freqs = palace_job(config; palace_build, np)
+    return freqs
+end
+
+"""
+    frequency_targeting_errfunc(targets_GHz, palace_build, np, solver_order, mesh_order, freq_log, param_log)
+
+Create an error function that can be minimized by an optimization routine to reach target frequencies.
+
+The returned function takes parameters that control the transmon's capacitor length and resonator's total length.
+It runs `compute_eigenfrequencies` and returns the mean squared relative error between the computed
+eigenfrequencies and `targets_GHz`.
+It also pushes frequencies and parameter values to the provided arrays `freq_log` and `param_log`.
+"""
+function frequency_targeting_errfunc(
+    targets_GHz,
+    palace_build,
+    np,
+    solver_order,
+    mesh_order,
+    freq_log,
+    param_log
+)
+    return function errfunc(x, p=()) # Many optimizer interfaces require a second argument for fixed parameters
+        freqs = compute_eigenfrequencies(
+            palace_build,
+            np;
+            cap_length=(1 / x[1]^2) * 620μm, # Transform so that frequency(x) is approximately linear
+            total_length=(1 / x[2]) * 5000μm, # ... and x_i are all ~1
+            solver_order,
+            mesh_order
+        )[1:2] # [1:2] because technically Palace can find more than two eigenfrequencies
+        push!(freq_log, freqs) # Log for later convenience
+        push!(param_log, copy(x)) # Copy because `x` gets reused
+        return sum((freqs .- targets_GHz) .^ 2 ./ targets_GHz .^ 2) / 2 # Mean squared relative error
+    end
+end
+
+"""
+    run_optimization(palace_build, np=1; targets_GHz=[3.0, 4.0], reltol=1e-2, solver_order=2, mesh_order=2)
+
+Optimize the transmon and resonator parameters to achieve target frequencies.
+
+The objective function generates a new schematic, SolidModel, and mesh;
+generates and validates a new Palace configuration file; runs Palace using the mesh
+and configuration file; and returns the mean squared relative error between the
+computed eigenfrequencies and `targets_GHz`.
+The optimization routine stops when the mean squared relative error is less than `reltol^2`.
+
+`np`, `solver_order`, and `mesh_order` are used for configuring and running Palace as in
+`compute_eigenfrequencies`.
+"""
+function run_optimization(
+    palace_build,
+    np=1;
+    targets_GHz=[3.0, 4.0],
+    reltol=1e-2,
+    solver_order=2,
+    mesh_order=2
+)
+    freq_log = []
+    param_log = []
+    errfunc = frequency_targeting_errfunc( # Create the objective function for optimization
+        targets_GHz,
+        palace_build,
+        np,
+        solver_order,
+        mesh_order,
+        freq_log,
+        param_log
+    )
+    final_params, info = prima( # Run the optimization
+        errfunc,
+        [1.0, 1.0]; # Initial parameters
+        ftarget=reltol^2, # Stop when `errfunc(x) < reltol^2`
+        xl=[0.6, 0.6], # Lower bounds
+        xu=[1.4, 1.4], # Upper bounds
+        rhobeg=0.2 # Initial trust region radius
+    )
+    println("""
+        Number of Palace runs: $(info.nf)
+        Initial parameters:
+            Transmon capacitor_length = 620.0μm
+            Resonator total_length = 5000.0μm
+        Initial frequencies: $(round.(first(freq_log), digits=3)) GHz
+        Final parameters:
+            Transmon capacitor_length = $(round(μm, 620.0μm/final_params[1]^2, digits=3))
+            Resonator total_length = $(round(μm, 5000.0μm/final_params[2], digits=3))
+        Final frequencies: $(round.(last(freq_log), digits=3)) GHz
+        """)
+    return final_params, info, freq_log, param_log
 end
 
 end # module