Taylor done

Author: Jon Schlipf

README.md

@@ -143,6 +143,13 @@ The CUDA implementation is heavily inspired by [this fork](https://github.com/al
 A major exception is the solution of the non-hermitian Eigenvalue problems, which is not implemented in CUDA. Here, a fallback CPU routine (`LinearAlgebra.eigen`) is used. 
 For a truncation order of `N=10`, GPU acceleration achieves a speedup of a factor of 5 in a system with dual AMD EPYC 7713 64-Core Processor and NVIDIA H100 accelerator.
 
+## Eigenvalue-free algorithm
+
+Xu and Charlton (references below) have introduced an alternative way, computing the matrix exponential by a Taylor series instead of an Eigentransform. This allows for efficient parallelism and thus improves performance drastically, about two orders of magnitude compared to ETM using high-end GPU. Only reflection and transmission are implemented so far, accessed via:
+
+```julia
+E,H=taylor_reftra(mdl,grd,xypoints,zpoints,λ,ste) #compute the electric and magnetic field
+```
 
 ## Mathematics
 
@@ -163,3 +170,5 @@ Marco Liscidini, Dario Gerace, Lucio Claudio Andreani, and J. E. Sipe, "Scatteri
 M. G. Moharam and T. K. Gaylord, "Rigorous coupled-wave analysis of planar-grating diffraction," J. Opt. Soc. Am. 71, 811-818 (1981) 
 
 Raymond Rumpf, "Improved formulation of scattering matrices for semi-analytical methods thatis consistent with convention," Progress In Electromagnetics Research B35, 241–261 (2011)
+
+J. Xu and M. D. B. Charlton, "Highly Efficient Rigorous Coupled-Wave Analysis Implementation Without Eigensystem Calculation," IEEE Access 12, 127966-127975 (2024)

src/Taylor/Taylor.jl

@@ -1,9 +1,18 @@
-using LinearAlgebra
 using LinearAlgebra,CUDA
 using ..Common
 
-export squarescale_reftra
-taylorx(dnx,dny,Kx,Ky,λ,l::SimpleLayer)=eigenmodes(dnx,dny,Kx,Ky,λ,l).X
+function taylorx(dnx,dny,Kx,Ky,λ,l::SimpleLayer)
+    #Xa=(eigenmodes(dnx,dny,Kx,Ky,λ,l).X
+    #return [Xa 0Xa;0Xa Xa]
+    IMa=Diagonal(Kx*0 .+1)
+    k0=2π/λ
+    za=0*IMa
+    ε=get_permittivity(l.material,λ)
+    m1=sqrt.(Complex.(Kx.*Kx+Ky.*Ky-ε*I))
+    m2=exp.(m1*k0*l.thickness)
+    lm=[m2 za za za; za m2 za za;za za m2 za;za za za m2]
+    return lm
+end
 function taylorx(dnx,dny,Kx,Ky,λ,l::AnisotropicLayer)
     # probably doable without taylor as well?
     IMa=Diagonal(Kx*0 .+1)
@@ -156,49 +165,6 @@ function taylorx(dnx,dny,Kx,Ky,λ,l::PatternedLayer)
     X=B2+(B3+A9)*A9
     return X^(2^m)
 end
-function squarescalex(dnx,dny,Kx,Ky,λ,l::PatternedLayer)
-    IMa=Diagonal(Kx*0 .+1)
-    IM=[IMa 0*IMa;0*IMa IMa]
-    k0=2π/real(λ)
-    #get the base permittivity
-    εxx=get_permittivity(l.materials[1],λ,1)*I
-    #add the permittivity for all inclusions
-    if minimum([typeof(m)<:Common.Isotropic for m in l.materials])
-        #all isotropic
-        εxx=get_permittivity(l.materials[1],λ)*I
-        for ct=1:length(l.geometries)
-            rec=reciprocal(l.geometries[ct],dnx,dny)
-            εxx+=rec*(get_permittivity(l.materials[ct+1],λ)-get_permittivity(l.materials[ct],λ))
-        end
-        εzz=εyy=εxx
-        εxy=εyx=0I
-    else
-        #anisotropic
-        εxx=get_permittivity(l.materials[1],λ,1)*I
-        εxy=get_permittivity(l.materials[1],λ,2)*I
-        εyx=get_permittivity(l.materials[1],λ,3)*I
-        εyy=get_permittivity(l.materials[1],λ,4)*I
-        εzz=get_permittivity(l.materials[1],λ,5)*I
-        for ct=1:length(l.geometries)
-            rec=reciprocal(l.geometries[ct],dnx,dny)
-            εxx+=rec*(get_permittivity(l.materials[ct+1],λ,1)-get_permittivity(l.materials[ct],λ,1))
-            εxy+=rec*(get_permittivity(l.materials[ct+1],λ,2)-get_permittivity(l.materials[ct],λ,2))
-            εyx+=rec*(get_permittivity(l.materials[ct+1],λ,3)-get_permittivity(l.materials[ct],λ,3))
-            εyy+=rec*(get_permittivity(l.materials[ct+1],λ,4)-get_permittivity(l.materials[ct],λ,4))
-            εzz+=rec*(get_permittivity(l.materials[ct+1],λ,5)-get_permittivity(l.materials[ct],λ,5))
-        end
-    end	 	
-    η=inv(εzz)
-    P=[Kx*η*Ky I-Kx*η*Kx;Ky*η*Ky-I -Ky*η*Kx]
-    Q=[Kx*Ky+εyx εyy-Kx*Kx;Ky*Ky-εxx -εxy-Ky*Kx]
-    A0=[0IM P;Q 0IM]*k0*l.thickness
-    nrm=maximum(sum(abs.(A0),dims=1))
-    m=Int(ceil(log2(nrm)))
-    m=0
-    A=A0*2.0^-m
-    X=exp(A)
-    return X^(2^m)
-end
 function taylor_reftra(ψin,m::RCWAModel,grd::RCWAGrid,λ)
     IMa=Diagonal(grd.Kx*0 .+1)
     IM=[IMa 0*IMa;0*IMa IMa]
@@ -221,25 +187,3 @@ function taylor_reftra(ψin,m::RCWAModel,grd::RCWAGrid,λ)
     return R,T
 
 end
-function squarescale_reftra(ψin,m::RCWAModel,grd::RCWAGrid,λ)
-    IMa=Diagonal(grd.Kx*0 .+1)
-    IM=[IMa 0*IMa;0*IMa IMa]
-    IMb=[IM 0*IM;0*IM IM]
-    X=[squarescalex(grd.dnx,grd.dny,grd.Kx,grd.Ky,λ,l) for l in m.layers]
-    Xp=IMb
-    for X in X
-        Xp*=X
-    end
-    ref=halfspace(grd.Kx,grd.Ky,m.εsup,λ) #superstrate and substrate
-    tra=halfspace(grd.Kx,grd.Ky,m.εsub,λ)
-    Y=Xp*[IM;-tra.V]
-    S=[Y [-IM;-ref.V]]\[IM;-ref.V]*ψin
-    to,ro=slicehalf(S)
-
-    kzin=grd.k0[3]#*real(sqrt(get_permittivity(m.εsup,λ)))
-    R=a2p(0ro,ro,ref.V,IM,kzin) #compute amplitudes to powers
-    T=-a2p(to,0to,tra.V,IM,kzin)
-
-    return R,T
-
-end

test/runtests.jl

@@ -59,7 +59,7 @@ n1=1+10rand()
 θ=88rand()+1
 α=360rand()
 grd=rcwagrid(0,0,100rand(),100rand(),θ,α,λ,ConstantPerm(n1^2))
-ste,stm=rcwasource(grd,real(n1))
+ste,stm=rcwasource(grd,(n1))
 ref=RigorousCoupledWaveAnalysis.halfspace(grd.Kx,grd.Ky,ConstantPerm(n1^2),λ)
 P1=-RigorousCoupledWaveAnalysis.a2p(ste,0ste,ref.V,I,grd.k0[3])
 P2=-RigorousCoupledWaveAnalysis.a2p(stm,0ste,ref.V,I,grd.k0[3])
@@ -75,7 +75,7 @@ eps4=10rand()+10im*rand()
 λ=1000rand()
 mdl=RCWAModel([PatternedLayer(100rand(),[ConstantPerm(eps2),ConstantPerm(eps3)],[Circle(rand())])],ConstantPerm(eps1),ConstantPerm(eps4))
 grd=rcwagrid(1,1,100rand(),100rand(),88rand()+.1,360rand(),λ,ConstantPerm(eps1))
-ste,stm=rcwasource(grd,real(√eps1))
+ste,stm=rcwasource(grd,(√eps1))
 Rte1,Tte1=etm_reftra(ste,mdl,grd,λ) 
 Rtm1,Ttm1=etm_reftra(stm,mdl,grd,λ) 
 Rte2,Tte2=srcwa_reftra(ste,mdl,grd,λ) 
@@ -95,7 +95,7 @@ eps4=10rand()
 λ=1000rand()
 mdl=RCWAModel([PatternedLayer(100rand(),[ConstantPerm(eps2),ConstantPerm(eps3)],[Circle(rand())])],ConstantPerm(eps1),ConstantPerm(eps4))
 grd=rcwagrid(1,1,100rand(),100rand(),88rand()+.1,360rand(),λ,ConstantPerm(eps1))
-ste,stm=rcwasource(grd,real(√eps1))
+ste,stm=rcwasource(grd,(√eps1))
 Rte1,Tte1=etm_reftra(ste,mdl,grd,λ) 
 Rtm1,Ttm1=etm_reftra(stm,mdl,grd,λ) 
 Rte2,Tte2=srcwa_reftra(ste,mdl,grd,λ) 
@@ -111,4 +111,24 @@ Rtm2,Ttm2=srcwa_reftra(stm,mdl,grd,λ)
 end
 
 
+@testset "RT_TaylorVSETM_simple" begin
+eps1=10rand()+10im*rand()
+eps2=10rand()+10im*rand()
+eps3=10rand()+10im*rand()
+eps4=10rand()+10im*rand()
+λ=1000rand()
+px=100rand()
+py=100rand()
+    mdl=RCWAModel([PatternedLayer(minimum([px,py])*rand(),[ConstantPerm(eps2),ConstantPerm(eps3)],[Circle(rand())])],ConstantPerm(eps1),ConstantPerm(eps4)) # thickness > pitch for sufficiently accurate Taylor
+grd=rcwagrid(1,1,px,py,88rand()+.1,360rand(),λ,ConstantPerm(eps1))
+ste,stm=rcwasource(grd,(√eps1))
+Rte1,Tte1=etm_reftra(ste,mdl,grd,λ) 
+Rtm1,Ttm1=etm_reftra(stm,mdl,grd,λ) 
+Rte2,Tte2=taylor_reftra(ste,mdl,grd,λ) 
+Rtm2,Ttm2=taylor_reftra(stm,mdl,grd,λ) 
+@test isapprox(Rte1,Rte2,atol=1E-5)
+@test isapprox(Rtm1,Rtm2,atol=1E-5)
+@test isapprox(Tte1,Tte2,atol=1E-5)
+@test isapprox(Ttm1,Ttm2,atol=1E-5)
+end