Tutorial: Problem from nanophotonics

One of the key challanges within the field of nanophotonics stems from the need to be able to compute the modes of frequency-despersive structures in a reliable and efficient way. The frequency (i.e. eigenvalue) dependency can be viewed as a nonlinearity and therefore naturally leads to nonlinear eigenvalue problems.

This tutorial is based on a research preprint where a model is set up and solved with SLEPc. An associated repository of code is available which should be used in combination with gmsh and onelab. Credit for the discretization and application should go to the authors of the paper, in particular Guillaume Demésy for providing the model files online.

We will use the model files and onelab for that simulation and reproduce computational results using NEP-PACK.

Part 1: Setup the matrices

You need to download the code with the model of the frequency-dispersive media here. (Note that we have fixed the specific version of the the model code in the URL.) You also need to install ONELAB from http://onelab.info/.

We need the matrices generated by this model, which are not saved to disk by default. A small modification of one of the project files is needed. You need to modify the NonLinearEVP.pro on the line before the EigenSolve specification for the res_NEP_E:

Print[M1]; // Add this line
EigenSolve[M1,neig,eig_target_re,eig_target_im,EigFilter[],
     { {1}, {-eps_oo_1,gam_1*eps_oo_1, -om_d_1^2,0}, {-1,0,0} },
     { {1}, {1,-gam_1},                              {1} } ];

That is, you should add the text Print[M1]; (in the current version) before line 354.

After you do this modification, load the model and click "run" in the Gmsh tool, you will obtain files in the current directory containing the FEM-discretizations needed to set up the problem (file_mat_MX.m.bin).

Part 2: Implementation in NEP-PACK

The NEP in this problem has the structure

The constants are given in the project file and we set them in our julia code:

julia> a_lat=50;
julia> cel=a_lat/(2*pi);
julia> nrm=a_lat/(2*pi*cel);
julia> om_d_1=1.1;
julia> gam_1=0.05;
julia> om_d_1=om_d_1/nrm;
julia> gam_1=gam_1/nrm;
julia> eps_oo_1=1.0;

The NEP in this example can be conveniently expressed as a SPMF_NEP, where the first function is constant, the second term is a rational function and the third is a quadratic term. We define them in a matrix function sense:

julia> f1=s-> one(s) #
julia> f2=s-> (s-gam_1*one(s))\(-eps_oo_1*s^3+gam_1*eps_oo_1*s^2-om_d_1^2*s)
julia> f3=s-> -s^2

If you have carried out Part 1, you should have the sparse discretization matrices available. They are stored in the PETSc-binary format. NEP-PACK contains functionality to load the generated matrices using the function Gallery.naive_petsc_read.

Suppose gmsh_files is a path to the bin-files generated in Part 1. These commands load the matrices

julia> A3=Gallery.naive_petsc_read(joinpath(gmsh_files,"file_mat_M15.m.bin")); # Switched order is intentional
julia> A2=Gallery.naive_petsc_read(joinpath(gmsh_files,"file_mat_M16.m.bin"));
julia> A1=Gallery.naive_petsc_read(joinpath(gmsh_files,"file_mat_M17.m.bin"));

The SPMF is created directly

julia> nep=SPMF_NEP([A1,A2,A3], [f1,f2,f3]);

We can now apply a method of choice. In this case we solve it with the tensor infinite Arnoldi method (tiar) with a shift/target the same as in the project files

julia> (λ,V)=tiar(nep,σ=0.00243604+0.366703im,logger=1,neigs=10,maxit=100);
-
--
=--
+---
+----
+-----
+------
+-------
+--------
+---------
+----------
+==---------
+++----------
+++=----------
++++-----------
++++------------
++++-------------
++++--------------
++++---------------
++++----------------
++++==---------------
++++++----------------
++++++-----------------
++++++=-----------------
++++++=------------------
+++++++-------------------
+++++++--------------------
+++++++---------------------
+++++++----------------------
+++++++==---------------------
+++++++==----------------------
+++++++++=----------------------
+++++++++=-----------------------
++++++++++------------------------
julia> λ
10-element Array{Complex{Float64},1}:
 0.002436044429913607 + 0.3667026531004412im
 0.001175454247612957 + 0.3748476897696621im
 0.006531655269175296 + 0.3736695833524133im
 0.013279531609677153 + 0.37899563698023087im
  0.04259449677920191 + 0.38730647107316im
  0.04388349538759248 + 0.39266023012447837im
 0.007774845736605659 + 0.25923501436468327im
  0.09285578050590135 + 0.41819550232817293im
  0.09374197960056296 + 0.41522753073470614im
  0.01578630143214552 + 0.49172628247971106im

Each row in the logger printout of tiar corresponds to an iteration. The sign - corresponds to an unconverged eigenvalue, + corresponds to a converged eigenvalue and = corresponds to an eigenvalue which is almost converged (interpreted as a factor 10 from the convergence criteria). The eigenvalue 0.007774845736605659 + 0.25923501436468327im which can be found in the printout above, is reported by the default setting in the gmsh tool applied to the model files.