From 890e65be576099209c374a4c6112436193a00b1a Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Thu, 28 Dec 2023 00:42:36 -0800 Subject: [PATCH 01/12] feat: add calculate_poynting to Mode --- femwell/maxwell/waveguide.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 68bffc52..3d766e6e 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -175,6 +175,26 @@ def E4(w): def calculate_propagation_loss(self, distance): return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance + def calculate_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 = self.basis.with_element(ElementDG(ElementTriP1())) + + # 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 + Px_proj = new_basis.project(Px, dtype=np.complex64) + Py_proj = new_basis.project(Py, dtype=np.complex64) + Pz_proj = new_basis.project(Pz, dtype=np.complex64) + + return new_basis, np.array([Px_proj, Py_proj, Pz_proj]) + def calculate_power(self, elements=None): if not elements: basis = self.basis From d5443de80acab1b901a9b7ae9c59258831bed51d Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Thu, 28 Dec 2023 01:29:03 -0800 Subject: [PATCH 02/12] doc: implement Poynting plot for effective_area example --- docs/photonics/examples/effective_area.py | 45 +++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py index 74ac9924..ca026c95 100644 --- a/docs/photonics/examples/effective_area.py +++ b/docs/photonics/examples/effective_area.py @@ -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,14 @@ neff_list.append(np.real(mode.n_eff)) aeff_list.append(mode.calculate_effective_area()) tm_list.append(mode.transversality) + p_list.append(mode.calculate_poynting()) 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 +188,44 @@ 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, P_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, P_fig_f = p_dict[str(h_fig_f)][idx_fig_f] + +fig, ax = plt.subplots(1, 2) + +basis_fig_h.plot(np.abs(P_fig_h[2]), 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(P_fig_f[2]), 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() From e7615e09303bf44e8c247a7b45fcb7ea8f8348ac Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Thu, 28 Dec 2023 01:54:02 -0800 Subject: [PATCH 03/12] fix: black refactor --- docs/photonics/examples/effective_area.py | 24 +++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py index ca026c95..fc99173b 100644 --- a/docs/photonics/examples/effective_area.py +++ b/docs/photonics/examples/effective_area.py @@ -208,24 +208,32 @@ fig, ax = plt.subplots(1, 2) basis_fig_h.plot(np.abs(P_fig_h[2]), 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_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') +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_h.mesh.restrict(subdomain).draw( + ax=ax[0], boundaries_only=True, color="k", linewidth=1.0 + ) basis_fig_f.plot(np.abs(P_fig_f[2]), 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_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') +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) + basis_fig_f.mesh.restrict(subdomain).draw( + ax=ax[1], boundaries_only=True, color="k", linewidth=1.0 + ) fig.tight_layout() plt.show() From bb320eed2e1949187590fa04be102c790a21960a Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Fri, 29 Dec 2023 00:09:43 -0800 Subject: [PATCH 04/12] feat: poynting calc with use the same element as solved Ez --- femwell/maxwell/waveguide.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 3d766e6e..b02d0662 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -181,6 +181,9 @@ def calculate_poynting(self): # 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 + (Et, Et_basis), (En, En_basis) = self.basis.split(self.E) new_basis = self.basis.with_element(ElementDG(ElementTriP1())) # Calculation of the Poynting vector @@ -350,7 +353,9 @@ def compute_modes( if order == 1: element = ElementTriN1() * ElementTriP1() elif order == 2: - element = ElementTriN2() * ElementTriP2() + # element = ElementTriN2() * ElementTriP2() + from skfem.element import ElementQuad2 + element = ElementTriN2() * ElementQuad2#ElementTriP2() else: raise AssertionError("Only order 1 and 2 implemented by now.") From dc4eb9361d2fb6f7cb61a522efadab2708c685a8 Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Fri, 29 Dec 2023 00:16:06 -0800 Subject: [PATCH 05/12] feat: add poynting as a Mode (cached) property --- docs/photonics/examples/effective_area.py | 2 +- femwell/maxwell/waveguide.py | 46 ++++++++++++----------- 2 files changed, 25 insertions(+), 23 deletions(-) diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py index fc99173b..226cc23d 100644 --- a/docs/photonics/examples/effective_area.py +++ b/docs/photonics/examples/effective_area.py @@ -97,7 +97,7 @@ neff_list.append(np.real(mode.n_eff)) aeff_list.append(mode.calculate_effective_area()) tm_list.append(mode.transversality) - p_list.append(mode.calculate_poynting()) + p_list.append(mode.poynting) break else: print(f"no TM mode found for {width}") diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index b02d0662..b771baee 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -69,6 +69,30 @@ 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 + (Et, Et_basis), (En, En_basis) = self.basis.split(self.E) + new_basis = self.basis.with_element(ElementDG(ElementTriP1())) + + # 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 + Px_proj = new_basis.project(Px, dtype=np.complex64) + Py_proj = new_basis.project(Py, dtype=np.complex64) + Pz_proj = new_basis.project(Pz, dtype=np.complex64) + + return new_basis, np.array([Px_proj, Py_proj, Pz_proj]) + @cached_property def te_fraction(self): """TE-fraction of the mode""" @@ -175,28 +199,6 @@ def E4(w): def calculate_propagation_loss(self, distance): return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance - def calculate_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 - (Et, Et_basis), (En, En_basis) = self.basis.split(self.E) - new_basis = self.basis.with_element(ElementDG(ElementTriP1())) - - # 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 - Px_proj = new_basis.project(Px, dtype=np.complex64) - Py_proj = new_basis.project(Py, dtype=np.complex64) - Pz_proj = new_basis.project(Pz, dtype=np.complex64) - - return new_basis, np.array([Px_proj, Py_proj, Pz_proj]) def calculate_power(self, elements=None): if not elements: From b12affefa8dbbdc53505b904e6a6d7755eb06575 Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Fri, 29 Dec 2023 00:16:36 -0800 Subject: [PATCH 06/12] refactor: black compliant --- femwell/maxwell/waveguide.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index b771baee..3d8d975c 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -199,7 +199,6 @@ def E4(w): def calculate_propagation_loss(self, distance): return -20 / np.log(10) * self.k0 * np.imag(self.n_eff) * distance - def calculate_power(self, elements=None): if not elements: basis = self.basis @@ -357,7 +356,8 @@ def compute_modes( elif order == 2: # element = ElementTriN2() * ElementTriP2() from skfem.element import ElementQuad2 - element = ElementTriN2() * ElementQuad2#ElementTriP2() + + element = ElementTriN2() * ElementQuad2 # ElementTriP2() else: raise AssertionError("Only order 1 and 2 implemented by now.") From 560f62745a2209ceca5cf92a8b3b5e18ecf797ec Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Fri, 29 Dec 2023 00:21:14 -0800 Subject: [PATCH 07/12] fix: change basis name to make it more clear --- femwell/maxwell/waveguide.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 3d8d975c..2625e7f6 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -79,7 +79,7 @@ def poynting(self): # New basis will be the discontinuous variant used by the solved Ez-field (Et, Et_basis), (En, En_basis) = self.basis.split(self.E) - new_basis = self.basis.with_element(ElementDG(ElementTriP1())) + poynting_basis = self.basis.with_element(ElementDG(ElementTriP1())) # Calculation of the Poynting vector Px = Ey * Hz - Ez * Hy @@ -87,11 +87,11 @@ def poynting(self): Pz = Ex * Hy - Ey * Hx # Projection of the Poynting vector on the new basis - Px_proj = new_basis.project(Px, dtype=np.complex64) - Py_proj = new_basis.project(Py, dtype=np.complex64) - Pz_proj = new_basis.project(Pz, dtype=np.complex64) + Px_proj = poynting_basis.project(Px, dtype=np.complex64) + Py_proj = poynting_basis.project(Py, dtype=np.complex64) + Pz_proj = poynting_basis.project(Pz, dtype=np.complex64) - return new_basis, np.array([Px_proj, Py_proj, Pz_proj]) + return poynting_basis, np.array([Px_proj, Py_proj, Pz_proj]) @cached_property def te_fraction(self): From 5ab12bdf67bda30f20498df763d607c9babf16f9 Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Fri, 29 Dec 2023 00:26:52 -0800 Subject: [PATCH 08/12] feat: add poynting component as Mode properties --- docs/photonics/examples/effective_area.py | 10 +++++----- femwell/maxwell/waveguide.py | 15 +++++++++++++++ 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py index 226cc23d..ddafb02a 100644 --- a/docs/photonics/examples/effective_area.py +++ b/docs/photonics/examples/effective_area.py @@ -97,7 +97,7 @@ neff_list.append(np.real(mode.n_eff)) aeff_list.append(mode.calculate_effective_area()) tm_list.append(mode.transversality) - p_list.append(mode.poynting) + p_list.append(mode.Pz) break else: print(f"no TM mode found for {width}") @@ -197,17 +197,17 @@ 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, P_fig_h = p_dict[str(h_fig_h)][idx_fig_h] +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, P_fig_f = p_dict[str(h_fig_f)][idx_fig_f] +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(P_fig_h[2]), ax=ax[0], aspect="equal") +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)" ) @@ -221,7 +221,7 @@ ax=ax[0], boundaries_only=True, color="k", linewidth=1.0 ) -basis_fig_f.plot(np.abs(P_fig_f[2]), ax=ax[1], aspect="equal") +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)" ) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 2625e7f6..7b4e0298 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -93,6 +93,21 @@ def poynting(self): return poynting_basis, np.array([Px_proj, Py_proj, Pz_proj]) + @cached_property + def Px(self): + basis, _P = self.poynting + return basis, _P[0] + + @cached_property + def Py(self): + basis, _P = self.poynting + return basis, _P[1] + + @cached_property + def Pz(self): + basis, _P = self.poynting + return basis, _P[2] + @cached_property def te_fraction(self): """TE-fraction of the mode""" From eaa5a36a0a136f4d4ae5230af53a08a910caf26e Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Sat, 30 Dec 2023 16:25:09 -0800 Subject: [PATCH 09/12] fix: restore ElementTriP2 instead of Quad2 for order=2 --- femwell/maxwell/waveguide.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 7b4e0298..10e383d4 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -369,10 +369,7 @@ def compute_modes( if order == 1: element = ElementTriN1() * ElementTriP1() elif order == 2: - # element = ElementTriN2() * ElementTriP2() - from skfem.element import ElementQuad2 - - element = ElementTriN2() * ElementQuad2 # ElementTriP2() + element = ElementTriN2() * ElementTriP2() else: raise AssertionError("Only order 1 and 2 implemented by now.") From 98cb85359c28d5cf17411b1c0d6874e5c3761cce Mon Sep 17 00:00:00 2001 From: lucasgrjn <18169686+lucasgrjn@users.noreply.github.com> Date: Sat, 30 Dec 2023 16:54:57 -0800 Subject: [PATCH 10/12] feat: for poynting, auto extraction of Ez elements --- femwell/maxwell/waveguide.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 10e383d4..5a720831 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -78,8 +78,7 @@ def poynting(self): (Hx, Hy), Hz = self.basis.interpolate(self.H) # New basis will be the discontinuous variant used by the solved Ez-field - (Et, Et_basis), (En, En_basis) = self.basis.split(self.E) - poynting_basis = self.basis.with_element(ElementDG(ElementTriP1())) + poynting_basis = self.basis.with_element(ElementDG(self.basis.elem.elems[1])) # Calculation of the Poynting vector Px = Ey * Hz - Ez * Hy From 9f0397a6b90e936d5da526daa389336e1c26cb32 Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Wed, 10 Jan 2024 11:01:57 -0800 Subject: [PATCH 11/12] use vectorbasis instead of three individual bases for the poynting vector --- femwell/maxwell/waveguide.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index 5a720831..d6e750e8 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -78,7 +78,9 @@ def poynting(self): (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(ElementDG(self.basis.elem.elems[1])) + poynting_basis = self.basis.with_element( + ElementVector(ElementDG(self.basis.elem.elems[1]), 3) + ) # Calculation of the Poynting vector Px = Ey * Hz - Ez * Hy @@ -86,26 +88,24 @@ def poynting(self): Pz = Ex * Hy - Ey * Hx # Projection of the Poynting vector on the new basis - Px_proj = poynting_basis.project(Px, dtype=np.complex64) - Py_proj = poynting_basis.project(Py, dtype=np.complex64) - Pz_proj = poynting_basis.project(Pz, dtype=np.complex64) + P_proj = poynting_basis.project(np.stack([Px, Py, Pz], axis=0), dtype=np.complex64) - return poynting_basis, np.array([Px_proj, Py_proj, Pz_proj]) + return poynting_basis, P_proj @cached_property def Px(self): basis, _P = self.poynting - return basis, _P[0] + return basis.split_bases()[0], _P[basis.split_indices()[0]] @cached_property def Py(self): basis, _P = self.poynting - return basis, _P[1] + return basis.split_bases()[1], _P[basis.split_indices()[1]] @cached_property def Pz(self): basis, _P = self.poynting - return basis, _P[2] + return basis.split_bases()[2], _P[basis.split_indices()[2]] @cached_property def te_fraction(self): From 464df19cbbd3fc3387a5c912b69fcd4b8463d2fa Mon Sep 17 00:00:00 2001 From: Helge Gehring <42973196+HelgeGehring@users.noreply.github.com> Date: Wed, 21 Feb 2024 01:32:26 -0800 Subject: [PATCH 12/12] add integrals for (1) and (3) --- docs/photonics/examples/effective_area.py | 32 +++++++++++++++++++++-- femwell/maxwell/waveguide.py | 9 +++---- 2 files changed, 33 insertions(+), 8 deletions(-) diff --git a/docs/photonics/examples/effective_area.py b/docs/photonics/examples/effective_area.py index ddafb02a..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 @@ -97,7 +97,35 @@ neff_list.append(np.real(mode.n_eff)) aeff_list.append(mode.calculate_effective_area()) tm_list.append(mode.transversality) - p_list.append(mode.Pz) + 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}") diff --git a/femwell/maxwell/waveguide.py b/femwell/maxwell/waveguide.py index d6e750e8..7d8548cb 100644 --- a/femwell/maxwell/waveguide.py +++ b/femwell/maxwell/waveguide.py @@ -92,18 +92,15 @@ def poynting(self): return poynting_basis, P_proj - @cached_property - def Px(self): + def Sx(self): basis, _P = self.poynting return basis.split_bases()[0], _P[basis.split_indices()[0]] - @cached_property - def Py(self): + def Sy(self): basis, _P = self.poynting return basis.split_bases()[1], _P[basis.split_indices()[1]] - @cached_property - def Pz(self): + def Sz(self): basis, _P = self.poynting return basis.split_bases()[2], _P[basis.split_indices()[2]]