Since there is already noticeable interest towards systems with several waveguides, would be nice to have native multi-waveguide geometry in TetraX, similarly to multilayer geometry. The simplest example would be directional coupler -like geometry, i.e. rectangular waveguides with the same thickness and different widths / spacings, coupled by dipole field only.
Hi Anton,
first of all, welcome to the TetraX forum and thank you for submitting your feature request. This is something easy to implement. We will do this shortly but it might not be lead to a new patch right away. In meantime, we will provide you with an example code to generate the mesh in a notebook within the next few days.
All the best,
Lukas
Hi Anton,
Here comes a general definition of a multilayer stack, with arbitrary widths and thicknesses and spacings between them.
Please make some more tests and let us know if this is what you expected to have.
Cheers,
Attila
import tetrax as tx
import numpy as np
import pygmsh
from typing import List, Union
from tetrax.sample.mesh.sample_mesh import GeometryType
from tetrax.sample.mesh.sample_mesh import MeshioMesh
def multilayer( # noqa: PLR0913
widths: List[float],
thicknesses: List[float],
spacings: List[float],
cell_sizes_width: Union[float, List[float]],
cell_sizes_thickness: Union[float, List[float]],
x0: float = 0,
y0: float = 0,
) -> MeshioMesh:
r"""
Cross-section magnetic multi-layer stack with arbitrary but finite widths.
Parameters
----------
widths : list(float)
Widths :math:`[w_1, w_2, w_3, ...]` of the magnetic layers (ordered with increasing :math:`y` coordinate).
List must be the same size as thicknesses and one element longer than ``spacings``.
thicknesses : list(float)
Thicknesses :math:`[d_1, d_2, d_3, ...]` of the magnetic layers (ordered with increasing :math:`y` coordinate).
List must be one element longer than ``spacings``.
spacings : list(float)
Spacings :math:`[s_1, s_2, s_3, ...]` between the magnetic layers, which is are the thicknesses of the non-magnetic interlayers
(ordered with increasing :math:`y` coordinate). List must be one element shorter than ``thicknesses``.
cell_size_width : list(float) or float
Characteristic lengths along the transversal direction the different magnetic layers.
If list, has to be same length as ``widths``. If not a list it will be same for all layers.
cell_sizes_thickness : list(float) or float
Characteristic lengths along the normal direction the different magnetic layers.
If list, has to be same length as ``thicknesses``. If
not a list it will be same for all layers.
x0 : float
x coordinate of the center (Default is 0).
y0 : float
y coordinate of the center (Default is 0).
Returns
-------
mesh : MeshioMesh
See Also
--------
rectangular, bilayer
"""
if len(widths) != len(thicknesses):
print(
f"Number of layer widths ({len(widths)}) and number of layer thicknesses ({len(thicknesses)}) do "
f"not match (has to be the same)."
)
cell_sizes_width = (
[cell_sizes_width for w in widths]
if not isinstance(cell_sizes_width, list)
else cell_sizes_width
)
cell_sizes_thickness = (
[cell_sizes_thickness for d in thicknesses]
if not isinstance(cell_sizes_thickness, list)
else cell_sizes_thickness
)
if len(spacings) != len(thicknesses) - 1:
print(
f"Number of layer thicknesses ({len(thicknesses)}) and number of spacings ({len(spacings)}) do "
f"not match (has to be off by one)."
)
if len(cell_sizes_width) != len(widths):
print(
f"Number of layer widths ({len(widths)}) and number of characteristic lengths ({len(cell_sizes_widths)}) do "
f"not match (has to be the same)."
)
if len(cell_sizes_thickness) != len(thicknesses):
print(
f"Number of layer thicknesses ({len(thicknesses)}) and number of characteristic lengths ({len(cell_sizes_thickness)}) do "
f"not match (has to be the same)."
)
d_tot = sum(spacings) + sum(thicknesses)
with pygmsh.geo.Geometry() as geom:
y1_i = -d_tot / 2 + y0
for width, d, s, lcw, lct in zip(widths, thicknesses, [*spacings, 0], cell_sizes_width, cell_sizes_thickness, strict=False):
# ny_i = int(d / lc)
# p1 = geom.add_point([0.0, y1_i, 0.0], mesh_size=lc)
# geom.extrude(p1, translation_axis=(0.0, d, 0.0), num_layers=ny_i)
nx = int(width / lcw)
ny = int(d / lct)
p = geom.add_point([-width / 2 + x0, y1_i, 0.0], mesh_size=lcw)
_, line, _ = geom.extrude(p, translation_axis=[width, 0, 0], num_layers=nx)
geom.extrude(line, translation_axis=[0, d, 0], num_layers=max(2, ny))
y1_i += d + s
mesh = geom.generate_mesh()
mesh.geometry_type = GeometryType.WAVEGUIDE
if len(spacings) == 0:
return mesh
xyz = mesh.points
upper = []
lower = []
y1_i = -d_tot / 2 + y0
for d, s in zip(thicknesses[:-1], spacings, strict=False):
y1_i += d
upper += list(np.where(xyz[:, 1] == y1_i)[0])
y1_i += s
lower += list(np.where(xyz[:, 1] == y1_i)[0])
mesh.interlayer_surface_nodes = np.concatenate((upper, lower))
mesh.interlayer_surface_partners = np.concatenate((lower, upper))
return mesh
# test with two waveguides of different widths
widths = [200, 50]
thicknesses = [10,10]
spacing = [5]
cell_size_width = 5
cell_size_thickness = 2
mlayer_msh = multilayer(widths,thicknesses, spacing,cell_size_width, cell_size_thickness)
mlayer = tx.Sample(mlayer_msh)
mlayer.show()
# let's apply some fields and relax it
mlayer.mag = (1,0,1)
mlayer.Bext = (100e-3,0,0)
relax = tx.experiments.relax(mlayer)
mlayer.show()
# computation of the spectrum
spectrum = tx.experiments.eigenmodes(mlayer)
# plotting the first 4 modes
spectrum.plot(n=[0,3], fscale='G', kscale='u')Hi Attila,
Thanks for the code! I was able to reproduce something similar to the Py-Py vertical directional coupler, but unfortunately not the original YIG-YIG horizontal directional coupler. For the YIG-YIG coupler, the frequencies seem to be about 1.5 times higher than in the paper. Here are the parameters I tried using:
w1 = 85
w2 = 85
t1 = 350
t2 = 350
s = 320
widths = [w1, w2]
thicknesses = [t1,t2]
spacing = [s]
cell_size_width = 5
cell_size_thickness = 10
mlayer_msh = multilayer(widths,thicknesses, spacing, cell_size_width, cell_size_thickness)
mlayer = tx.Sample(mlayer_msh)
# mlayer.show()
# let's apply some fields and relax it
BextVal = 56
BextAxis = 2
magList = [0,0,0]
bextList = [0,0,0]
magList[BextAxis] = 1
bextList[BextAxis] = BextVal*1e-3
mlayer.mag = magList
mlayer.Bext = bextList
wgmat = "YIG"
mlayer.material = tx.materials.yig
relax = tx.experiments.relax(mlayer)
# mlayer.show()
# computation of the spectrum
kval = 15e6
spectrum = tx.experiments.eigenmodes(mlayer, num_modes=4, kmin=-kval, kmax=kval, num_k=201)
Hi Anton,
Thanks for the tests!
I’m convinced that if you would have used the material parameters given in the manuscript you linked in (Micromagnetic simulations section), you would have gotten a much better agreement.
The YIG material parameters implemented into TetraX are a little different from those in the paper.
Cheers,
Attila
Indeed, it seems to be much closer after taking the values from the paper (though there is still some difference, maybe it could be attributed to magnetic field or anisotropy differences). I didn’t expect 1.5× difference from the minor differences in material parameters, to be honest 
Here are the plots, by the way:
Happy to hear! Indeed, sometimes small differences in material parameters can result in large frequency shifts. In case you want a one to one correspondence (or as close as possible to it), you might need to get the mumax^3 script and check carefully all material parameters, including gamma, anisotropies (if there are any considered), discretization and so on. I often face this problem, that in the manuscript - for whatever reasons - people publish different parameters to what they actually used.
Cheers,
Attila