Merge branch 'main' into update-julia

.github/workflows/build.yml

@@ -28,7 +28,7 @@ jobs:
       - name: Install JuliaFormatter and format
         run: |
           julia  -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))'
-      - uses: pre-commit/[email protected].0
+      - uses: pre-commit/[email protected].1
 
   flake8-lint:
     runs-on: ubuntu-latest

.github/workflows/docs.yml

@@ -43,8 +43,6 @@ jobs:
           mamba install julia
           pip install jupyter-book
           pip install --upgrade jupytext
-          pip install sphinx=="5.3.0"
-          pip install docutils=="0.17.1"
           pip install gdsfactory gplugins
           pip list
       
@@ -60,7 +58,7 @@ jobs:
         run: |
           pip install juliacall
           julia -e 'using Pkg; Pkg.add("PythonCall")'
-      - uses: actions/cache@v3
+      - uses: actions/cache@v4
         with:
           path: |
             docs/_build/.jupyter_cache

docs/_toc.yml

@@ -52,6 +52,7 @@ parts:
     - file: julia/lithium_niobate_phase_shifter.jl
     - file: julia/thermal_simple.jl
     - file: julia/heater_3d.jl
+    - file: julia/thermal_fill.jl
 
   - caption: Benchmarks
     chapters:

docs/julia/thermal_fill.jl

@@ -0,0 +1,212 @@
+# ---
+# jupyter:
+#   jupytext:
+#     custom_cell_magics: kql
+#     formats: jl:percent,ipynb
+#     text_representation:
+#       extension: .jl
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.11.2
+#   kernelspec:
+#     display_name: base
+#     language: julia
+#     name: julia-1.9
+# ---
+
+# %% [markdown]
+# # Effective thermal conductivity
+
+# Usually, after the design of the circuitry, the empty space on chips gets filled with equally spaced small rectangles in each layer.
+# Optically, these structures are supposed to be far enough away that their influence on the structures can be neglected.
+# But for thermal considerations, those fill structures can have an impact on the temperature distribution on the chip and thus e.g. on the crosstalk between thermal phase shifters.
+# As it's computationally challenging to include all the small cuboids in the model (which is especially for the meshing a major challenge),
+# a preferable approach is to consider the filled area as a homogenous area of higher thermal conductivity.
+# For this, we calculate the effective thermal conductivity of the filled area by examining a single unit cell.
+
+# To have an intuitively understandable problem, we consider half of the unit cell to be filled with a highly thermally conductive material (metal/silicon) surrounded by a material with low thermal conductance (e.g. silicon dioxide)
+# Let's start with the mesh!
+
+# %% 
+
+unitcell_width = 2
+unitcell_length = 2
+unitcell_thickness = 2
+
+fill_width = 2
+fill_length = 1
+fill_thickness = 2
+
+# %% tags=["hide-input", "thebe-init", "remove-stderr"]
+using PyCall
+np = pyimport("numpy")
+shapely = pyimport("shapely")
+shapely.affinity = pyimport("shapely.affinity")
+OrderedDict = pyimport("collections").OrderedDict
+meshwell = pyimport("meshwell")
+Prism = pyimport("meshwell.prism").Prism
+Model = pyimport("meshwell.model").Model
+
+model = Model()
+
+fill_polygon = shapely.box(
+    xmin = unitcell_length / 2 - fill_length / 2,
+    ymin = unitcell_width / 2 - fill_width / 2,
+    xmax = unitcell_length / 2 + fill_length / 2,
+    ymax = unitcell_width / 2 + fill_width / 2,
+)
+
+fill = Prism(
+    polygons = fill_polygon,
+    buffers = Dict(0 => 0.0, fill_thickness => 0.0),
+    model = model,
+    physical_name = "fill",
+    mesh_order = 1,
+    resolution = Dict("resolution" => 0.2, "SizeMax" => 1.0, "DistMax" => 1.0),
+)
+
+unitcell_polygon =
+    shapely.box(xmin = 0, ymin = 0, xmax = unitcell_length, ymax = unitcell_width)
+unitcell_buffer = Dict(0 => 0.0, unitcell_thickness => 0.0)
+
+unitcell = Prism(
+    polygons = unitcell_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "unitcell",
+    mesh_order = 2,
+)
+
+"""
+BOUNDARIES
+"""
+
+right_polygon = shapely.box(
+    xmin = unitcell_length,
+    ymin = 0,
+    xmax = unitcell_length + 1,
+    ymax = unitcell_width,
+)
+right = Prism(
+    polygons = right_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "right",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+left_polygon = shapely.box(xmin = -1, ymin = 0, xmax = 0, ymax = unitcell_width)
+left = Prism(
+    polygons = left_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "left",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+up_polygon = shapely.box(
+    xmin = 0,
+    ymin = unitcell_width,
+    xmax = unitcell_length,
+    ymax = unitcell_width + 1,
+)
+up = Prism(
+    polygons = up_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "up",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+down_polygon = shapely.box(xmin = 0, ymin = -1, xmax = unitcell_length, ymax = 0)
+down = Prism(
+    polygons = down_polygon,
+    buffers = unitcell_buffer,
+    model = model,
+    physical_name = "down",
+    mesh_bool = false,
+    mesh_order = 0,
+)
+
+"""
+ASSEMBLE AND NAME ENTITIES
+"""
+entities = [unitcell, fill, up, down, left, right]
+
+mesh = model.mesh(
+    entities_list = entities,
+    verbosity = 0,
+    filename = "mesh.msh",
+    default_characteristic_length = 0.2,
+)
+
+
+
+# %% tags=["hide-input", "thebe-init", "remove-stderr"]
+using Gridap
+using GridapGmsh
+using Gridap.Geometry
+using GridapMakie, CairoMakie
+
+using Femwell.Maxwell.Electrostatic
+using Femwell.Thermal
+
+# %% [markdown]
+# For simplicity, we define the thermal conductivity of the unitcell to be 1 and the thermal conductivity of the fill to be 100, which is almost negligible in comparison.
+
+# %%
+thermal_conductivities = ["unitcell" => 1, "fill" => 100.0]
+
+model = GmshDiscreteModel("mesh.msh")
+Ω = Triangulation(model)
+dΩ = Measure(Ω, 1)
+
+labels = get_face_labeling(model)
+tags = get_face_tag(labels, num_cell_dims(model))
+τ = CellField(tags, Ω)
+
+thermal_conductivities =
+    Dict(get_tag_from_name(labels, u) => v for (u, v) in thermal_conductivities)
+ϵ_conductivities(tag) = thermal_conductivities[tag]
+
+# %% [markdown]
+# We define the temperatures at both sides of the unitcell to define the direction in which we want to estimate the thermal conductivity.
+# In other directions the boundaries are considered to be insulating/mirroring.
+
+# %% tags=["remove-stderr", "hide-output"]
+boundary_temperatures = Dict("left___unitcell" => 100.0, "right___unitcell" => 0.0)
+
+
+# %% [markdown]
+# We start with calculating the temperature distribution.
+# From this, we calculate, analog to how we calculate resistances for electrical simulations first the integral
+# $\int σ \left|\frac{\mathrm{d}T}{\mathrm{d}\vec{x}}\right|^2 dA which gives an equivalent of the power
+# and calculate from this the effective thermal conductivity.
+# We expect the thermal conductivity almost twice as high as the thermal conductivity of the material with lower thermal conductivity.
+
+# %% tags=["remove-stderr"]
+T0 = calculate_temperature(
+    ϵ_conductivities ∘ τ,
+    CellField(x -> 0, Ω),
+    boundary_temperatures,
+    order = 2,
+)
+temperature_difference = abs(sum(values(boundary_temperatures) .* [-1, 1]))
+power = abs(
+    sum(
+        ∫(
+            (ϵ_conductivities ∘ τ) *
+            (norm ∘ gradient(temperature(T0))) *
+            (norm ∘ gradient(temperature(T0))),
+        )dΩ,
+    ),
+)
+
+println(
+    "Effective thermal conductivity: ",
+    (power / temperature_difference^2) /
+    ((unitcell_thickness * unitcell_width) / unitcell_length),
+)

docs/photonics/examples/effective_area.py

@@ -22,7 +22,7 @@
 import shapely.affinity
 from matplotlib.ticker import MultipleLocator
 from shapely.ops import clip_by_rect
-from skfem import Basis, ElementTriP0
+from skfem import Basis, ElementTriP0, Functional
 from skfem.io.meshio import from_meshio
 
 from femwell.maxwell.waveguide import compute_modes
@@ -48,10 +48,12 @@
 neff_dict = dict()
 aeff_dict = dict()
 tm_dict = dict()
+p_dict = dict()
 for h in h_list:
     neff_list = []
     aeff_list = []
     tm_list = []
+    p_list = []
     for width in w_list:
         width = width * 1e-3
         nc = shapely.geometry.box(
@@ -95,12 +97,42 @@
                 neff_list.append(np.real(mode.n_eff))
                 aeff_list.append(mode.calculate_effective_area())
                 tm_list.append(mode.transversality)
+                p_list.append(mode.Sz)
+
+                @Functional
+                def I(w):
+                    return 1
+
+                @Functional
+                def Sz(w):
+                    return w["Sz"]
+
+                @Functional
+                def Sz2(w):
+                    return w["Sz"] ** 2
+
+                Sz_basis, Sz_vec = mode.Sz
+
+                print(
+                    "int(Sz)", Sz.assemble(Sz_basis, Sz=Sz_basis.interpolate(Sz_vec))
+                )  # 1 as it's normalized
+                print("int_core(1)", I.assemble(Sz_basis.with_elements("core")))  # area of core
+                print(
+                    "int_core(Sz)",
+                    Sz.assemble(
+                        Sz_basis.with_elements("core"),
+                        Sz=Sz_basis.with_elements("core").interpolate(Sz_vec),
+                    ),
+                )
+                print("int(Sz^2)", Sz2.assemble(Sz_basis, Sz=Sz_basis.interpolate(Sz_vec)))
+
                 break
         else:
             print(f"no TM mode found for {width}")
     neff_dict[str(h)] = neff_list
     aeff_dict[str(h)] = aeff_list
     tm_dict[str(h)] = tm_list
+    p_dict[str(h)] = p_list
 # %% [markdown]
 # Plot the result
 
@@ -184,3 +216,52 @@
 ax4.legend()
 
 plt.show()
+
+# %% [markdown]
+# We can also plot the modes to compare them with the reference article
+
+# %%
+# Data figure h
+w_fig_h = 500
+idx_fig_h = min(range(len(w_list)), key=lambda x: abs(w_list[x] - w_fig_h))
+h_fig_h = 0.7
+basis_fig_h, Pz_fig_h = p_dict[str(h_fig_h)][idx_fig_h]
+
+# Data figure f
+w_fig_f = 600
+idx_fig_f = min(range(len(w_list)), key=lambda x: abs(w_list[x] - w_fig_f))
+h_fig_f = 0.5
+basis_fig_f, Pz_fig_f = p_dict[str(h_fig_f)][idx_fig_f]
+
+fig, ax = plt.subplots(1, 2)
+
+basis_fig_h.plot(np.abs(Pz_fig_h), ax=ax[0], aspect="equal")
+ax[0].set_title(
+    f"Poynting vector $S_z$\nfor h = {h_fig_h}μm & w = {w_fig_h}nm\n(Reproduction of Fig.1.h)"
+)
+ax[0].set_xlim(-w_fig_h * 1e-3 / 2 - 0.1, w_fig_h * 1e-3 / 2 + 0.1)
+ax[0].set_ylim(capital_h - 0.1, capital_h + h_fig_h + 0.1)
+# Turn off the axis wrt to the article figure
+ax[0].axis("off")
+# Add the contour
+for subdomain in basis_fig_h.mesh.subdomains.keys() - {"gmsh:bounding_entities"}:
+    basis_fig_h.mesh.restrict(subdomain).draw(
+        ax=ax[0], boundaries_only=True, color="k", linewidth=1.0
+    )
+
+basis_fig_f.plot(np.abs(Pz_fig_f), ax=ax[1], aspect="equal")
+ax[1].set_title(
+    f"Poynting vector $S_z$\nfor h = {h_fig_f}μm & w = {w_fig_f}nm\n(Reproduction of Fig.1.f)"
+)
+ax[1].set_xlim(-w_list[idx_fig_f] * 1e-3 / 2 - 0.1, w_list[idx_fig_f] * 1e-3 / 2 + 0.1)
+ax[1].set_ylim(capital_h - 0.1, capital_h + h_fig_f + 0.1)
+# Turn off the axis wrt to the article figure
+ax[1].axis("off")
+# Add the contour
+for subdomain in basis_fig_f.mesh.subdomains.keys() - {"gmsh:bounding_entities"}:
+    basis_fig_f.mesh.restrict(subdomain).draw(
+        ax=ax[1], boundaries_only=True, color="k", linewidth=1.0
+    )
+
+fig.tight_layout()
+plt.show()