Last active
February 15, 2023 22:50
-
-
Save simonp0420/5d1896af5a7db8262bd3496c5111a9f9 to your computer and use it in GitHub Desktop.
Demonstrate bug when using MKL.jl and MKLSparse.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
calcpml3d(NGRID, NPML) -> (sx, sy, sz) | |
Calculate PML Parameters | |
# Input Arguments | |
- `NGRID`: Iterable containing `(Nx, Ny, Nz)`, number of cells on grid | |
- `NPML`: Iterable containing (NXLO, NXHI, NYLO, NYHI, NZLO, NZHI), the size of PML | |
# Return Values | |
- `sx`, `sy`, `sz`: Complex matrices of size `(Nx, Ny, Nz)`, resp., containing the | |
complex coordinate stretching parameters for each dimension of the grid. | |
""" | |
function calcpml3d(NGRID,NPML) | |
######################################################### | |
## HANDLE INPUT ARGUMENTS | |
######################################################### | |
# DEFINE PML PARAMETERS | |
amax = 4 | |
cmax = 1 | |
p = 3 | |
# EXTRACT GRID SIZE | |
Nx, Ny, Nz = (n for n in NGRID) | |
# EXTRACT PML SIZE | |
NXLO, NXHI, NYLO, NYHI, NZLO, NZHI = (n for n in NPML) | |
######################################################### | |
## CALCULATE PML PARAMETERS | |
######################################################### | |
# INITIALIZE PML PARAMETERS TO PROBLEM SPACE | |
sx, sy, sz = (ones(ComplexF64, Nx, Ny, Nz) for _ in 1:3) | |
# CALCULATE XLO PML | |
for nx in 1:NXLO | |
ax = 1 + (amax - 1) * (nx/NXLO)^p | |
cx = cmax * sinpi(nx/(2*NXLO))^2; | |
sx[NXLO - nx + 1, :, :] .= ax * (1 - im * 60cx) | |
end | |
# CALCULATE XHI PML | |
for nx in 1:NXHI | |
ax = 1 + (amax - 1) * (nx/NXHI)^p | |
cx = cmax * sinpi(nx/(2*NXHI))^2 | |
sx[Nx - NXHI + nx, :, :] .= ax * (1 - im * 60cx) | |
end | |
# CALCULATE YLO PML | |
for ny in 1:NYLO | |
ay = 1 + (amax - 1) * (ny/NYLO)^p | |
cy = cmax * sinpi(ny/(2*NYLO))^2 | |
sy[:, NYLO - ny + 1, :] .= ay * (1 - im * 60cy) | |
end | |
# CALCULATE YHI PML | |
for ny in 1:NYHI | |
ay = 1 + (amax - 1) * (ny/NYHI)^p | |
cy = cmax * sinpi(ny/(2*NYHI))^2 | |
sy[:, Ny - NYHI + ny, :] .= ay * (1 - im * 60cy) | |
end | |
# CALCULATE ZLO PML | |
for nz in 1:NZLO | |
az = 1 + (amax - 1) * (nz/NZLO)^p | |
cz = cmax * sinpi(nz/(2*NZLO))^2 | |
sz[:, :, NZLO - nz + 1] .= az * (1 - im * 60cz) | |
end | |
# CALCULATE ZHI PML | |
for nz in 1:NZHI | |
az = 1 + (amax - 1) * (nz/NZHI)^p | |
cz = cmax * sinpi(nz/(2*NZHI))^2 | |
sz[:, :, Nz - NZHI + nz] .= az * (1 - im * 60cz) | |
end | |
return (sx, sy, sz) | |
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
fdfd3d(dev::NamedTuple, src::NamedTuple) -> dat::NamedTuple | |
Three-dimensional FDFD for periodic structures | |
# INPUT ARGUMENTS | |
- `dev`: Device parameters, with the following fields: | |
- `ER2` or | |
`ER2xx` `ER2xy` `ER2xz` | |
`ER2yx` `ER2yy` `ER2yz` relative permittivity on 2x grid | |
`ER2zx` `ER2zy` `ER2zz` | |
- `UR2`: or | |
`UR2xx` `UR2xy` `UR2xz` | |
`UR2yx` `UR2yy` `UR2yz` relative permeability on 2x grid | |
`UR2zx` `UR2zy` `UR2zz` | |
- `RES`: `[dx, dy]` grid resolution on Yee grid | |
- `NPML`: `[NZLO, NZHI]` size of SCPML on Yee grid | |
- `src`: Source parameters with the following fields: | |
- `λ₀`: free-space wavelength | |
- `θ`: polar angle of incidence | |
- `ϕ`: azimuthal angle of incidence | |
- `pte`: complex amplitude of the TE polarization | |
- `ptm`: complex amplitude of the TM polarization | |
# RETURN VALUES | |
- `dat`: Output data with the following fields: | |
- `fx`, `fy`, `fz`: simulated field | |
- `m`, `n`: diffraction order indices | |
- `RDE`: diffraction efficiencies of reflected waves | |
- `REF`: overall reflectance | |
- `TDE`: diffraction efficiencies of transmitted waves | |
- `TRN`: overall transmittance | |
- `s11`: complex ref. coefficients of diffraction orders | |
- `s21`: complex tran. coefficients of diffraction orders | |
- `CON`: Power conservation parameter, equal to `REF + TRN` | |
""" | |
function fdfd3d(dev,src) | |
# DEFINE SOLVER PARAMETERS | |
tol = 1e-6 | |
maxit = 15000 | |
# UTILITY FUNCTIONS | |
diagonalize(x) = Diagonal(vec(x)) | |
######################################################### | |
## HANDLE INPUT ARGUMENTS | |
######################################################### | |
# GET SIZE OF GRID | |
if haskey(dev, :ER2) | |
Nx2, Ny2, Nz2 = size(dev.ER2) | |
else | |
Nx2, Ny2, Nz2 = size(dev.ER2xx) | |
end | |
# CALCULATE GRID | |
Nx, Ny, Nz = (Nx2, Ny2, Nz2) .÷ 2 | |
dx, dy, dz = dev.RES | |
Sx, Sy, Sz = (Nx, Ny, Nz) .* (dx, dy, dz) | |
# GRID AXES | |
xa2 = [1:Nx2;] * dx/2; xa2 .-= mean(xa2) | |
ya2 = [1:Ny2;] * dy/2; ya2 .-= mean(ya2) | |
za2 = [0.5:Nz2-0.5;] * dz/2 | |
#[Y2,X2,Z2] = meshgrid(ya2,xa2,za2) | |
# WAVE NUMBER | |
k0 = 2π / src.λ₀ | |
# SPECIAL MATRICES | |
M = Nx * Ny * Nz | |
ZZ = spzeros(M, M) | |
#I0 = speye(M,M); | |
# PERMITTIVITY TENSOR ELEMENTS | |
if haskey(dev, :ER2) | |
ERxx = diagonalize(dev.ER2[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
ERxy = ZZ | |
ERxz = ZZ | |
ERyx = ZZ | |
ERyy = diagonalize(dev.ER2[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
ERyz = ZZ | |
ERzx = ZZ | |
ERzy = ZZ | |
ERzz = diagonalize(dev.ER2[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
else | |
ERxx = diagonalize(dev.ER2xx[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
ERxy = diagonalize(dev.ER2xy[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
ERxz = diagonalize(dev.ER2xz[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
ERyx = diagonalize(dev.ER2yx[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
ERyy = diagonalize(dev.ER2yy[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
ERyz = diagonalize(dev.ER2yz[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
ERzx = diagonalize(dev.ER2zx[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
ERzy = diagonalize(dev.ER2zy[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
ERzz = diagonalize(dev.ER2zz[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
end | |
# PERMEABILITY TENSOR ELEMENTS | |
if haskey(dev, :UR2) | |
URxx = diagonalize(dev.UR2[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
URxy = ZZ | |
URxz = ZZ | |
URyx = ZZ | |
URyy = diagonalize(dev.UR2[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
URyz = ZZ | |
URzx = ZZ | |
URzy = ZZ | |
URzz = diagonalize(dev.UR2[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
else | |
URxx = diagonalize(dev.UR2xx[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
URxy = diagonalize(dev.UR2xy[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
URxz = diagonalize(dev.UR2xz[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
URyx = diagonalize(dev.UR2yx[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
URyy = diagonalize(dev.UR2yy[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
URyz = diagonalize(dev.UR2yz[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
URzx = diagonalize(dev.UR2zx[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
URzy = diagonalize(dev.UR2zy[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
URzz = diagonalize(dev.UR2zz[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
end | |
######################################################### | |
## PERFORM FDFD ANALYS | |
######################################################### | |
# EXTRACT MATERIAL PROPERTIES IN EXTERNAL REGIONS | |
ur_ref = URxx[1, 1] | |
ur_trn = URxx[M, M] | |
er_ref = ERxx[1, 1] | |
er_trn = ERxx[M, M] | |
n_ref = sqrt(ur_ref * er_ref) | |
n_trn = sqrt(ur_trn * er_trn) | |
# CALCULATE PML PARAMETERS | |
sx2, sy2, sz2 = calcpml3d([Nx2 Ny2 Nz2], 2 .* [0, 0, 0, 0, dev.NPML...]) | |
sx_ey = diagonalize(1 ./ sx2[1:2:Nx2,2:2:Ny2,1:2:Nz2]) | |
sx_ez = diagonalize(1 ./ sx2[1:2:Nx2,1:2:Ny2,2:2:Nz2]) | |
sy_ex = diagonalize(1 ./ sy2[2:2:Nx2,1:2:Ny2,1:2:Nz2]) | |
sy_ez = diagonalize(1 ./ sy2[1:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
sz_ex = diagonalize(1 ./ sz2[2:2:Nx2, 1:2:Ny2,1:2:Nz2]) | |
sz_ey = diagonalize(1 ./ sz2[1:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
sx_hy = diagonalize(1 ./ sx2[2:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
sx_hz = diagonalize(1 ./ sx2[2:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
sy_hx = diagonalize(1 ./ sy2[1:2:Nx2, 2:2:Ny2,2:2:Nz2]) | |
sy_hz = diagonalize(1 ./ sy2[2:2:Nx2, 2:2:Ny2,1:2:Nz2]) | |
sz_hx = diagonalize(1 ./ sz2[1:2:Nx2, 2:2:Ny2,2:2:Nz2]) | |
sz_hy = diagonalize(1 ./ sz2[2:2:Nx2, 1:2:Ny2,2:2:Nz2]) | |
# CALCULATE WAVE VECTORS | |
kinc = k0 * n_ref .* [sin(src.θ) * cos(src.ϕ), | |
sin(src.θ) * sin(src.ϕ), | |
cos(src.θ)] | |
m = [-(Nx ÷ 2):((Nx-1) ÷2);] | |
n = [-(Ny ÷ 2):((Ny-1) ÷ 2);] | |
kx = kinc[1] .- m * 2π / Sx | |
ky = kinc[2] .- n * 2π / Sy | |
kz_ref = sqrt.(complex(((k0*n_ref)^2 .- kx.^2 .- ky.^2))) | |
kz_trn = sqrt.(complex(((k0*n_trn)^2 .- kx.^2 .- ky.^2))) | |
# BUILD DERIVATIVE MATRICES | |
NS = [Nx, Ny, Nz] | |
RES = [dx, dy, dz] | |
BC = [1, 1, 0] | |
DEX, DEY, DEZ, DHX, DHY, DHZ = yeeder3d(NS, k0*RES, BC, kinc/k0) | |
# CALCULATE INTERPOLATION MATRICES (assumes src.θ == 0) | |
iszero(src.θ) || error("Nonzero θ=$θ not supported for RX, RY, and RZ calculation") | |
RX = (0.5 * k0 * dx) .* abs.(DEX) | |
RY = (0.5 * k0 * dy) .* abs.(DEY) | |
RZ = (0.5 * k0 * dz) .* abs.(DEZ) | |
# FORM MATERIALS TENSORS | |
ER = [ERxx RX*RY'*ERxy RX*RZ'*ERxz | |
RY*RX'*ERyx ERyy RY*RZ'*ERyz | |
RZ*RX'*ERzx RZ*RY'*ERzy ERzz] | |
UR = [URxx RX*RY'*URxy RX*RZ'*URxz | |
RY*RX'*URyx URyy RY*RZ'*URyz | |
RZ*RX'*URzx RZ*RY'*URzy URzz] | |
# BUILD THE WAVE MATRIX A | |
CE = [ZZ -sz_hx*DEZ sy_hx*DEY | |
sz_hy*DEZ ZZ -sx_hy*DEX | |
-sy_hz*DEY sx_hz*DEX ZZ] | |
CH = [ZZ -sz_ex*DHZ sy_ex*DHY | |
sz_ey*DHZ ZZ -sx_ey*DHX | |
-sy_ez*DHY sx_ez*DHX ZZ] | |
A = CH * (UR \ CE) - ER | |
# CALCULATE POLARIZATION VECTOR P | |
az = [0.0, 0.0, 1.0] | |
if abs(src.θ) < 1e-6 | |
ate = [0.0, 1.0, 0.0] | |
else | |
ate = az × kinc | |
ate /= norm(ate) | |
end | |
atm = kinc × ate | |
atm /= norm(atm) | |
P = src.pte * ate + src.ptm * atm | |
# CALCULATE SOURCE FIELD fsrc | |
fphase = [cis(-(kinc[1]*x + kinc[2]*y + kinc[3]*z)) for y in ya2, x in xa2, z in za2] | |
fx = P[1] * fphase[2:2:Nx2, 1:2:Ny2, 1:2:Nz2] | |
fy = P[2] * fphase[1:2:Nx2, 2:2:Ny2, 1:2:Nz2] | |
fz = P[3] * fphase[1:2:Nx2, 1:2:Ny2, 2:2:Nz2] | |
fsrc = [vec(fx); vec(fy); vec(fz)] | |
# CALCULATE SCATERED FIELD MASKING MATRIX | |
nz = dev.NPML[1] + 3 | |
Q = zeros(Nx, Ny, Nz) | |
Q[:, :, 1:nz] .= 1 | |
Q = Diagonal(repeat(vec(Q), 3)) | |
#Q = diagonalize(Q) | |
#Q = [Q ZZ ZZ; ZZ Q ZZ; ZZ ZZ Q] | |
# CALCULATE SOURCE VECTOR B | |
QAAQ = Q*A - A*Q | |
b = QAAQ * fsrc | |
@show typeof.((QAAQ, fsrc, b)) | |
@show size.((QAAQ, fsrc, b)) | |
@show nnz(QAAQ) | |
@show maximum(abs ∘ imag, QAAQ) | |
@show maximum(abs, QAAQ - real(QAAQ)) | |
@show maximum(abs, b - real(QAAQ) * fsrc) | |
#= | |
# COMPUTE JACOBI PRECONDITIONER | |
Pl = Diagonal(inv.(Array(diag(A)))) | |
# SOLVE FOR FIELD | |
Agpu = CuSparseMatrixCSC(A) | |
bgpu = CuArray(b) | |
Plgpu = CuArray(Pl) | |
f, hist = idrs(A, b; Pl=Pl, log=true, verbose=true, maxiter=2000, abstol=1e-8, reltol=1e-6) | |
@time fgpu, histgpu = idrs(Agpu, bgpu; Pl=Plgpu, log=true, maxiter=2000, abstol=1e-8, reltol=1e-6) | |
#f = bicg(A,b,tol,maxit,[],[],f); | |
#f = gather(f); | |
# EXTRACT FIELD COMPONENTS | |
fx = reshape(f[1:M], Nx, Ny, Nz) | |
fy = reshape(f[M+1:2*M], Nx, Ny, Nz) | |
fz = reshape(f[2*M+1:3*M], Nx, Ny, Nz) | |
######################################################### | |
## ANALYZE REFLECTION AND TRANSMISSION | |
######################################################### | |
# EXTRACT FIELD FROM RECORD PLANES | |
nz_ref = dev.NPML[1] + 2 | |
fx_ref = fx[:, :, nz_ref] | |
fy_ref = fy[:, :, nz_ref] | |
fz_ref = fz(:, :, nz_ref-1:nz_ref) | |
nz_trn = Nz - dev.NPML[1] | |
fx_trn = fx[:, :, nz_trn] | |
fy_trn = fy[:, :, nz_trn] | |
fz_trn = fz[:, :, nz_trn-1:nz_trn] | |
# INTERPOLATE FIELDS TO ORIGIN OF YEE CELLS | |
Ex_ref = copy(fx_ref) | |
Ex_ref[1, :] = (fx_ref[Nx, :] + fx_ref[1, :]) / 2 | |
Ex_ref[2:Nx, :] = (fx_ref[1:Nx-1, :] + fx_ref[2:Nx, :]) / 2 | |
Ey_ref = copy(fy_ref) | |
Ey_ref[:,1] = (fy_ref[:, Ny] + fy_ref[:, 1]) / 2 | |
Ey_ref[:, 2:Ny] = (fy_ref[:, 1:Ny-1] + fy_ref[:, 2:Ny]) / 2 | |
Ez_ref = (fz_ref[:,:,1] + fz_ref[:, :, 2]) / 2 | |
Ex_trn = copy(fx_trn) | |
Ex_trn[1,:] = (fx_trn[Nx, :] + fx_trn[1, :]) / 2 | |
Ex_trn[2:Nx, :] = (fx_trn[1:Nx-1, :] + fx_trn[2:Nx, :]) / 2 | |
Ey_trn = copy(fy_trn) | |
Ey_trn[:, 1] = (fy_trn[:, Ny] + fy_trn[:, 1]) / 2 | |
Ey_trn[:, 2:Ny] = (fy_trn[:, 1:Ny-1] + fy_trn[:, 2:Ny]) / 2 | |
Ez_trn = (fz_trn[:, :, 1] + fz_trn[:, :, 2]) / 2 | |
# REMOVE PHASE TILT | |
Ex_ref ./= fphase[2:2:Nx2, 1:2:Ny2, 1] | |
Ey_ref ./= fphase[1:2:Nx2, 2:2:Ny2, 1] | |
Ez_ref ./= fphase[1:2:Nx2, 1:2:Ny2, 2] | |
Ex_trn ./= fphase[2:2:Nx2, 1:2:Ny2, 1] | |
Ey_trn ./= fphase[1:2:Nx2, 2:2:Ny2, 1] | |
Ez_trn ./= fphase[1:2:Nx2, 1:2:Ny2, 2] | |
# CALCULATE AMPLITUDES OF DIFFRACTION ORDERS | |
NxNy = Nx * Ny | |
ax_ref = fftshift(fft(Ex_ref)) / NxNy | |
ay_ref = fftshift(fft(Ey_ref)) / NxNy | |
az_ref = fftshift(fft(Ez_ref)) / NxNy | |
a_ref = abs2.(ax_ref) + abs2.(ay_ref) + abs2.(az_ref) | |
s11 = sqrt.(ax_ref.^2 + ay_ref.^2 + az_ref.^2) | |
ax_trn = fftshift(fft(Ex_trn)) / NxNy | |
ay_trn = fftshift(fft(Ey_trn)) / NxNy | |
az_trn = fftshift(fft(Ez_trn)) / NxNy | |
a_trn = abs2.(ax_trn) + abs2.(ay_trn) + abs2.(az_trn) | |
s21 = sqrt(ax_trn.^2 + ay_trn.^2 + az_trn.^2) | |
# CALCULATE DIFFRACTION EFFICIENCIES | |
RDE = a_ref .* real(kz_ref / kinc[3]) | |
TDE = a_trn .* real(ur_ref / ur_trn * kz_trn / kinc[3]) | |
# Calculate Overall Response | |
REF = sum(RDE) | |
TRN = sum(TDE) | |
CON = REF + TRN; | |
return (; fx, fy, fz, m, n, RDE, REF, TDE, TRN, s11, s21, CON) | |
=# | |
end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using LinearAlgebra | |
using SparseArrays | |
#using MKL, MKLSparse | |
using Statistics: mean | |
include("yeeder3d.jl") | |
include("calcpml3d.jl") | |
include("fdfd3d.jl") | |
# DEFINE UNITS | |
const meters = 1 | |
const centimeters = 1e-2meters | |
const seconds = 1 | |
const hertz = 1 / seconds | |
const gigahertz = 1e9 * hertz | |
const degrees = π / 180 | |
# CONSTANTS | |
const c0 = 299792458 * meters/seconds | |
######################################################### | |
## DASHBOARD | |
######################################################### | |
# DEFINE SOURCE PARAMETERS | |
f1 = 4.5gigahertz | |
f2 = 5.5gigahertz | |
NFREQ = 500 | |
FREQ = LinRange(f1, f2, NFREQ) | |
θ = ϕ = 0degrees | |
pte = 1 / sqrt(2) | |
ptm = im / sqrt(2) | |
# DEFINE GMRF PARAMETERS | |
er1 = 1.0 | |
er2 = 2.0 | |
erL = 3.0 | |
erH = 5.0 | |
a = 3.7centimeters | |
t = 1.5centimeters | |
frac = 0.5 | |
# DEFINE GRID PARAMETERS | |
NRES = 20 | |
NPML = [10, 10] | |
ermax = max(er1, er2, erL, erH) | |
nmax = sqrt(ermax) | |
lam_max = c0 / minimum(FREQ) | |
SPACER = lam_max * [1, 1] | |
######################################################### | |
## CALCULATE OPTIMIZED GRID | |
######################################################### | |
# CALCULATE PRELIMINARY RESOLUTION | |
lam_min = c0 / maximum(FREQ) | |
dx = dy = dz = lam_min/nmax/NRES | |
# SNAP GRID TO CRITICAL DIMENSIONS | |
snap(w, Δ) = begin n = ceil(Int, w/Δ); δ = w/n; return (n, δ) end | |
nx, dx = snap(a, dx) | |
ny, dy = snap(a, dy) | |
nz, dz = snap(t, dz) | |
# CALCULATE GRID SIZE | |
Nx = ceil(Int, a / dx) | |
Sx = Nx * dx | |
Ny = ceil(Int, a / dy) | |
Sy = Ny * dy | |
Sz = SPACER[1] + t + SPACER[2] | |
Nz = NPML[1] + ceil(Int, Sz / dz) + NPML[2] | |
Sz = Nz * dz | |
# 2x GRID | |
Nx2, dx2 = (2Nx, dx/2) | |
Ny2, dy2 = (2Ny, dy/2) | |
Nz2, dz2 = (2Nz, dz/2) | |
# GRID AXES | |
xa = [1:Nx;] * dx; xa .-= mean(xa) | |
ya = [1:Ny;] * dy; ya .-= mean(ya) | |
za = [0.5:Nz-0.5;] * dz | |
xa2 = [1:Nx2;] * dx2; xa2 .-= mean(xa2) | |
ya2 = [1:Ny2;] * dy2; ya2 .-= mean(ya2) | |
za2 = [0.5:Nz2-0.5;] * dz2 | |
# MESHGRIDS | |
#[Y,X,Z] = meshgrid(ya,xa,za); | |
#[Y2,X2,Z2] = meshgrid(ya2,xa2,za2); | |
######################################################### | |
## BUILD MATERIALS ARRAYS | |
######################################################### | |
# INITIALIZE TO FREE SPACE | |
ER2 = ones(Nx2, Ny2, Nz2) | |
UR2 = ones(Nx2, Ny2, Nz2) | |
# BUILD GRATING LAYER | |
r = a * sqrt(frac / π) | |
r² = r * r | |
ER2 = [x^2 + y^2 ≤ r² ? erH : erL for y in ya2, x in xa2, z in za2] | |
# ADD SUPERSTRATE AND SUBSTRATE | |
nz1 = 2 * NPML[1] + round(Int, SPACER[1] / dz2) + 1 | |
nz2 = nz1 + round(Int, t / dz2) - 1 | |
ER2[:, :, 1:nz1-1] .= er1 | |
ER2[:, :, nz2+1:Nz2] .= er2; | |
######################################################### | |
## PERFORM FREQUENCY SWEEP | |
######################################################### | |
# dev and src | |
RES = [dx, dy, dz] | |
dev = (; NPML, ER2, UR2, RES) | |
f0 = f1 | |
λ₀ = c0 / f0 | |
src = (; θ, ϕ, pte, ptm, λ₀) | |
# Call FDFD3D() | |
dat = fdfd3d(dev, src) | |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using LinearAlgebra, SparseArrays | |
""" | |
yeeder3d(NS,RES,BC,kinc) -> (DEX,DEY,DEZ,DHX,DHY,DHZ) | |
Compute derivative matrices on a 3D Yee grid | |
# Input Arguments | |
- `NS`: `(Nx Ny Nz)` Grid size. | |
- `RES`: `(dx, dy, dz)` Grid resolution. For normalized grids contains `(k0*dx, k0*dy, k0*dz)`. | |
- `BC`: `(xbc, ybc, zbc)` Boundary conditions, where 0(1) implies Dirichlet(periodic) boundary conditions. | |
- `kinc`: `[kx ky kz]` Incident wave vector (needed only for PBCs). For normalized grids contains `kinc/k0`. | |
# Return Values | |
- `DEX`: Derivative matrix wrt x for electric fields | |
- `DEY`: Derivative matrix wrt y for electric fields | |
- `DEZ`: Derivative matrix wrt z for electric fields | |
- `DHX`: Derivative matrix wrt x for magnetic fields | |
- `DHY`: Derivative matrix wrt y for magnetic fields | |
- `DHZ`: Derivative matrix wrt z for magnetic fields | |
""" | |
function yeeder3d(NS, RES, BC, kinc=zeros(3)) | |
######################################################### | |
## HANDLE INPUT ARGUMENTS | |
######################################################### | |
# EXTRACT GRID PARAMETERS | |
Nx, Ny, Nz = NS | |
dx, dy, dz = RES | |
# DETERMINE MATRIX SIZE | |
M = Nx * Ny * Nz | |
# Determine if real or complex matrices are needed | |
if (Nx == 1 && kinc[1] ≠ 0) || (Ny == 1 && kinc[2] ≠ 0) || (Nz == 1 && kinc[3] ≠ 0) || any(≠(0), BC) | |
T = ComplexF64 | |
else | |
T = Float64 | |
end | |
######################################################### | |
## BUILD DEX | |
######################################################### | |
if Nx==1 | |
# HANDLE IF SIZE IS 1 CELL | |
DEX = sparse(-im * kinc[1] * I, M, M) | |
else | |
# HANDLE ALL OTHER CASES | |
d0 = fill(T(-1/dx), M) # Center Diagonal | |
d1 = fill(T(1/dx), M) # Upper Diagonal | |
d1[Nx+1:Nx:M] .= zero(T) | |
d1 = d1[2:end] | |
# Build Derivative Matrix with Dirichlet BCs | |
DEX = spdiagm(M, M, 0 => d0, 1 => d1) | |
# Incorporate Periodic Boundary Conditions | |
if BC[1] == 1 | |
d1 = zeros(T, M - (Nx-1)) | |
d1[1:Nx:M] .= cis(-kinc[1] * (Nx*dx)) / dx | |
DEX += spdiagm(M, M, 1-Nx => d1) | |
end | |
end | |
######################################################### | |
## BUILD DEY | |
######################################################### | |
# HANDLE IF SIZE IS 1 CELL | |
if Ny==1 | |
DEY = sparse(-im * kinc[2] * I, M, M) | |
else | |
# HANDLE ALL OTHER CASES | |
d0 = fill(T(-1/dy), M) # Center Diagonal | |
# Upper Diagonal | |
d1 = [fill(T(1/dy), (Ny-1)*Nx); zeros(T, Nx)] | |
d1 = repeat(d1, Nz-1) | |
d1 = [zeros(T, Nx); d1; fill(T(1/dy), (Ny-1)*Nx)] | |
d1 = d1[Nx+1:end] | |
# Build Derivative Matrix with Dirichlet BCs | |
DEY = spdiagm(M, M, 0 => d0, Nx => d1) | |
# Incorporate Periodic Boundary Conditions | |
if BC[2] == 1 | |
ph = cis(-kinc[2] * (Ny*dy)) / dy | |
d1 = [fill(ph, Nx); zeros(T, (Ny-1)*Nx)] | |
d1 = repeat(d1, Nz-1) | |
d1 = [d1; fill(ph, (Ny-1)*Nx); zeros(T, Nx)] | |
d1 = d1[1:M-Nx*(Ny-1)] | |
DEY += spdiagm(M, M, -Nx*(Ny-1) => d1) | |
end | |
end | |
######################################################### | |
## BUILD DEZ | |
######################################################### | |
# HANDLE IF SIZE IS 1 CELL | |
if Nz==1 | |
DEZ = sparse(-im * kinc[3] * I, M, M) | |
else | |
# HANDLE ALL OTHER CASES | |
d0 = fill(T(1/dz), M) # Center Diagonal | |
# Build Derivative Matrix with Dirichlet BCs | |
DEZ = spdiagm(M, M, 0 => -d0, Nx*Ny => d0[Nx*Ny+1:end]) | |
# Incorporate Periodic Boundary Conditions | |
if BC[3] == 1 | |
d0 = fill(cis(-kinc[3] * (Nz*dz)) / dz, M - Nx*Ny*(Nz-1)) | |
DEZ += spdiagm(M, M, -Nx*Ny*(Nz-1) => d0) | |
end | |
end | |
######################################################### | |
## BUILD DHX, DHY AND DHZ | |
######################################################### | |
DHX = sparse(-DEX') | |
DHY = sparse(-DEY') | |
DHZ = sparse(-DEZ') | |
return (DEX, DEY, DEZ, DHX, DHY, DHZ) | |
end # function |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment