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!
nMat
and nVak` for example. Write a function that gets all parameters as arguments and that returns your equation with them filled in. Try a simpler example first and add more stuff later. The FEM message means that no all PDE coefficients could be compiled and for the complex value case you'd need 10.4. With the tensor product gird you might get away with an earlier version. $\endgroup$NDSolve
. Do you mean that the problem comes from theSetDelayed
(:=
) definition of my refractive index function? $\endgroup$NDSolve::mxsst:
that you report here. For the FEM code, is this a stationary problem - then why MethodOfLines? If it's transient which is the time variable? For the FEM code, those warnings will disappear if you remove the constants from your definitions, as those, when compiled will need to call back. So that for example?n
just depends onx
andy
. $\endgroup$