I am using FreeFEM++ to simulate the cylindrical waveguide, but no reasonable solution is returned. Here I’m using the 94 GHz light (𝜆 = 3.17mm) as input, and radius of the waveguide is 1.27mm.

If I uncomment the `savevtk`

lines to seed the solutions, nothing appeared:

Need I provide more information to help debugging?

Here is the code:

```
// heavily copied from https://github.com/FreeFem/FreeFem-sources/blob/master/examples/plugin/waveguide.edp
load "MUMPS"
load "iovtk"
load "medit"
load "msh3"
load "Element_Mixte3d"
// The boundary value problem:
// (sigma = 0, here k is the wavenumber)
// -k^2*E + curl(curl E) = 0 in Omega
// E x n = 0 on the side
// Curl(E) x n + i*beta n x (E x n) = Ginc on z = 0
// Curl(E) x n + i*beta n x (E x n) = 0 on z = c
bool debug = true;
real r = 0.00127;
real d = r;
int nx = 20;
int nz = 20;
real zmin = 0.0;
real zmax = 0.01;
// build the mesh
border C(t=0.0, 2.0*pi) {
x=r*cos(t);
y=r*sin(t);
label=1;
};
mesh Thx = buildmesh(C(nx));
int [int] rup = [0, 1], rdown = [0, 2], rmid = [1, 3];
int guide = 3, in = 1, out = 2;
mesh3 Th = buildlayers(Thx, nz, zbound=[zmin, zmax]
, coef = 1.
, labelup = rup
, labelmid = rmid
, labeldown = rdown
);
//medit("Cyl", Th, wait=debug);
//plot(Th, wait=debug);
cout << "labels(Th) = " << labels(Th) << endl;
real f = 94*10^9;
real er = 1;
real c0 = 299792458;
real k = 2*pi*f*sqrt(er) / c0;
int m = 1, n = 1;
real beta = sqrt(k^2 - (m*pi/d)^2 - (n*pi/d)^2);
real Z = sqrt(er) * 120 * pi;
real ukb = 1 / (k^2 - beta^2);
func expbz = exp(-1i*beta*z);
func ExTE = (1i*k*Z) * ukb * (n*pi) / d * cos(m*pi*x/d) * sin(n*pi*y/d) * expbz;
func EyTE = -(1i*k*Z) * ukb * (n*pi) / d * sin(m*pi*x/d) * cos(n*pi*y/d) * expbz;
func Gix = +2*1i*beta*ExTE;
func Giy = +2*1i*beta*EyTE;
func Giz = 0;
fespace Nh(Th, Edge13d);
Nh<complex> [Ex, Ey, Ez], [vx, vy, vz];
// Macros
macro Curl(ux,uy,uz) [dy(uz)-dz(uy),dz(ux)-dx(uz),dx(uy)-dy(ux)] // EOM
macro Nvec(ux,uy,uz) [uy*N.z-uz*N.y,uz*N.x-ux*N.z,ux*N.y-uy*N.x] // EOM //uxN
macro Curlabs(ux,uy,uz) [abs(dy(uz)-dz(uy)),abs(dz(ux)-dx(uz)),abs(dx(uy)-dy(ux))] //EOM
problem waveguide([Ex, Ey, Ez], [vx, vy, vz], solver=sparsesolver) =
int3d(Th)(Curl(Ex,Ey,Ez)' * Curl(vx,vy,vz))
- int3d(Th)(k^2*[Ex,Ey,Ez]' * [vx, vy, vz])
+ int2d(Th,in,out)(1i * beta * Nvec(Ex,Ey,Ez)' * Nvec(vx,vy,vz))
- int2d(Th,in)([vx,vy,vz]' * [Gix,Giy,Giz])
+ on(guide,Ex=0,Ey=0,Ez=0);
waveguide;
cout << "r = " << r << " Integrated Ex = " << int3d(Th)(abs(Ex)^2) << endl;
medit("real", Th, [real(Ex),real(Ey),real(Ez)]);
//cout << (abs(Ex)^2) << end;
//cout << (abs(Ey)^2) << end;
//cout << (abs(Ez)^2) << end;
//int[int] Order = [1];
//savevtk("result_x.vtu", Th, abs(real(Ex)), dataname="RealEx", order=Order);
//savevtk("result_y.vtu", Th, abs(real(Ey)), dataname="RealEy", order=Order);
//savevtk("result_z.vtu", Th, abs(real(Ez)), dataname="RealEz", order=Order);
//medit("real", Th, [real(Ex),real(Ey),real(Ez)]);
```

and the log