Skip to content

Commit

Permalink
Update ModeSolvers-Tutorial.md
Browse files Browse the repository at this point in the history 
Update PCF tutorial to take advantage of the functions neff_approx_PCF and add_cylindrical_PML
  • Loading branch information
@ovanvincq
ovanvincq authored Jan 12, 2024
1 parent 56d647f commit 6aa61a5
Showing 1 changed file with 21 additions and 124 deletions.
145 changes: 21 additions & 124 deletions docs/src/ModeSolvers-Tutorial.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,7 +168,7 @@ save("FEM1_Pz.png",fig); nothing #hide
![Pz for FM computed with FEM](FEM1_Pz.png)

## FEM: Photonic Crystal Fiber
This example is more complex since, in a PCF, the modes are not guided modes but leaky modes so that the computation requires a PML. The fiber is constituted of three rings of air hole (n=1) inserted in silica (n=1.45). The pitch is 2 µm, the hole diameter is 1.5 µm and the PML begin at 8 µm from the core and its thickness is 2 µm.
This example is more complex since, in a PCF, the modes are not guided modes but leaky modes so that the computation requires a PML. The fiber is constituted of three rings of air hole (n=1) inserted in silica (n=1.45). The pitch is 2 µm, the hole diameter is 1.5 µm and the PML begins at 8 µm from the fiber center and its thickness is 2 µm.
![PCF mesh](./assets/PCF.png)
First, the mesh is loaded:
```@example 4
Expand All @@ -180,136 +180,33 @@ using GridapMakie
using GLMakie
model = GmshDiscreteModel("../../models/PCF.msh");
```
Then we define the functions useful in the construction of the permittivity and permeability tensors:
Then we define the permittivity function and add a PML to obtain the final permittivity and permeability tensors:
```@example 4
const alpha=10.0;
const r_pml=8.0;
const d_pml=2.0;
const r_hole=0.75;
const eps3=1.45^2;
const Pitch=2;
x1,y1=ring(1);
x2,y2=ring(2);
x3,y3=ring(3);
xc=[x1;x2;x3]*Pitch;
yc=[y1;y2;y3]*Pitch;
function eps_xx(x)
r=hypot(x[1],x[2]);
if r<=r_pml
for i=1:length(xc)
rc=hypot(x[1]-xc[i],x[2]-yc[i]);
if rc<r_hole
return ComplexF64(1.0);
end
function eps_PCF(x)
Pitch=2;
r_hole=0.75;
x1,y1=ring(1);
x2,y2=ring(2);
x3,y3=ring(3);
xc=[x1;x2;x3]*Pitch;
yc=[y1;y2;y3]*Pitch;
for i=1:length(xc)
rc=hypot(x[1]-xc[i],x[2]-yc[i]);
if rc<r_hole
return 1.0;
end
return ComplexF64(eps3);
else
phi=atan(x[2],x[1]);
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return eps3*(rt/(r*sr)*(cos(phi))^2+(r*sr/rt)*(sin(phi))^2);
end
return 1.45^2;
end
function eps_yy(x)
r=hypot(x[1],x[2]);
if r<=r_pml
for i=1:length(xc)
rc=hypot(x[1]-xc[i],x[2]-yc[i]);
if rc<r_hole
return ComplexF64(1.0);
end
end
return ComplexF64(eps3);
else
phi=atan(x[2],x[1]);
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return eps3*(rt/(r*sr)*(sin(phi))^2+(r*sr/rt)*(cos(phi))^2);
end
end
function eps_zz(x)
r=hypot(x[1],x[2]);
if r<=r_pml
for i=1:length(xc)
rc=hypot(x[1]-xc[i],x[2]-yc[i]);
if rc<r_hole
return ComplexF64(1.0);
end
end
return ComplexF64(eps3);
else
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return eps3*(rt/r)*sr;
end
end
function eps_xy(x)
r=hypot(x[1],x[2]);
if r<=r_pml
return ComplexF64(0.0);
else
phi=atan(x[2],x[1]);
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return eps3*sin(phi)*cos(phi)*(rt/(r*sr)-r*sr/rt);
end
end
function mu_xx(x)
r=hypot(x[1],x[2]);
if r<=r_pml
return ComplexF64(1.0);
else
phi=atan(x[2],x[1]);
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return (rt/(r*sr)*(cos(phi))^2+(r*sr/rt)*(sin(phi))^2);
end
end
function mu_yy(x)
r=hypot(x[1],x[2]);
if r<=r_pml
return ComplexF64(1.0);
else
phi=atan(x[2],x[1]);
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return (rt/(r*sr)*(sin(phi))^2+(r*sr/rt)*(cos(phi))^2);
end
end
function mu_zz(x)
r=hypot(x[1],x[2]);
if r<=r_pml
return ComplexF64(1.0);
else
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return (rt/r)*sr;
end
end
function mu_xy(x)
r=hypot(x[1],x[2]);
if r<=r_pml
return ComplexF64(0.0);
else
phi=atan(x[2],x[1]);
rt=r-im*alpha/3.0*(r-r_pml)^3/(d_pml)^2;
sr=1.0-im*alpha*(r-r_pml)^2/(d_pml)^2;
return sin(phi)*cos(phi)*(rt/(r*sr)-r*sr/rt);
end
end
function zf(x)
return ComplexF64(0.0)
end
epsilon=tensor3(eps_xx,eps_xy,zf,eps_xy,eps_yy,zf,zf,zf,eps_zz);
mu=tensor3(mu_xx,mu_xy,zf,mu_xy,mu_yy,zf,zf,zf,mu_zz);
epsilon=add_cylindrical_PML(eps_PCF,8,2,10);
mu=add_cylindrical_PML(x->1.0,8,2,10);
```
Then we can compute ten modes whose effectif index is close to 1.41 at the wavelength of 1.3 µm and plot the z-component of the Poynting Vector of the interesting mode:
Then we can compute ten modes whose effective indices are close to the approximate value calculated for the fundamental mode at the wavelength of 1.3 µm:
```@example 4
m=FEM(1.3,10,epsilon,mu,model,1.41,field=true,order=2,solver=:MUMPS)
neff_approx=approx_neff_PCF(1.3,1.5,2);
m=FEM(1.3,10,epsilon,mu,model,neff_approx,field=true,order=2,solver=:MUMPS)
```
The last two modes are fundamental modes. We can compute and plot the z-component of the Poynting vector of the last mode:
```@example 4
Px,Py,Pz=PoyntingVector(m[end]);
fig,ax,plot_obj=GLMakie.plot(m[end].Ω,Pz,axis=(aspect=DataAspect(),),colormap=:jet)
Expand Down

0 comments on commit 6aa61a5

Please sign in to comment.