diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index fd89037d..b30eeffc 100644 --- a/.github/workflows/build.yml +++ b/.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/action@v3.0.0 + - uses: pre-commit/action@v3.0.1 flake8-lint: runs-on: ubuntu-latest diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 3d2d17f8..110c3e45 100644 --- a/.github/workflows/docs.yml +++ b/.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 diff --git a/docs/_toc.yml b/docs/_toc.yml index 00ea7c7b..3c27db04 100644 --- a/docs/_toc.yml +++ b/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: diff --git a/docs/julia/thermal_fill.jl b/docs/julia/thermal_fill.jl new file mode 100644 index 00000000..af329b21 --- /dev/null +++ b/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), +) diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py index 74ac9924..13349303 100644 --- a/docs/photonics/examples/effective_area.py +++ b/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() diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 68bffc52..7d8548cb 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -69,6 +69,41 @@ def n_eff(self): """Effective refractive index of the mode""" return self.k / self.k0 + @cached_property + def poynting(self): + """Poynting vector of the mode""" + + # Extraction of the fields + (Ex, Ey), Ez = self.basis.interpolate(self.E) + (Hx, Hy), Hz = self.basis.interpolate(self.H) + + # New basis will be the discontinuous variant used by the solved Ez-field + poynting_basis = self.basis.with_element( + ElementVector(ElementDG(self.basis.elem.elems[1]), 3) + ) + + # Calculation of the Poynting vector + Px = Ey * Hz - Ez * Hy + Py = Ez * Hx - Ex * Hz + Pz = Ex * Hy - Ey * Hx + + # Projection of the Poynting vector on the new basis + P_proj = poynting_basis.project(np.stack([Px, Py, Pz], axis=0), dtype=np.complex64) + + return poynting_basis, P_proj + + def Sx(self): + basis, _P = self.poynting + return basis.split_bases()[0], _P[basis.split_indices()[0]] + + def Sy(self): + basis, _P = self.poynting + return basis.split_bases()[1], _P[basis.split_indices()[1]] + + def Sz(self): + basis, _P = self.poynting + return basis.split_bases()[2], _P[basis.split_indices()[2]] + @cached_property def te_fraction(self): """TE-fraction of the mode"""