Frequency domain Maxwell equations with PML boundary conditions

Question

I'm trying to solve a full-vectorial wave equation for an arbitrarily shaped wave guide, by using NDSolve and perfectly matched layer (PML) conditions.

The PML conditions can be stated as a coordinate transformation of the form $\text{$\partial $/$\partial $x $\to $ $\alpha $x(x) $\partial $/$\partial $x}$ and so on. As the functionality of the Mathematica built-in Curl doesn't extend so far I set up a function that acts as intended:

generalizedCurl3D[coordTransfVector_,applicationVector_,coordNamesVector_]:={
  coordTransfVector[[2]] D[applicationVector[[3]],coordNamesVector[[2]]] -
   coordTransfVector[[3]] D[applicationVector[[2]],coordNamesVector[[3]]],
  coordTransfVector[[3]] D[applicationVector[[1]],coordNamesVector[[3]]] - 
   coordTransfVector[[1]] D[applicationVector[[3]],coordNamesVector[[1]]],
  coordTransfVector[[1]] D[applicationVector[[2]],coordNamesVector[[1]]] - 
   coordTransfVector[[2]] D[applicationVector[[1]],coordNamesVector[[2]]]
}

If the coordinate transformation is set to unity $\text{$\{$1,1,1$\}$}$ the function is identical to the normal Curl, as the comparison yields True:

generalizedCurl3D[{1, 1, 1}, {ψx[x, y, z], ψy[x, y, z], ψz[x, y, z]}, 
    {x, y, z}] == {-Derivative[0, 0, 1][ψy][x, y, z] + 
    Derivative[0, 1, 0][ψz][x, y, z], Derivative[0, 0, 1][ψx][x, y, z] - 
    Derivative[1, 0, 0][ψz][x, y, z], -Derivative[0, 1, 0][ψx][x, y, z] + 
    Derivative[1, 0, 0][ψy][x, y, z]}

Now I derive the wave equation I need to solve. In the following I will consider only the magnetic field $\text{$\psi $[x,y,z]}$ because it is continuous at the boundary of the wave guide. The refractive index profile $n$ depends on $x$ and $y$ ($n=n[x,y]$) and is constant along the $z$ direction. The functions for PML will be set up later, but at this moment it is only important that the coordinate transformation functions $\text{$\alpha $x}$ and $\text{$\alpha $y}$ depend only on $x$ and $y$ ($\text{$\alpha $x[x]}$ and $\text{$\alpha $y[y]}$) while $\text{$\alpha $z}=1$.

I start with the Helmholtz equation, assuming that the time separation ansatz works. This way I only have to assume that the field will have fast oscillations in z direction with the effective refractive index $n0$ and the wave number $k$ ($e^{-i \cdot k \cdot \text{n0}\cdot z}$). The rest should be only slow oscillations of the field components in $x$ and $y$ direction. I also tend to set the $\text{$\psi $z}$ component of the $\psi$-vector to zero, because probably the cross-talk of $x$,$y$ components to the $z$ field component should be negligible.

Edit: In case you try to check my calculation, please make sure that the 3rd component of $\psi $ is defined to be zero. To keep it general for 3dimensional analysis (with easy "switching on or off") I define it nevertheless, but multiply it with 0.

cTV = {αx[x], αy[y], 1};
cNV = {x, y, z};
ψ = {ψx[x, y, z] E^(-I k n0 z), ψy[x, y, z] E^(-I k n0 z), 0 ψz[x, y, z] E^(-I k n0 z)};

Plugging now the still analytical expressions into the Helmholtz equation

($\text{$\nabla \times \nabla \times \psi $== - }\mu \text{$\epsilon $ }\partial ^2\left/\partial t^2\right.\text{ $\psi $ = + }k^2\text{ $\psi $}$)

and dividing all by the fast changing $z$-term ($e^{-i \cdot k \cdot \text{n0}\cdot z}$) gives me the following :

eqs = FullSimplify[(generalizedCurl3D[cTV, 1/n[x, y]^2 generalizedCurl3D[cTV, ψ, cNV], cNV]
    - k^2 ψ ) 1/E^(-I k n0 z)]

Now I replace the analytical expressions by numerical values and functions and define the size of the boundary and the PML-layer (all in SI units). (Note that $n0$ is effectively a propagation constant that has to be found by trial and error: if it is chosen wrong, then the field should spread outside the waveguide.)

Edit: in the meanwhile I found out that I should use numerical quantities close to 1 instead of SI units for micrometer length scales, because obviously Mathematica performs numerical integration with machine precision and if the numbers are too close to machine precision the numerical "signal-to-noise-ratio" can cause singularities in the solution which causes in turn NDSolve to get stuck or to blow up the required memory.

So now I am using the following values (instead of 10^-6):

xBound = 12; yBound = 12; zBound = 5;
nVak = 1.; nMat = 1.5; λ = 0.8; k = (2 π)/λ;
n0 = 1.48 (*this is the "arbitrarily chosen" because yet unknown propagation constant*);
waveGuideR = 1; (*wave guide radius*)
n[x_, y_] := If[x^2 + y^2 <= waveGuideR^2, nMat, nVak];
(*let's have a circular waveguide with refractive index of 1.5*)
theoReflCoeff = 10.^-2;(*1/theoReflCoeff is the theoretical damping
coefficient of the PML layer*)
pmlWidth = 1;(*size of the PML-layer*)

αx[x_] := Piecewise[{{1 - I 3 λ (x + (xBound - pmlWidth))^2/(4 π nVak pmlWidth^3) 
    Log[1/theoReflCoeff], x < -(xBound - pmlWidth)}, {nVak, -(xBound - pmlWidth)
    <= x <= (xBound - pmlWidth)}, {1 - I 3 λ (x - (xBound - pmlWidth))^2/(4 π nVak 
    pmlWidth^3) Log[1/theoReflCoeff], x > (xBound - pmlWidth)}}];
αy[y_] := Piecewise[{{1 - I 3 λ (y + (yBound - pmlWidth))^2/(4 π nVak pmlWidth^3) 
    Log[1/theoReflCoeff], y < -(yBound - pmlWidth)}, {nVak, -(yBound - pmlWidth) 
    <= y <= (yBound - pmlWidth)}, {1 - I 3 λ (y - (yBound - pmlWidth))^2/(4 π nVak 
    pmlWidth^3) Log[1/theoReflCoeff], y > (yBound - pmlWidth)}}];

The PML - coordinate transformation is formed such that outside the PML layer the derivatives are multiplied by one and inside the PML are multiplied by a complex value (so that they should cause damping on the wave). For clarity the imaginary part and the absolute part of the PML layer is shown in the following:

Of course the Helmholtz equation in input line 6 is actually a vector in 3 dimensions. However, I am interested only in the evolution of the field in $x$ and $y$ direction. That is why I will neglect the calculation for the field in $z$-dimension. (Here I am not totally sure that what I am doing is physically correct, since I chose for $\text{$\psi $z}$ to be zero in the initial calculation, but due to the Curl-operation I still get the contribution in the z-dimension. Anyway, if plug eqs[[3]] into NDSolve Mathematica tells me that the system is overdetermined).

The 2 coupled differential equations that have to be solved are then:

diffEq1 = eqs[[1]] == 0;
diffEq2 = eqs[[2]] == 0;

Now I define at the edges (after the PML layers) periodic boundary conditions for $x$ and $y$ directions and fields.

boundDirectionXfieldX = ψx[-xBound, y, z] == ψx[xBound, y, z];
boundDirectionYfieldX = ψx[x, -yBound, z] == ψx[x, yBound, z];
boundDirectionXfieldY = ψy[-xBound, y, z] == ψy[xBound, y, z];
boundDirectionYfieldY = ψy[x, -yBound, z] == ψy[x, yBound, z];

Furthermore I need to define the starting field that is launched into the waveguide. In this case I assume the y-component of the field to be zero, because due to the coupling of the equations also higher and more complicated transversal field modes should develop as the they indeed do in real experiments. Since I am looking only for solutions that will be guided it is unimportant what kind of initial field I chose. (Even if it is the wrong field distribution for the correct $n0$ a stable distribution in the waveguide should evolve.)

startCondDirectionZfieldX = ψx[x, y, 0] == E^(-((x^2 + y^2)/waveGuideR^2));
startCondDirectionZfieldY = ψy[x, y, 0] == 0;

Because there is a second derivative of $\psi$ in $z$ - direction, I need additionally a Neumann - boundary condition for the field launched into the waveguide at $z=0$. I think this is physically correct, because if the waveguide is infinite in z-direction and I manage somehow to generate a field with the above shape along an extended length of the waveguide then all the physical effects should still occur after this. (Correct me if I'm wrong here)

neumannCondFieldX = Derivative[0, 0, 1][ψx][x, y, 0] == 0;
neumannCondFieldY = Derivative[0, 0, 1][ψy][x, y, 0] == 0;

With this I can set up a equation system for NDSolve to solve. Additionally I used EvaluationMonitor to have at least an inkling where NDSolve currently is as well as a MemoryConstrained evaluation limited to 13 GB of RAM (but I don't it works the way I implemented it)

MemoryConstrained[
 Monitor[Fkt = {ψx, ψy} /. First@NDSolve[{diffEq1, diffEq2,
   boundDirectionXfieldX, boundDirectionYfieldX,
   boundDirectionXfieldY, boundDirectionYfieldY,
   startCondDirectionZfieldX, startCondDirectionZfieldY,
   neumannCondFieldX, neumannCondFieldY}, {ψx, ψy}, {z, 
   0., zBound}, {x, -xBound, xBound}, {y, -yBound, yBound}, 
  Method -> {"MethodOfLines","SpatialDiscretization" -> {"TensorProductGrid",
  "DifferenceOrder" -> "Pseudospectral"}}, 
  EvaluationMonitor :> (stepz = z)], AngularGauge[Dynamic[stepz/zBound], {0, 1}, 
  GaugeLabels -> Automatic]], 13 2^30]

The most puzzling part during is the error message:

NDSolve::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable y. 

This is where the strange things happen and I would like to know how to arrive at a proper well-behaved solution:

• Mathematica gives me an answer in spite of the above error message. The solution is not well behaved because it changes, depending on how large I chose the size of the boundaries. Even more troubling is the fact that for some combinations of the boundary sizes Mathematica Kernel hangs itself up as it uses almost the full available memory.
• The obvious thing to do would be to change the number of grid points. However if I chose smaller or larger number of "MaxPoints" than the default 100, NDSolve runs always into a memory limit. I wonder how this can be if I chose a smaller size? I would have thought that in such case the solution would be just less precise?
• The setting "DifferenceOrder"->"Pseudospectral" seems to be paramount. Anything else runs into the above memory problems. Only thanks to the postings of Complex valued 2+1D PDE Schrödinger equation, numerical method for `NDSolve`? I was able to get any result at all.

Here is the result calculated and displayed:

zSlices = Table[Plot3D[Abs[Fkt[[1]][x, y, z]], {x, -xBound, xBound}, 
    {y, -yBound, yBound}, PlotRange -> {0, All}, PlotPoints -> 200], 
    {z, 0, zBound, zBound/6.}];
Export["UnstableSolution.gif", zSlices, AnimationRepetitions -> Infinity, 
    "DisplayDurations" -> .4]

Edit: Another culprit in my calculations seems to be the PML itself. If I discard the PML and use the normal Helmholtz-equation without any coordinate transformations and just pure periodic boundary condition I get a physically plausible solution:

But my joy is slightly marred by the fact that since I want to find stable propagating solutions in the waveguide, I must somehow get rid off the "reflected noisy waves". I would greatly appreciate if somebody could help me with the perfectly matched layer conditions. Thank you!

Edit: you can probably skip the next part, because as long as my university doesn't give me access to Mathematica 10.4 I will be stuck on the Finite Element Method. But it would be still cool to know if the Finite Element method makes better work of the PML-conditions than the pseudospectral decomposition I'm forced to use above.. ;-)

Ok, the next thing I thought to use were the new capabilities of Mathematica considering Finite Elements for calculation, because probably the observed issues are caused by the discontinuities on the waveguide or the PML transformation of cartesian coordinates. In such a case the Method in NDSolve must be changed as in the following, if I understand Mathematica help correctly. Although I had to change the boundary conditions to be zero instead of periodic and starting condition to be 0 at the boundaries to be consistent everywhere...

Fkt = {ψx, ψy} /. First@NDSolve[{diffEq1, diffEq2,
 ψx[-xBound, y, z] == ψx[xBound, y, z] == 0,
 ψx[x, -yBound, z] == ψx[x, yBound, z] == 0,
 ψy[-xBound, y, z] == ψy[xBound, y, z] == 0,
 ψy[x, -yBound, z] == ψy[x, yBound, z] == 0,
 ψx[x, y, 0] == If[x^2 + y^2 <= waveGuideR^2, 1, 0],
 ψy[x, y, 0] == 0, 
 Derivative[0, 0, 1][ψx][x, y, 0] == 0, 
 Derivative[0, 0, 1][ψy][x, y, 0] == 0},
 {ψx, ψy}, {z, 0., 20 zBound}, {x, -xBound, xBound}, {y, -yBound, yBound},     
 Method -> {"PDEDiscretization" -> {"MethodOfLines", 
    "SpatialDiscretization" -> "FiniteElement"}}]

Now I get two new errors:

 CompiledFunction::cfex: Could not complete external evaluation at instruction 10; proceeding with uncompiled evaluation.
 NDSolve::femdpop: The FEMStiffnessElements operator failed. 

where one of them has already been asked about in Error message for FEMStiffnessElements. @ilian mentioned that this problem has been resolved in Mathematica 10.4.

The problem is that the newest Mathematica version my university is able to provide is 10.3 (I don't know how long it will take for the data processing center to distribute the 10.4 version).

So I currently I have no opportunity to see whether the 10.4 can solve my problem so that I just wait a little bit or if I have to find a different way, e.g. to use COMSOL, or something similar, with which I'm not really familiar...

Thanks for bearing with me and my rantings up to the end!

Answer 1

I'm not that familiar with electromagnetism, either, but I think there're at least 4 issues in your solving process:

1. There's no need to "consider only the magnetic field", because electric field is also continuous at the boundary of the wave guide in your case.

2. Your definition for curl in stretched coordinate is wrong. When the scale factor is $(s_x(x), s_y(y), s_z(z))$, the generalized curl should be $(\frac{\psi _z{}^{(0,1,0)}(x,y,z)}{s_{\text{y}}(y)}-\frac{\psi _y{}^{(0,0,1)}(x,y,z)}{s_{\text{z}}(z)},\frac{\psi _x{}^{(0,0,1)}(x,y,z)}{s_{\text{z}}(z)}-\frac{\psi _z{}^{(1,0,0)}(x,y,z)}{s_{\text{x}}(x)},\frac{\psi _y{}^{(1,0,0)}(x,y,z)}{s_{\text{x}}(x)}-\frac{\psi _x{}^{(0,1,0)}(x,y,z)}{s_{\text{y}}(y)})$.

3. It's not clear to me if you're dealing with a 2D problem or 3D problem, but whatever it is, the deduction for Helmholtz equation seems to be wrong. For 2D case, there should be no derivative of z; for 3D case, "fast oscillations in z direction" should be introduced as a nonhomogeneous term of the equation.

4. In the final step, you're trying to solve initial-boundary value problem of Helmholtz equation, but it's a ill-posed problem. There exists techniques for dealing with this problem of course, but the standard approach is to set up a boundary value problem.

…Well, I'm not sure if the material you're referring to is improper, but, to have a better understanding for PML, you can refer to e.g. this, this or this. To help you understand PML better, in the rest part of this answer, I'll show you my implementation for SC-PML in 2D case. ($\text{TE}^\text{z}$ mode. )

First, we need the equation. The specific formulas can be found in numerous materials, but here I'll deduce the governing equation

$$\nabla _s\times {\mu^{-1} }{\nabla _s\times E}- \omega ^2 \epsilon E =-i \omega J$$

with Mathematica to make the code more instructive and elegant. The key point is implementing the generalized curl $\nabla _s\times$. There're many possible solutions, for example:

Cross[{d@x/sx, d@y/sy, d@z/sz}, {f, g, h}[x, y, z] // Through] /. d[v_] h_ :> D[h, v]

Or make use of DChange:

DChange[Curl[{f[x, y, z], g[x, y, z], h[x, y, z]}, {x, y, z}], {x == xx s["x"], 
  y == yy s@"y", z == zz s@"z"}, {x, y, z}, {xx, yy, zz}, {f[x, y, z], g[x, y, z], 
  h[x, y, z]}]

But I think the simplest approach is the one mentioned in tutorial/VectorAnalysis:

inde = {x, y, z};    
sf = s[ToString@#]@# & /@ inde    
vf = Times @@ sf    
curlS = (sf Curl[sf #, {x, y, z}]/vf) &

The remaining part is straightforward:

Εlst = Ε[ToString@#] @@ Most@inde & /@ inde /. Ε["z"] -> (0 &)    
jlst = j[ToString@#] @@ Most@inde & /@ inde /. j["z"] -> (0 &)    
eqnS = Simplify[
  curlS[curlS@Εlst/mu0] - omega^2 e0 Εlst == -I omega jlst // Thread // Most, 
    {mu0 > 0, s[ToString@#]@# != 0 & /@ inde} // Flatten]

Next, define the formula for scale factor $s$. Here I simply follow the one in the paper linked above, as far as I can tell, this is also the most popular way to define $s$:

sigmamax[thick_] = -(((m + 1) Log@R)/(2 eta0 thick))    
sigma[l_, thick_] = sigmamax[thick] (l/thick)^m    
sgenerator[x_, {lb_, rb_}, thick_] = 
 Piecewise[{{1 - I sigma[x - (rb - thick), thick]/(omega e0) , 
    x > rb - thick}, {1 - I sigma[(lb + thick) - x, thick]/(omega e0) , x < lb + thick}},
   1]

{sfunc["x"][x_], sfunc["y"][y_]} = 
 MapThread[sgenerator, {{x, y}, {{lb@#, rb@#}, {lb@#2, rb@#2}}, {th@#, th@#2}} &["x", 
    "y"]] // Simplify

Then substitute specific values into the equation:

lam = 532 10^-9;

domain = {-2 lam, 2 lam};
thickness = 4 lam/10;
pdeS = Block[{e0 = 8854/10^3*10^-12, mu0 = 1257/10^3*10^-6, conduct = 10^7}, 
   Block[{omega = (2 Pi)/(lam Sqrt[mu0 e0]), m = 4, R = E^-16, eta0 = Sqrt[mu0/e0], j, 
     lb, rb, th}, ({{lb@#, rb@#}, {lb@#2, rb@#2}, {th@#, th@#2}, {j[#][x, y], 
          j[#2][x, y]}} = {domain, 
         domain, {thickness, thickness}, {conduct Exp[10 (-(x/lam)^2 - (y/lam)^2)], 
          0}}) &["x", "y"]; eqnS /. s -> sfunc]];
bc = Function[{x, 
      y}, {Ε["x"][x, y] == 0, Ε["y"][x, y] == 
       0}] @@@ {{domain[[1]], y}, {domain[[2]], y}, {x, domain[[1]]}, {x, domain[[2]]}} //
    Flatten;

Last step is to solve the equation set. We can solve {pdeS, bc} with NDSolve directly if $VersionNumber >= 11.1:

<< NDSolve`FEM`  
mesh = ToElementMesh[FullRegion[2], {domain, domain}, "MaxCellMeasure" -> lam/10^9]; 
bcfem = DirichletCondition[{Ε["x"][x, y] == 0, Ε["y"][x, y] == 0}, True]; 
sol = NDSolveValue[{pdeS, bcfem}, {Ε["x"], Ε["y"]}, 
  Element[{x, y}, mesh]]
GraphicsGrid[
   With[{domain = Sequence @@ (domain + {thickness, -thickness})}, 
     Outer[Plot3D[#1[First[#2][x, y]], Evaluate[{x, domain}], 
           Evaluate[{y, domain}], PlotRange -> All, 
           PlotLabel -> #1[Last[#2]]] & , {Re, Im}, 
       {{sol[[1]], "Ex"}, {sol[[2]], "Ey"}}, 1]], 
   ImageSize -> Large]

If you're still in or before v9 (where "FiniteElement" isn't introduced yet), or between v10.0 and v11.0 (where a bug isn't fixed yet), finite difference method (FDM) can be used for solving the problem. I'll use pdetoae for discretizing:

points = 50;
grid = Array[# &, points, domain];
difforder = 2;
(*Definition of pdetoae isn't included in this code piece,
  please find it in the link above.*)
ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder];
del = Most@Rest@# &;

ae = del /@ del@# & /@ ptoa@Simplify`PWToUnitStep@pdeS;
aebc = MapAt[del, ptoa@bc, List /@ Range@4];
{b, mat} = CoefficientArrays[{ae, aebc} // Flatten, 
   Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten];

sollst = LinearSolve[-mat, N@b];
(* Alternatively, if you're confused about del: *)
(*
fullae = ptoa@Simplify`PWToUnitStep@pdeS;
fullaebc = ptoa@bc;
{b, mat} = CoefficientArrays[{fullae, fullaebc} // Flatten, 
   Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten];
sollst = LeastSquares[-mat, N@b, Method -> Direct]; // AbsoluteTiming
 *)
solmat = ArrayReshape[sollst, {2, points, points}];
{solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat;

GraphicsGrid[
 With[{domain = Sequence @@ (domain + {thickness, -thickness})}, 
  Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], 
     PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", 
     PlotPoints -> 50] &, {Re, 
    Im}, {{solΕx, 
     "\!\(\*SubscriptBox[\(Ε\), \(x\)]\)"}, {solΕy, 
     "\!\(\*SubscriptBox[\(Ε\), \(y\)]\)"}}, 1]], ImageSize -> Large]