The Issue:
I'm trying to build 0.1.5 inside of gentoo and the pyproject.toml is missing a 'build-backend', making it difficult to determine what build system to use for this repo.

Here is the reported error from the gentoo build system.
Compiling source in /var/tmp/portage/dev-python/femwell-0.1.5/work/femwell-0.1.5 ...
python3_11: running distutils-r1_run_phase distutils-r1_python_compile
ERROR: dev-python/femwell-0.1.5::localrepo failed (compile phase):
Unable to obtain build-backend from pyproject.toml

some interesting links that may allow expanding femwell use cases and integrations

http://onelab.info/photonics/
https://gitlab.onelab.info/gmsh

When I pip install femwell, I get the following: No module named 'femwell.maxwell'

I think the examples use more recent code than the pypi package.

Many projects provide benchmarks to analytical calculations, literature, or other codes

It would be good to include a section in the docs dedicated to showing agreement between femwell and other sources

Internally, we could also build tests around these benchmarks

In for example https://helgegehring.github.io/femwell/julia/waveguide_anisotropic.html, the reference shows both a URL and a doi link, which point to the same thing

I think it's fine for the bibtex file to contain the URLs, but would be nice if the bibtex processor could have logic to only show the doi if it exists and otherwise show the URL.

@HelgeGehring @simbilod @mdecea @joamatab
Hi, I have tried to use the live launch on example code but they fail. Also, I tried to run the codes on Google colab but its seems it fails to import from "OrderedDict".
The following is the error I get. Could you help? Thanks

**_"---------------------------------------------------------------------------

OSError Traceback (most recent call last)

in <cell line: 9>()
7 from skfem.io import from_meshio
8
----> 9 from femwell.mesh import mesh_from_OrderedDict
10 from femwell.thermal import solve_thermal
11

3 frames

/usr/lib/python3.10/ctypes/init.py in init(self, name, mode, handle, use_errno, use_last_error, winmode)
372
373 if handle is None:
--> 374 self._handle = _dlopen(self._name, mode)
375 else:
376 self._handle = handle

OSError: libGLU.so.1: cannot open shared object file: No such file or directory"_**

Adding the possibility to ]have symmetry planes to save on mesh cells to increase computational speed would be a great add-on when fine meshes are required. Perhaps I missed it on the documentation, but I could only spot examples with full 2D mesh even when symmetry is present.

I encountered an issue when solving for 1st quasi TE&TM modes for different widths. The lines represents different widths. There are jumps in neff (flipping TE fraction). But if we follow the TE fraction guidance, the behavior is monotonic.

maybe somehow similar issue in the example https://helgegehring.github.io/femwell/photonics/examples/vary_width.html

Hello,
I'm new to femwell and I was wondering if there's a particular example that I could use to compute the supported modes in an air filled capillary waveguide with a fused silica cladding. The internal diameter of the capillary is 700 microns and outer diameter is 850 micron. I appreciate any help!

-Timothy

I have noticed an issue when using mesh_from_Dict(), stemming from the following lines:

       # Add overall bounding box
        total = unary_union(list(shapes_dict.values()))
        all_polygons = [total, *list(shapes_dict.values())]
        # Break up all shapes so that plane is tiled with non-overlapping layers, get the maximal number of fragments
        # Equivalent to BooleanFragments
        np.seterr(invalid="ignore")
        listpoly = [a.intersection(b) for a, b in combinations(all_polygons, 2)]
        np.seterr(**initial_settings)
        rings = [
            LineString(list(object.exterior.coords))
            for object in listpoly
            if not (object.geom_type in ["Point", "LineString"] or object.is_empty)
        ]

This is due to the fact that listpoly will have a MultiLineString when intersecting a LinearString with total. The change that could (probably) solve this is:

       # Add overall bounding box
        total = unary_union(list(shapes_dict.values()))
        all_polygons = [total, *list(shapes_dict.values())]
        # Break up all shapes so that plane is tiled with non-overlapping layers, get the maximal number of fragments
        # Equivalent to BooleanFragments
        np.seterr(invalid="ignore")
        listpoly = [a.intersection(b) for a, b in combinations(all_polygons, 2)]
        np.seterr(**initial_settings)
        rings = [
            LineString(list(object.exterior.coords))
            for object in listpoly
            if not (object.geom_type in ["Point", "LineString", "MultiLineString"] or object.is_empty)
        ]

Minimal working example

from collections import OrderedDict

import numpy as np

from shapely.geometry import box, LineString, MultiLineString
from skfem.io.meshio import from_meshio

from femwell.mesh import mesh_from_OrderedDict,mesh_from_Dict

box1 = box(-1,-1,1,1)
box2 = box(-1,5,1,6)

surface = box1.buffer(0.5).exterior
polygons = dict(box1 = box1,
               box2 = box2,
               surface = surface)

resolutions = dict(box1 = {'resolution': 0.2, 'distance':1},
                   box2 =  {'resolution': 0.2, 'distance':1},
                   surface =  {'resolution': 0.2, 'distance':1})

mesh = from_meshio(mesh_from_Dict(polygons, resolutions))

mesh.draw()

Hi @HelgeGehring

I am trying to install femwell on a Windows VM. I installed Python 3.11.5, VS Code, and the Python plugin.

In the terminal prompt for Python in VS Code, I tried pip install femwell, and get the following error:

I presume femwell should work on Windows? I was able to successfully pip install gmsh and pygmsh, so at least I have some of the requirements installed.

Thank you
Lukas

The new scipy backend fails with

scipy.sparse.linalg._eigen.arpack.arpack.ArpackError: ARPACK error -9999: Could not build an Arnoldi factorization. IPARAM(5) returns the size of the current Arnoldi factorization. The user is advised to check that enough workspace and array storage has been allocated.

the same problem used to work well with slepc

We could add https://doi.org/10.1038/ncomms8565 to the docs, should be easy to reproduce :)

2015 J Notaros M Popovic - Finite dif complex wavevector band structure solver for periodic radiative phot.pdf

Sharing paper on 2d bloch solver implementation

error: ModuleNotFoundError: No module named 'petsc4py'
I tried -
pip install petsc4py (version - petsc4py 3.18.5), but this package has issues to install (for any version)

Thank you for sharing this nice package for mode analysis!

When running the demo https://helgegehring.github.io/femwell/julia/waveguide.html, when running
modes = calculate_modes(model, ε ∘ τ, λ = 1.55, num = 2, order = 1), I encounter the following error.

ERROR: ┌ Error: XYAUPD_Exception: Maximum number of iterations taken. All possible eigenvalues of OP has been found.
│ IPARAM(5) returns the number of wanted converged Ritz values.
│ info = 1
└ @ Arpack ~/.julia/packages/Arpack/FCvNd/src/libarpack.jl:47

To reproduce it is to clone the repo, instantiate it in julia, install the femwell in python, run the Femwell.jl in REPL to compile the module, change line 63 in docs/julia/waveguide.jl from using Femwell.Maxwell.Waveguide to using Main.Femwell.Maxwell.Waveguide, and run waveguide.jl.

Hi,

I was trying to recreate figure 3b(W = 6 um, H = 0.8 um) from this paper. Strangely, GVD is doest not match when wavelength is greater than 3100nm. This happen when i am trying to recreate other curves in figure 3b as well. I wonder if i missed any parameters for the solver. Any help would be appreciate!

Here is my code:

from femwell.maxwell.waveguide import compute_modes
from skfem import Basis, ElementTriP0
from tqdm import tqdm
from femwell.visualization import plot_domains
from skfem.io import from_meshio
import shapely
import matplotlib.pyplot as plt
from refractive_index import n_MgF2, n_Si3N4, n_Air
from collections import OrderedDict
from numpy.polynomial import Polynomial
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline
import pandas as pd
from scipy.constants import speed_of_light
from femwell.mesh import mesh_from_OrderedDict

# waveguide parameters
width = 6  # um
height = 0.8  # um

n2 = 2.5e-19  # m^2/W n2 is the nonlinear refractive index at the center
Alpha = 0.7  # loss (dB/cm)

wavelength_range = [310, 5100]
wavelegnth_step = 70  #  steps

n_core = n_Si3N4
n_lower_cladding = lambda w: n_MgF2(w, ray="o") 
n_air = n_Air

# Construct waveguide geometry
core = shapely.geometry.box(-width / 2, 0, +width / 2, height)
lower_cladding = shapely.geometry.box(-6, -6, 6, 0)
air = shapely.geometry.box(-6, 0, 6, 6)
polygons = OrderedDict(
    core=core,
    lower_cladding=lower_cladding,
    air= air
)

# Define material property and resolution of waveguide
resolutions = dict(core={"resolution": 0.04, "distance": 0.1},
                   lower_cladding={"resolution": 0.15, "distance": 0.2},
                   air={"resolution": 0.2, "distance": 0.2})

n_dict = {"core": n_core, "lower_cladding": n_lower_cladding, "air": n_air}

# Calculate dispersion and gamma
mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions))
mesh.draw().show()
plot_domains(mesh)
plt.show()

print("start")
# Calculate dispersion and gamma
basis0 = Basis(mesh, ElementTriP0())
epsilon = basis0.zeros()
wavelength_list = np.linspace(wavelength_range[0], wavelength_range[1], wavelegnth_step)
neff_list = []
aeff_list = []
for wavelength in tqdm(wavelength_list):
    wavelength = wavelength * 1e-3
    for subdomain, n in n_dict.items():
        epsilon[basis0.get_dofs(elements=subdomain)] = n(wavelength) ** 2
    modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=3, order=1)
    modes_sorted = modes.sorted(key=lambda mode: -np.real(mode.n_eff))
    mode = modes_sorted[0]
    #mode.show(mode.E.real, direction ="x")
    neff_list.append(np.real(mode.n_eff))
    aeff_list.append(mode.calculate_effective_area())

neff_list = np.array(neff_list)
aeff_list = np.array(aeff_list)
wls = np.array(wavelength_list)

#--------calculate GVD-----------------
fig, axs = plt.subplots(3, 1, figsize=(9, 20))

y_spl = UnivariateSpline(wls, neff_list, s=0, k=3)
x_range = np.linspace(wls[0], wls[-1], 1000)
print(y_spl(1000))

axs[0].set_xlabel("Wavelength / µm")
axs[0].set_ylabel("neff")
axs[0].set_title(" neff vs wavelength fit")
axs[0].semilogy(x_range, y_spl(x_range))
axs[0].semilogy(wls, neff_list, 'ro', label='data')
axs[0].legend()

y_spl_2d = y_spl.derivative(n=2)
axs[1].set_xlabel("Wavelength / µm")
axs[1].set_ylabel("neff''")
axs[1].set_title("second derivative fit")
axs[1].plot(x_range, y_spl_2d(x_range))

ref_gvd = pd.read_csv("../reference_data/GVD.csv", dtype=np.float64)
ref_gvd_x, ref_gvd_y = np.split(ref_gvd.values, 2, axis=1)
axs[2].plot(ref_gvd_x, ref_gvd_y, c="green", label="paper")

GVD = (-wls / (2.99792e-7) * y_spl_2d(wls))
axs[2].plot(wls, GVD, label="calculated", c="red")

axs[2].set_ylabel("GVD")
axs[2].set_ylim(-1500, 250)
axs[2].set_xlim(500, 6000)
axs[2].set_title("GVD paramter")
axs[2].legend()

plt.tight_layout()
plt.show()

The examples are great, but it would be helpful to have better docstrings for functions such as compute_modes

@EliDoris

Code:

%load_ext autoreload
%autoreload 2

import matplotlib.pyplot as plt
import gdsfactory as gf
from tqdm.auto import tqdm
import numpy as np
from gdsfactory.technology import LayerStack



import sys
import logging
from rich.logging import RichHandler
import gdsfactory as gf

from amf_pdk_pkg import PDK, LAYER_STACK, strip_SiN, strip_SiN_only

from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import shapely
import shapely.affinity
from scipy.constants import epsilon_0, speed_of_light
from shapely.ops import clip_by_rect
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio

from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict
from femwell.visualization import plot_domains


gf.config.rich_output()
PDK.activate()

logger = logging.getLogger()
logger.removeHandler(sys.stderr)
logging.basicConfig(level="WARNING", datefmt="[%X]", handlers=[RichHandler()])

filtered_layerstack = LayerStack(
    layers={
        k: LAYER_STACK.layers[k]
        for k in (
            "sin",
            "clad",
            "box",
        )
    }
)

resolutions = {}
resolutions["sin"] = {"resolution": 0.02, "distance": 2}
resolutions['clad'] = {"resolution": 0.2, "distance": 1}
resolutions['box'] = {"resolution": 0.2, "distance": 1}
c = gf.components.straight(cross_section=strip_SiN)
mesh = c.to_gmsh(
    type="uz",
    xsection_bounds=[(5, -5), (5, 5)],
    layer_stack=filtered_layerstack,
    filename=f"mesh.msh",
    resolutions=resolutions
)
mesh = from_meshio(mesh)
mesh.draw().show()

filtered_layerstack = LayerStack(
    layers={
        k: LAYER_STACK.layers[k]
        for k in (
            "sin",
            "clad",
            "box",
        )
    }
)

resolutions = {}
resolutions["sin"] = {"resolution": 0.02, "distance": 2}
resolutions['clad'] = {"resolution": 0.2, "distance": 1}
resolutions['box'] = {"resolution": 0.2, "distance": 1}
c = gf.components.straight(cross_section=strip_SiN)
mesh = c.to_gmsh(
    type="uz",
    xsection_bounds=[(5, -5), (5, 5)],
    layer_stack=filtered_layerstack,
    filename=f"mesh.msh",
    resolutions=resolutions
)
mesh = from_meshio(mesh)
mesh.draw().show()

Mesh:

Epsilon Plot SiN:

Epsilon Plot Box:

When trying to create the epsilon plot for the cladding I get the following error:
Traceback (most recent call last): File "/home/st4eve/Desktop/purple/mask-design/purple_headers/tristan_sim/mode_solving.py", line 65, in <module> epsilon[basis0.get_dofs(elements=subdomain)] = n**2 File "/home/st4eve/miniconda3/envs/layout/lib/python3.10/site-packages/skfem/assembly/dofs.py", line 172, in __array__ return self.flatten() File "/home/st4eve/miniconda3/envs/layout/lib/python3.10/site-packages/skfem/assembly/dofs.py", line 94, in flatten (self.obj IndexError: index 9014 is out of bounds for axis 1 with size 9014

Hi everyone,

I was trying to get the overlap integral for both quasi-TE/TM modes as provided by the example but I am not sure if I did something wrong. Here is the result I got for the same example.

from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
from shapely import box
from skfem import Basis, ElementTriP0, ElementTriP1
from skfem.io import from_meshio
from tqdm import tqdm

from femwell.fiber import e_field_gaussian, overlap
from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict




wl=1.55
core_w,core_h=0.2,0.3
clad_w,clad_h=15,15
num_modes=2

core = box(-core_w/2, -core_h/2, core_w/2, core_h/2)
polygons = OrderedDict(core=core, clad=core.buffer(clad_w, resolution=12))

resolutions = dict(core={"resolution": 0.02, "distance": 0.5})
mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=5))
mesh.draw().show()

basis0 = Basis(mesh, ElementTriP0(), intorder=4)
epsilon = basis0.zeros().astype(complex)
epsilon[basis0.get_dofs(elements="core")] = 1.9963**2
epsilon[basis0.get_dofs(elements="clad")] = 1.444**2
basis0.plot(np.real(epsilon)**0.5, colorbar=True).show()

modes = compute_modes(basis0, epsilon, wavelength=wl, mu_r=1, num_modes=num_modes)

mfds = np.linspace(2, 20, 2**5)
TE_eff=np.zeros_like(mfds)
TM_eff=np.zeros_like(mfds)


#sort the modes as TE&TM
for i,mode in enumerate(modes):
    if mode.te_fraction >= 0.7:
        TE_mode=mode
        print(f"TE mode: {mode.te_fraction*100}%")
        modes[i].plot_intensity()
        plt.title("TE")
        plt.show()
        fig, axs = modes[i].plot(modes[i].E.real, direction="x")
            
    elif mode.te_fraction <= 0.3:
        
        TM_mode=mode
        print(f"TM mode: {100-mode.te_fraction*100}%")
        modes[i].plot_intensity()
        plt.title("TM")
        plt.show()
        fig, axs = modes[i].plot(modes[i].E.real, direction="x")
    else:
        print("Neither quasi-TE/TM modes")
        
for j,mfd in enumerate(tqdm(mfds)):
    basis_fiber = basis0.with_element(ElementTriP1())
    x_fiber = basis_fiber.project(
        lambda x: e_field_gaussian(np.sqrt(x[0] ** 2 + x[1] ** 2), 0, mfd / 2, 1, wavelength=wl),
        dtype=complex,
    )

    TE_eff[j] = overlap(
        basis_fiber, TE_mode.basis.interpolate(TE_mode.E)[0][0], basis_fiber.interpolate(x_fiber)
    )
    TM_eff[j] = overlap(
        basis_fiber, TM_mode.basis.interpolate(TM_mode.E)[0][1], basis_fiber.interpolate(x_fiber)
    )

plt.figure()
plt.plot(mfds, TE_eff,label="TE")
plt.plot(mfds, TM_eff,label="TM")
plt.xlabel("Mode field diameter / um")
plt.ylabel("Coupling efficiency")
plt.legend()
plt.tight_layout()
plt.show()

from collections import OrderedDict

from skfem import *
import shapely

from mesh import mesh_from_OrderedDict

height = 3
a = .330
b = .7
c = .2

box = shapely.box(0,-height,a,height)
structure = shapely.box(0,-b/2,a,b/2)
hole = shapely.box(0,-c/2,a/2,c/2)

resolutions = {
    'box': {'resolution':.1, 'distance':1}
}

mesh = mesh_from_OrderedDict(OrderedDict(
    hole=hole,
    structure=structure,
    box=box,
), resolutions=resolutions, filename='mesh.msh')

doesn't mesh :/

@simbilod

In electrostatic and thermal simulations it is fairly common to sweep a variable, say voltage or temperature or something else.

In a sweep, using the solution computed in the previous step can speed things up and sometimes even allow the simulation to converge when it would not if you initialized the solver from scratch. Is there a possibility to provide an initial solution to the solvers?

Could we allow the meshing to do periodic meshes? @simbilod

Luckily, adjoint methods for (most) eigenvalue problems tend to be rather trivial if you are optimizing the eigenvalue directly. So, for example, if you want to match the effective index of two modes by optimizing the dimensions of the waveguide, you can use classical gradient-descent methods to do so.

Since there's only a few parameters (eg width, height), it's reasonable to just use a derivative-free method. However, in this case, it's still easy to compute the gradient using an adjoint method. Specifically, you can leverage the Hellmann-Feynman theorem, which only requires the mode fields themselves (ie there's no adjoint solve needed).

If you wanted to also optimize a quantity that depends on the eigenvectors (like matching the group index of two modes) then there's an extra perturbation term involving an additional linear solve (though since you are assembling an FEM matrix, you can probably perform this linear solve really easily). This sort of optimization is especially useful for modulator design, for example.

The hardest part of all of this is computing the $A_u$ term, which is the derivative of your system matrix wrt the DOF. But that can be automated with autodiff if done right.

But you can imagine some interesting dispersion engineering, especially when wanting to engineer RF and optical modes simultaneously. Especially if you can start to introduce some interesting DOF.

(cc @joelslaby, @joamatab)

It would be nice to be able to plot only Ex, Ey, or Ez, rather than all three as displayed by the mode.show() method.

I looked at the plot_mode function, and I see that plotting all three is hard coded. Adding an optional parameter to the function would be an easy addition.

https://github.com/HelgeGehring/femwell/blob/main/femwell/maxwell/waveguide.py

https://github.com/HelgeGehring/femwell/blob/main/docs/references.bib is a bit messy I think. Maybe good to run it through https://github.com/FlamingTempura/bibtex-tidy and/or https://github.com/Pezmc/BibLatex-Check

Does not seem to have slepc

Are anisotropic materials supported?
I couldn't find anything in the documentation or the examples.

It would be great to add the boundary conditions for TE-modes and TM-modes at the outer boundaries of the domain to the math-docs

I am running into issues installing slepc4py on Windows which requires petsc.
Petsc seems to be requiring a UNIX version of cygwin.

https://petsc.org/main/install/windows/

I am not sure how to pass this step, can you help me with this?

See description in #141 and for example notebook

Can you release a new version?

femwell depends on scikit-fem@ git+https://github.com/requests/requests.git

gdsfactory/gdsfactory#1032

Would it make sense to have dedicated runners for this project?

https://github.com/antmicro/runner should make it easy to deploy in a existing Google Cloud environment.

It would be great to have a feature to mode solve in a subdomain, or only keep the modes that have X amount of overlap with some spatial region automatically (an example might be enough)

Would be especially useful to deal with PML modes

@HelgeGehring @mdecea

Would be awesome if the meshing would support 3D meshing! @simbilod

Hi guys,
The mode class makes accessing results much easier now, but I am still puzzled by the fields and their bases.
I noticed in waveguide.py that the transversal and the normal bases differ (e.g.: line 398), so that Ex and Ey are sampled on a different grid with respect to Ez.
Is there a specific reason for this? In FDE (MPB, Lumerical, etc.) there is only one grid, so I would argue that all components could in principle be associated to one basis only.
Am I missing anything here? Thanks in advance!

How can we add python3.7 support?

that way we can use it in binder

I am calculating the dispersion for different geometries. And mainly use loops for now. Is it possible to implement the calculations on arrays to avoid the loops? I think i need to redefine the mesh for every geometry but at least the wavelength sweep?
https://helgegehring.github.io/femwell/photonics/examples/vary_wavelength.html#

Could we allow non-triangular meshes, e.g. quadrilateral meshes? @simbilod

See #68

@simbilod

At the moment the units in the labels are a bit mixed in their format, e.g. for seconds:

• t [s] (which I mostly used up to now)
• t (s)
• t / s

The last one is according to https://academia.stackexchange.com/questions/18357/are-there-any-guidelines-for-labeling-axes-in-plots-graphs the preferred formatting according to the International System of Units

should we just stick to the last one?

@simbilod

meshing code originally contained in the repo has been refactored to be pure GMSH and extended to 3D, see https://github.com/simbilod/meshwell

TODO

• rewrite examples to use the new API (very similar)
• make meshwell an (optional?) dependency
• remove legacy code

Hello @HelgeGehring
Thank you for this great optical mode calculating tool.
Since I'm trying to find modes of an anisotropic dielectric waveguide, I would like to install a Julia version. I've already install through pip, but I'm not able to run the Julia example code. Am I supposed to install a Julia package of femwell? If so, how can I install it? (It seems that "add Femwell" doesn't work.)

Thank you.

See docs https://helgegehring.github.io/femwell/photonics/examples/selecting_modes.html, nitride mode is not being sorted correctly

Pickling Modes currently leads to ~100MB objects on disk

What would be a best practice to save/load modes to disk?

In gdsfactory with the old femwell API I had something that saved the mesh, element, an field values, and would reconstruct the basis at runtime, which led to much smaller files (~1MB or less instead of 100)

@mdecea and I just had a discussion and we think it might be good if compute_modes would return a list of objects instead of the eigenvalues.
Attributes:

• effective refractive index
• propagation constant
• wavelength
• frequency
• E
• H
• power
• energy_current_density
• basis
• epsilon
• tefrac
• tmfrac

And methods like:

• overlap_with(another_mode)
• scalar_product_with(another_mode)
• coupling_coefficient_with(another_mode)
• confinement_factor(within_elements)
• plot

@simbilod @mdecea as this is probably the largest refactoring up to now, would be great to have some feedback before I start doing this

The next step might be an object which replaces the list with all these modes and offers functions to sort the modes by

• effective refractive index - real part
• effective refractive index - imag part
• power in a certain area
  and of course the attributes
• epsilon
• wavelength / ...

that might be a reason for version 0.1

issues with scipy backend

when running waveguide_modes.py

AttributeError                            Traceback (most recent call last)
Cell In[1], line 129
    124 @LinearForm
    125 def unit_load(v, w):
    126     return v
--> 129 M = unit_load.assemble(basis.with_elements("core"))
    131 times = np.array([dt * i for i in range(steps + 1)])
    132 plt.xlabel("Time [/](https://jupytext+.vscode-resource.vscode-cdn.net/) us")

AttributeError: 'CellBasis' object has no attribute 'with_elements'

Hello,

Just as @elizaleung830 I'm trying to calculate the effective area, but for Spontaneous Four-wave Mixing (SFWM), which is a third-order non-linear effect and the Aeff needs to take into account the interaction between the electric field distributions of three modes (we consider a degenerate pump, which means there are two identical pumps).

The goal is to replicate the results found in the Ansys Lumerical tutorial on Spontaneous Four-Wave Mixing (SFWM) in a Microring Resonator Photon Source, as documented at https://optics.ansys.com/hc/en-us/articles/15100783091731, following the methodology outlined in the publication https://doi.org/10.1070/QEL16511:

Where $u_p$, $u_s$, $u_i$ are the mode function describing the transverse spatial distribution of the field and normalised $\int|u(x, y)|^{2} \mathrm{~d} x \mathrm{~d} y=1$.

So, for a mode in femwell, does mode.E correspond to the unnormalized electric field distribution? I added normalization in femwell.maxwell.waveguide and calculated the method for Aeff, but I obtained a rather strange result.

Function added (for testing I unfreeze the class to add a E0 for the normalized electric field):

def calculate_sfwm_Aeff(
    basis: Basis,
    mode_p,
    mode_s,
    mode_i,
) -> np.complex64:
    """
    Calculates the overlap integral for SFWM process by considering the interactions
    between pump, signal, and idler modes in the xy plane.

    Args:
        basis (Basis): Common basis used for all modes in the same geometric structure.
        mode_p, mode_s, mode_i: Mode instances for pump, signal, and idler, respectively.

    Returns:
        np.complex64: The overlap integral result for the SFWM process.
    """
    def normalize_mode(mode):
        @Functional
        def E2(w):
            return (w["E"][0] * np.conj(w["E"][0]))


        E_xy, _ = mode.basis.interpolate(mode.E) #xy, z
        E_squared_integral = E2.assemble(mode.basis, E=E_xy)
        normalization_factor = 1 / np.sqrt(E_squared_integral)
        mode.E0 = mode.E * normalization_factor
    
    normalize_mode(mode_p)
    normalize_mode(mode_s)
    normalize_mode(mode_i)

    #Normalization needed for uv(x,y)!
    E_p, _ = mode_p.basis.interpolate(mode_p.E0) #xy, z
    E_s, _ = mode_s.basis.interpolate(mode_s.E0)
    E_i, _ = mode_i.basis.interpolate(mode_i.E0)
    
    @Functional(dtype=np.complex64)
    def sfwm_overlap(w):
        overlap_SFWM = w["E_p"][0] * w["E_p"][0] * np.conj(w["E_s"][0]) * np.conj(w["E_i"][0])
        return (overlap_SFWM)

    # Assemble the integral over the basis to compute the overlap
    overlap_result = sfwm_overlap.assemble(
        basis,
        E_p=E_p,
        E_s=E_s,
        E_i=E_i,
    )

    return  1/overlap_result

And for the main code:

from collections import OrderedDict
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import speed_of_light
from shapely.geometry import box
from shapely.ops import clip_by_rect
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio
from tqdm import tqdm
from femwell.maxwell.waveguide import compute_modes
from femwell.maxwell.waveguide import calculate_sfwm_Aeff
from femwell.mesh import mesh_from_OrderedDict
from scipy.constants import c, epsilon_0 ,pi

def n_X(wavelength):
    x = wavelength
    return (1 
            + 2.19244563 / (1 - (0.20746607 / x) ** 2)
            + 0.13435116 / (1 - (0.3985835 / x) ** 2)
            + 2.20997784 / (1 - (0.20747044 / x) ** 2)
    ) ** 0.5

# Box
def n_silicon_dioxide(wavelength):
    x = wavelength
    return (
        1
        + 0.6961663 / (1 - (0.0684043 / x) ** 2)
        + 0.4079426 / (1 - (0.1162414 / x) ** 2)
        + 0.8974794 / (1 - (9.896161 / x) ** 2)
    ) ** 0.5

Clad=1


core = box(0, 0, 0.5, 0.39)
polygons = OrderedDict(
    core=core,
    box=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, -np.inf, np.inf, 0),
    clad=clip_by_rect(core.buffer(1.5, resolution=4), -np.inf, 0, np.inf, np.inf),
)

resolutions = {"core": {"resolution": 0.025, "distance": 2.}}

mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6))

num_modes = 2

lambda_p0 = 1.4
lambda_s0 = 1.097
lambda_i0 = 1.686

basis0 = Basis(mesh, ElementTriP0())

epsilon_p = basis0.zeros(dtype=complex)
epsilon_s = basis0.zeros(dtype=complex)
epsilon_i = basis0.zeros(dtype=complex)


for wavelength, epsilon in zip([lambda_p0, lambda_s0, lambda_i0], [epsilon_p, epsilon_s, epsilon_i]):
    for subdomain, n_func in {
        "core": n_X,
        "box": n_silicon_dioxide,
        "clad": lambda _: Clad,
    }.items():
        n = n_func(wavelength)
        epsilon[basis0.get_dofs(elements=subdomain)] = n**2


modes_p = compute_modes(basis0, epsilon_p, wavelength=lambda_p0, num_modes=num_modes, order=1)
modes_s = compute_modes(basis0, epsilon_s, wavelength=lambda_s0, num_modes=num_modes, order=1)
modes_i = compute_modes(basis0, epsilon_i, wavelength=lambda_i0, num_modes=num_modes, order=1)


mode_p = max(modes_p, key=lambda mode: mode.te_fraction)
mode_s = max(modes_s, key=lambda mode: mode.te_fraction)
mode_i = max(modes_i, key=lambda mode: mode.te_fraction)


A_eff = calculate_sfwm_Aeff(basis0, mode_p, mode_s, mode_i)
print(A_eff)

chi_3 = 5e-21  # m^2/V^2  #7e-20?
lambda_p0_m = lambda_p0 * 1e-6  #to m
n_p0 = np.real(mode_p.n_eff)
A_eff_m2 = A_eff * 1e-12  #to m^2

omega_p0 = 2 * pi * c / lambda_p0_m

gamma = (3 * chi_3 * omega_p0) / (4 * epsilon_0 * c**2 * n_p0**2 * A_eff_m2)

print("gamma:",gamma)

After confirming that the electric field distribution was normalized, I obtained a complex number with a negative real part for

?According to the comparison with Lumerical's results, the Aeff should be on the order of $0.469 \ \text{um}^2$.

For FEFD, we had "issues" due to source directionality

Is this even a problem? In https://nbviewer.org/github/fancompute/workshop-invdesign/blob/master/01_First_simulation.ipynb they just compute the Poynting vector somewhere away from the source.

Could we simply compare the amount of propagating energy and/or the fields close to the source vs at the target outputs modes/positions to get the S-parameters?

@HelgeGehring

Does it support the calculation of effective mode area or nonlinear coefficent of waveguide? I find the equation effective area in the maxwell.md but i can't find it in the example.
Thanks in advance.