Hello everyone,
I need some guidance on how to achieve proper results with this simulation. I want to simulate the light propagation through a medium with an inclusion. I assume sourceless, non-magnetic, time-harmonic conditions and want to include birefringence, leading to the Maxwell equations. The inclusion in the considered inhomogeneous material leads to a higher refractive index in its center and thus a tensorial quantity for the permittivity of the whole medium. For simplicity and starting, I just simulate a block of material with a dummy permittivity tensor. However, the light wave is not propagating through the medium at all and is somehow damped.
The differential equation I’m using is
With that, I’m assuming trivial Neumann boundary conditions on the boundary. However, the surface integral is implemented, though the factor f is set to zero, so its removed as in the variational formulation as above. I’m also applying periodic boundary conditions to neglect effects of the boundary on the light propagation. My implemented Code looks as follows:
load "msh3"
real wavelength = 1;
real k0=2*pi/wavelength;
real thickness = 500;
real height=thickness*3;
int n=5;
int[int] labellist= [94,95,96,97,98,99];
mesh3 Th = cube(
height/thickness*n,
height/thickness*n,
thickness/thickness*n,
//[x*height,y*height,z*thickness],
[x*height-height/2,y*height-height/2,z*thickness-thickness/2],
label=labellist,
flags=3
);
fespace Vh(Th, P1, periodic=[[94,z,x],[96,z,x],[95,z,y],[97,z,y]]);
Vh<complex> Ex,Ey,Ez,vx,vy,vz;
func real epsxx(real a, real b, real c) {(2.0)^2;}; // ε_xx
func real epsyy(real a, real b, real c) {(2.0)^2;}; // ε_yy
func real epszz(real a, real b, real c) {(2.0)^2;}; // ε_zz
func real epsxy(real a, real b, real c) {(1.5)^2;}; // ε_xy
func real epsxz(real a, real b, real c) {(1.5)^2;}; // ε_xz
func real epsyz(real a, real b, real c) {(1.5)^2;}; // ε_yz
func real epsyx(real a, real b, real c) {epsxy;}; // ε_yx
func real epszx(real a, real b, real c) {epsxz;}; // ε_zx
func real epszy(real a, real b, real c) {epsyz;}; // ε_zy
macro epsDotEDotv(E1,E2,E3,v1,v2,v3,a,b,c) (
epsxx(a,b,c)*E1*v1 + epsyx(a,b,c)*E2*v1 + epszx(a,b,c)*E3*v1 +
epsxy(a,b,c)*E1*v2 + epsyy(a,b,c)*E2*v2 + epszy(a,b,c)*E3*v2 +
epsxz(a,b,c)*E1*v3 + epsyz(a,b,c)*E2*v3 + epszz(a,b,c)*E3*v3
) //
real f = 0;
problem maxwell([Ex,Ey,Ez],[vx,vy,vz]) =
int3d(Th)(
dx(Ex)*dx(vx)+dy(Ex)*dy(vx)+dz(Ex)*dz(vx)+
dx(Ey)*dx(vy)+dy(Ey)*dy(vy)+dz(Ey)*dz(vy)+
dx(Ez)*dx(vz)+dy(Ez)*dy(vz)+dz(Ez)*dz(vz)
)
- int3d(Th)(
k0^2*epsDotEDotv(Ex,Ey,Ez,vx,vy,vz,x,y,z)
)
- int2d(Th)(
f*vx+f*vy+f*vz
)
+ on(98,Ey=1)
;
maxwell;
Vh Eim = imag(Ey);
Vh EInt = imag(Ex)^2+real(Ex)^2+imag(Ey)^2+real(Ey)^2+imag(Ez)^2+real(Ez)^2;
Vh EIntlog = log(imag(Ey)^2+real(Ey)^2);
plot(Ex,cmm="Ex",wait=1,value=0,fill=1,nbiso=65);
plot(Ey,cmm="Ey",wait=1,value=0,fill=1,nbiso=65);
plot(Ez,cmm="Ez",wait=1,value=0,fill=1,nbiso=65);
plot(Eim,cmm="Im{Ey}",wait=1,value=0,fill=1,nbiso=65);
plot(EInt,cmm="|E|",wait=1,value=0,fill=1,nbiso=65);
with the possibility of location depending epsxx, epsyy,… for later investigations. Hoewever, its not even possible to see the waves propagating trough the medium. With the boundary condition Ey=1 on the plane, the light should travel through the thickness of the block. With a wave length of 1000, I can see some light waves propagating throug the medium, but that corresponds to unphysically thin material. Am I missing something or are there mistakes in my equations?
Another part I can’t wrap my head around is the possible extension for the phase evaluation. Here, I would choose an imaginary ansatz for Ex, Ey and Ez as follows: Ex=Ex exp(i phi), Ey=Ey exp(i phi), Ez=Ez exp(i phi) with and equal phase. Hence, I could have the dependencies on four real variables, Ex, Ey, Ez and phi but would require and complex solver. Is this somehow possible? My approach would look like following
Vh Exc,Eyc,Ezc,vxc,vyc,vzc;
//Vh<complex> Exc,Eyc,Ezc,vxc,vyc,vzc;
problem maxwellIM([Exc,Eyc,Ezc,phi],[vxc,vyc,vzc,psi]) =
int3d(Th)(
((dx(Exc)+1i*Exc*dx(phi))*(dx(vxc)+1i*vxc*dx(psi)) + (dx(Eyc)+1i*Eyc*dx(phi))*(dx(vyc)+1i*vyc*dx(psi)) +
(dx(Ezc)+1i*Ezc*dx(phi))*(dx(vzc)+1i*vzc*dx(psi)) + (dy(Exc)+1i*Exc*dy(phi))*(dy(vxc)+1i*vxc*dy(psi)) +
(dy(Eyc)+1i*Eyc*dy(phi))*(dy(vyc)+1i*vyc*dy(psi)) + (dy(Ezc)+1i*Ezc*dy(phi))*(dy(vzc)+1i*vzc*dy(psi)) +
(dz(Exc)+1i*Exc*dz(phi))*(dz(vxc)+1i*vxc*dz(psi)) + (dz(Eyc)+1i*Eyc*dz(phi))*(dz(vyc)+1i*vyc*dz(psi)) +
(dz(Ezc)+1i*Ezc*dz(phi))*(dz(vzc)+1i*vzc*dz(psi)))*exp(1i*phi)*exp(1i*psi)
)
- int3d(Th)(
k0^2 * epsDotEDotv(Exc,Eyc,Ezc,vxc,vyc,vzc) *exp(1i*phi)*exp(1i*psi)
)
- int2d(Th)(
(f*vxc + f*vyc + f*vzc)*exp(1i*psi)
)
+ on(98, Eyc=1,phi=0)
;
maxwellIM;
Vh Eimc = imag(Eyc);
Vh EIntc = imag(Exc)^2+real(Exc)^2+imag(Eyc)^2+real(Eyc)^2+imag(Ezc)^2+real(Ezc)^2;
Vh EIntlogc = log(imag(Eyc)^2+real(Eyc)^2);
plot(Exc,cmm="Ex",wait=1,value=0,fill=1,nbiso=65);
plot(Eyc,cmm="Ey",wait=1,value=0,fill=1,nbiso=65);
plot(Ezc,cmm="Ez",wait=1,value=0,fill=1,nbiso=65);
plot(Eimc,cmm="Im{Ey}",wait=1,value=0,fill=1,nbiso=65);
plot(EIntc,cmm="|E|",wait=1,value=0,fill=1,nbiso=65);
but is yielding an error and not allowing the product with imaginary i. Independent of how Exc,Eyc,… are defined, as real or complex.
Thanks in advance!
Edit: Typo, imaginary approach
