From d72169bbf4ca7375a7ec9dc0ec4fe4ecf6c911ff Mon Sep 17 00:00:00 2001 From: hochunlin Date: Sat, 16 Dec 2023 19:45:34 -0800 Subject: [PATCH] Update v0.4.0 Former-commit-id: 98cb3b6d09417e5b34035a97731376ed34bad97d --- .../README.md | 68 +++++++----- .../build_epsilon_disorder.jl | 0 ...ide_disorder_with_phase_conjugation.ipynb} | 48 +++++---- ...inside_disorder_with_phase_conjugation.jl} | 60 ++++++----- .../intensity_comparison.png | Bin ...ase_conjugated_focusing.gif.REMOVED.git-id | 0 .../regular_focusing.gif.REMOVED.git-id | 0 .../README.md | 90 ++++++++++++++-- .../open_channel_through_disorder.ipynb | 102 ++++++++++++++++-- .../open_channel_through_disorder.jl | 83 +++++++++++--- .../README.md | 2 +- .../reflection_matrix_Gaussian_beams.ipynb | 2 +- .../reflection_matrix_Gaussian_beams.jl | 2 +- src/MESTI.jl | 6 +- src/mesti2s.jl | 44 ++++---- src/mesti_build_channels.jl | 24 ++--- src/mesti_build_fdfd_matrix.jl | 2 +- ... => mesti_build_transverse_function_1d.jl} | 10 +- src/mesti_main.jl | 28 ++--- src/mesti_matrix_solver.jl | 7 +- ...mal_PML.jl => mesti_optimal_pml_params.jl} | 6 +- ...tudinal.jl => mesti_setup_longitudinal.jl} | 6 +- src/mesti_subpixel_smoothing.jl | 21 ++-- src/mumps3_interface.jl | 32 ++++-- src/mumps3_printing.jl | 25 ++++- test/interface_t_r_test.jl | 2 +- test/unitary_test.jl | 2 +- 27 files changed, 482 insertions(+), 190 deletions(-) rename examples/{2d_focusing_phase_conjugated_light_through_disorder => 2d_focusing_inside_disorder_with_phase_conjugation}/README.md (78%) rename examples/{2d_focusing_phase_conjugated_light_through_disorder => 2d_focusing_inside_disorder_with_phase_conjugation}/build_epsilon_disorder.jl (100%) rename examples/{2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.ipynb => 2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.ipynb} (83%) rename examples/{2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.jl => 2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.jl} (78%) rename examples/{2d_focusing_phase_conjugated_light_through_disorder => 2d_focusing_inside_disorder_with_phase_conjugation}/intensity_comparison.png (100%) rename examples/{2d_focusing_phase_conjugated_light_through_disorder => 2d_focusing_inside_disorder_with_phase_conjugation}/phase_conjugated_focusing.gif.REMOVED.git-id (100%) rename examples/{2d_focusing_phase_conjugated_light_through_disorder => 2d_focusing_inside_disorder_with_phase_conjugation}/regular_focusing.gif.REMOVED.git-id (100%) rename src/{build_transverse_function_1d.jl => mesti_build_transverse_function_1d.jl} (95%) rename src/{get_optimal_PML.jl => mesti_optimal_pml_params.jl} (89%) rename src/{setup_longitudinal.jl => mesti_setup_longitudinal.jl} (94%) diff --git a/examples/2d_focusing_phase_conjugated_light_through_disorder/README.md b/examples/2d_focusing_inside_disorder_with_phase_conjugation/README.md similarity index 78% rename from examples/2d_focusing_phase_conjugated_light_through_disorder/README.md rename to examples/2d_focusing_inside_disorder_with_phase_conjugation/README.md index 09a4c92..a7ddc0d 100644 --- a/examples/2d_focusing_phase_conjugated_light_through_disorder/README.md +++ b/examples/2d_focusing_inside_disorder_with_phase_conjugation/README.md @@ -1,6 +1,6 @@ -# Focusing Phase Conjugated Light Through Disorder +# Focusing Inside Disorder With Phase Conjugation -In this example, we show how to use mesti() to compute the field profile of a point source in a scattering disordered medium, do phase conjugation to determine an incident wavefront that can focus on the the disorder, and then use mesti2s() again to compute the field profile to show its focus. +In this example, we show how to use mesti() to project field generated from a point source inside disorder onto propagating channels through APF method, do phase conjugation to determine an incident wavefront that can focus inside the disorder, and then use mesti2s() again to compute the field profile to show its focus. ```julia # Call necessary packages @@ -43,7 +43,7 @@ build_epsilon_disorder(W, L, r_min, r_max, min_sep, no_scatterer_center = true) ``` -# Compute the field profile of a point source +# Projecting field generated from a point source onto propagating channels through APF ```julia syst = Syst() @@ -63,28 +63,48 @@ Bx.pos = [[m0_focus,l0_focus+pml_npixels+1,1,1]] Bx.data = [ones(1,1)] # put PML along z-direction -pml = get_optimal_PML(syst.wavelength/syst.dx) +pml = mesti_optimal_pml_params(syst.wavelength/syst.dx) pml.npixels = pml_npixels pml.direction = "z" syst.PML = [pml] -# field profile: input from a point source -Ex_field, _ = mesti(syst, [Bx]) +# build channels for the low side (air) +channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low) +N_prop_low = channels_low.N_prop # number of propagating channels on the low side + +# build projection profiles C on the low side +Cx = Source_struct() +ny_Ex = Int(W/dx) +Cx.pos = [[1,pml_npixels+1,ny_Ex,1]] +# sqrt(nu)*conj(u_Ex(m)) serves as a projection profile for a propagating channel, where nu = sin(kzdx) +# here we project the result onto all propagating channels as a output basis +C_low = (conj(channels_low.u_x_m(channels_low.kydx_prop))).*reshape(channels_low.sqrt_nu_prop,1,:).*reshape(exp.((-1im*0.5)*channels_low.kzdx_prop),1,:) # 0.5 pixel backpropagation indicates that the projection(detect) region is half a pixel away from z = 0 +Cx.data = [C_low] + +# Calculate output channel amplitude coefficients (w) for the propagating channels from the point source through APF +w, _ = mesti(syst, [Bx], [Cx]) + +# We can also compute the field profile Ex_field and project the field on the detect region onto the projection profiles. +# It would generate the same output channel amplitude coefficients +# Ex_field, _ = mesti(syst, [Bx]) +# C = transpose(C_low) +# w = C_low*Ex_field[:,pml_npixels+1] ``` ```text:Output ===System size=== ny_Ex = 5400; nz_Ex = 1381 for Ex(y,z) UPML on -z +z sides; ; yBC = periodic; zBC = PEC -Building B,C... elapsed time: 0.230 secs -Building A ... elapsed time: 9.228 secs -< Method: factorize_and_solve using MUMPS in single precision with AMD ordering > -Analyzing ... elapsed time: 3.977 secs -Factorizing ... elapsed time: 64.983 secs -Solving ... elapsed time: 4.065 secs - Total elapsed time: 76.477 secs +Building B,C... elapsed time: 1.664 secs +Building A ... elapsed time: 8.106 secs +< Method: APF using MUMPS in single precision with AMD ordering > +Building K ... elapsed time: 2.822 secs +false +Analyzing ... elapsed time: 5.076 secs +Factorizing ... elapsed time: 76.478 secs + Total elapsed time: 99.805 secs ``` -# Compute the regular focusing and phase-conjugated focusing profile +# Compute the regular focusing and phase-conjugated focusing profiles ```julia # Specify the system for mesti2s() and mesti_build_channels() @@ -101,19 +121,12 @@ syst.zPML = [pml] # equivalent average epsilon for this disordered system epsilon_ave = mean(epsilon) -# build channels for the equivalent average epsilon and low side (air) +# build channels for the equivalent average epsilon channels_ave_epsilon = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_ave) -channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low) N_prop_ave_epsilon = channels_ave_epsilon.N_prop # number of propagating channels on the equivalent average epsilon -N_prop_low = channels_low.N_prop # number of propagating channels on the low side -dn = 0.5 # regular focus wavefront -wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-dn)) - -# build projection matrix C on the low side -C_low = channels_low.sqrt_nu_prop.*exp.((-1im*dn)*channels_low.kzdx_prop).*convert(Matrix, adjoint(channels_low.u_x_m(channels_low.kydx_prop))) -proj_coefficient = C_low*Ex_field[:,pml_npixels+1] +wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-0.5)) # 0.5 pixel indicates that the source is half a pixel away from z = 0 # specify two input incident wavefronts: # (1) regular focusing wavefront @@ -123,12 +136,11 @@ input.v_low = zeros(ComplexF64, N_prop_low, 2) input.v_low[:, 1] = wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]/norm(wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]) # for the phased conjugated input: -# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky) +# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky) # perm() means permute a vector that switches one propagating channel with one # having a complex-conjugated transverse profile -# for periodic boundary, this flips the sign of ky. -input.v_low[:, 2] = conj(proj_coefficient)[channels_low.ind_prop_conj]/norm(proj_coefficient) - +# for the periodic boundary, this flips the sign of ky. +input.v_low[:, 2] = conj(w)[channels_low.ind_prop_conj]/norm(w) # we will also get the field profile in the free spaces on the two sides, for # plotting purpose. @@ -252,4 +264,4 @@ println("I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = ", @sprintf("%d", round(a ``` ```text:Output I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = 1800 -``` +``` \ No newline at end of file diff --git a/examples/2d_focusing_phase_conjugated_light_through_disorder/build_epsilon_disorder.jl b/examples/2d_focusing_inside_disorder_with_phase_conjugation/build_epsilon_disorder.jl similarity index 100% rename from examples/2d_focusing_phase_conjugated_light_through_disorder/build_epsilon_disorder.jl rename to examples/2d_focusing_inside_disorder_with_phase_conjugation/build_epsilon_disorder.jl diff --git a/examples/2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.ipynb b/examples/2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.ipynb similarity index 83% rename from examples/2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.ipynb rename to examples/2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.ipynb index 21b2ea1..57f7c9e 100644 --- a/examples/2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.ipynb +++ b/examples/2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.ipynb @@ -5,7 +5,7 @@ "id": "ae0e7aae-9838-44df-9330-8a1d8f9fc85d", "metadata": {}, "source": [ - "# Focusing Phase Conjugated Light Through Disorder" + "# Focusing Inside Disorder With Phase Conjugation" ] }, { @@ -13,7 +13,7 @@ "id": "96a5a5eb-94ba-49be-9f3f-3e6d71cea482", "metadata": {}, "source": [ - "In this example, we show how to use mesti() to compute the field profile of a point source in a scattering disordered medium, do phase conjugation to determine an incident wavefront that can focus on the the disorder, and then use mesti2s() again to compute the field profile to show its focus." + "In this example, we show how to use mesti() to project field generated from a point source inside disorder onto propagating channels through APF method, do phase conjugation to determine an incident wavefront that can focus inside the disorder, and then use mesti2s() again to compute the field profile to show its focus." ] }, { @@ -80,7 +80,7 @@ "id": "be175504-fd46-4922-9262-527fa6389963", "metadata": {}, "source": [ - "# Compute the field profile of a point source " + "# Projecting field generated from a point source onto propagating channels through APF" ] }, { @@ -107,13 +107,32 @@ "Bx.data = [ones(1,1)]\n", "\n", "# put PML along z-direction\n", - "pml = get_optimal_PML(syst.wavelength/syst.dx)\n", + "pml = mesti_optimal_pml_params(syst.wavelength/syst.dx)\n", "pml.npixels = pml_npixels\n", "pml.direction = \"z\"\n", "syst.PML = [pml]\n", "\n", - "# field profile: input from a point source\n", - "Ex_field, _ = mesti(syst, [Bx]);" + "# build channels for the low side (air)\n", + "channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)\n", + "N_prop_low = channels_low.N_prop # number of propagating channels on the low side \n", + "\n", + "# build projection profiles C on the low side\n", + "Cx = Source_struct()\n", + "ny_Ex = Int(W/dx)\n", + "Cx.pos = [[1,pml_npixels+1,ny_Ex,1]]\n", + "# sqrt(nu)*conj(u_Ex(m)) serves as a projection profile for a propagating channel, where nu = sin(kzdx)\n", + "# here we project the result onto all propagating channels as a output basis \n", + "C_low = (conj(channels_low.u_x_m(channels_low.kydx_prop))).*reshape(channels_low.sqrt_nu_prop,1,:).*reshape(exp.((-1im*0.5)*channels_low.kzdx_prop),1,:) # 0.5 pixel backpropagation indicates that the projection(detect) region is half a pixel away from z = 0\n", + "Cx.data = [C_low]\n", + "\n", + "# Calculate output channel amplitude coefficients (w) for the propagating channels from the point source through APF\n", + "w, _ = mesti(syst, [Bx], [Cx])\n", + "\n", + "# We can also compute the field profile Ex_field and project the field on the detect region onto the projection profiles. \n", + "# It would generate the same output channel amplitude coefficients\n", + "# Ex_field, _ = mesti(syst, [Bx])\n", + "# C = transpose(C_low)\n", + "# w = C_low*Ex_field[:,pml_npixels+1]" ] }, { @@ -121,7 +140,7 @@ "id": "21c79756-b0be-4e62-b998-ce405190a866", "metadata": {}, "source": [ - "# Compute the field profile of a point source " + "# Compute the regular focusing and phase-conjugated focusing profiles" ] }, { @@ -145,19 +164,12 @@ "# equivalent average epsilon for this disordered system\n", "epsilon_ave = mean(epsilon)\n", "\n", - "# build channels for the equivalent average epsilon and low side (air)\n", + "# build channels for the equivalent average epsilon\n", "channels_ave_epsilon = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_ave)\n", - "channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low)\n", "N_prop_ave_epsilon = channels_ave_epsilon.N_prop # number of propagating channels on the equivalent average epsilon \n", - "N_prop_low = channels_low.N_prop # number of propagating channels on the low side \n", "\n", - "dn = 0.5\n", "# regular focus wavefront\n", - "wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-dn))\n", - "\n", - "# build projection matrix C on the low side\n", - "C_low = channels_low.sqrt_nu_prop.*exp.((-1im*dn)*channels_low.kzdx_prop).*convert(Matrix, adjoint(channels_low.u_x_m(channels_low.kydx_prop)))\n", - "proj_coefficient = C_low*Ex_field[:,pml_npixels+1]\n", + "wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-0.5)) # 0.5 pixel indicates that the source is half a pixel away from z = 0\n", "\n", "# specify two input incident wavefronts:\n", "# (1) regular focusing wavefront\n", @@ -167,11 +179,11 @@ "input.v_low[:, 1] = wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]/norm(wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)])\n", "\n", "# for the phased conjugated input: \n", - "# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)\n", + "# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky)\n", "# perm() means permute a vector that switches one propagating channel with one\n", "# having a complex-conjugated transverse profile\n", "# for the periodic boundary, this flips the sign of ky. \n", - "input.v_low[:, 2] = conj(proj_coefficient)[channels_low.ind_prop_conj]/norm(proj_coefficient)\n", + "input.v_low[:, 2] = conj(w)[channels_low.ind_prop_conj]/norm(w)\n", "\n", "# we will also get the field profile in the free spaces on the two sides, for\n", "# plotting purpose.\n", diff --git a/examples/2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.jl b/examples/2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.jl similarity index 78% rename from examples/2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.jl rename to examples/2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.jl index 070ccd9..f69e7b3 100644 --- a/examples/2d_focusing_phase_conjugated_light_through_disorder/focusing_phase_conjugated_light_through_disorder.jl +++ b/examples/2d_focusing_inside_disorder_with_phase_conjugation/focusing_inside_disorder_with_phase_conjugation.jl @@ -1,9 +1,10 @@ -## Focusing Phase Conjugated Light Through Disorder +## Focusing Inside Disorder With Phase Conjugation -# In this example, we show how to use mesti() to compute the field profile of -# a point source in a scattering disordered medium, do phase conjugation to -# determine an incident wavefront that can focus on the the disorder, and then -# use mesti2s() again to compute the field profile to show its focus. +# In this example, we show how to use mesti() to project field generated from +# a point source inside disorder onto propagating channels through APF method, +# do phase conjugation to determine an incident wavefront that can focus inside +# the disorder, and then use mesti2s() again to compute the field profile to +# show its focus. # Call necessary packages using MESTI, GeometryPrimitives, LinearAlgebra, Statistics, Printf @@ -42,7 +43,7 @@ build_epsilon_disorder(W, L, r_min, r_max, min_sep, epsilon_scat, epsilon_bg, build_TM; no_scatterer_center = true) -## Compute the field profile of a point source +## Projecting field generated from a point source onto propagating channels through APF syst = Syst() pml_npixels = 15 syst.length_unit = "lambda_0" @@ -60,15 +61,34 @@ Bx.pos = [[m0_focus,l0_focus+pml_npixels+1,1,1]] Bx.data = [ones(1,1)] # put PML along z-direction -pml = get_optimal_PML(syst.wavelength/syst.dx) +pml = mesti_optimal_pml_params(syst.wavelength/syst.dx) pml.npixels = pml_npixels pml.direction = "z" syst.PML = [pml] -# field profile: input from a point source -Ex_field, _ = mesti(syst, [Bx]) - -## Compute the regular focusing and phase-conjugated focusing profile +# build channels for the low side (air) +channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low) +N_prop_low = channels_low.N_prop # number of propagating channels on the low side + +# build projection profiles C on the low side +Cx = Source_struct() +ny_Ex = Int(W/dx) +Cx.pos = [[1,pml_npixels+1,ny_Ex,1]] +# sqrt(nu)*conj(u_Ex(m)) serves as a projection profile for a propagating channel, where nu = sin(kzdx) +# here we project the result onto all propagating channels as a output basis +C_low = (conj(channels_low.u_x_m(channels_low.kydx_prop))).*reshape(channels_low.sqrt_nu_prop,1,:).*reshape(exp.((-1im*0.5)*channels_low.kzdx_prop),1,:) # 0.5 pixel backpropagation indicates that the projection(detect) region is half a pixel away from z = 0 +Cx.data = [C_low] + +# Calculate output channel amplitude coefficients (w) for the propagating channels from the point source through APF +w, _ = mesti(syst, [Bx], [Cx]) + +# We can also compute the field profile Ex_field and project the field on the detect region onto the projection profiles. +# It would generate the same output channel amplitude coefficients +# Ex_field, _ = mesti(syst, [Bx]) +# C = transpose(C_low) +# w = C_low*Ex_field[:,pml_npixels+1] + +## Compute the regular focusing and phase-conjugated focusing profiles # Specify the system for mesti2s() and mesti_build_channels() syst = Syst() syst.epsilon_xx = epsilon @@ -83,19 +103,12 @@ syst.zPML = [pml] # equivalent average epsilon for this disordered system epsilon_ave = mean(epsilon) -# build channels for the equivalent average epsilon and low side (air) +# build channels for the equivalent average epsilon channels_ave_epsilon = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_ave) -channels_low = mesti_build_channels(Int(W/dx), yBC, 2*pi*syst.wavelength*dx, epsilon_low) N_prop_ave_epsilon = channels_ave_epsilon.N_prop # number of propagating channels on the equivalent average epsilon -N_prop_low = channels_low.N_prop # number of propagating channels on the low side -dn = 0.5 # regular focus wavefront -wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-dn)) - -# build projection matrix C on the low side -C_low = channels_low.sqrt_nu_prop.*exp.((-1im*dn)*channels_low.kzdx_prop).*convert(Matrix, adjoint(channels_low.u_x_m(channels_low.kydx_prop))) -proj_coefficient = C_low*Ex_field[:,pml_npixels+1] +wf_reg_focus = exp.(-1im*channels_ave_epsilon.kydx_prop*(m0_focus)) .* exp.(-1im*channels_ave_epsilon.kzdx_prop*(l0_focus-0.5)) # 0.5 pixel indicates that the source is half a pixel away from z = 0 # specify two input incident wavefronts: # (1) regular focusing wavefront @@ -105,12 +118,11 @@ input.v_low = zeros(ComplexF64, N_prop_low, 2) input.v_low[:, 1] = wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]/norm(wf_reg_focus[Int((N_prop_ave_epsilon-N_prop_low)/2+1):Int(end-(N_prop_ave_epsilon-N_prop_low)/2)]) # for the phased conjugated input: -# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky) +# conj(coefficient*u) = conj(coefficient)*conj(u) = conj(coefficient)*perm(u(ky)) = perm(conj(coefficient))*u(ky) # perm() means permute a vector that switches one propagating channel with one # having a complex-conjugated transverse profile # for the periodic boundary, this flips the sign of ky. -input.v_low[:, 2] = conj(proj_coefficient)[channels_low.ind_prop_conj]/norm(proj_coefficient) - +input.v_low[:, 2] = conj(w)[channels_low.ind_prop_conj]/norm(w) # we will also get the field profile in the free spaces on the two sides, for # plotting purpose. @@ -202,4 +214,4 @@ intensity_plot = plot(plt3, plt4, layout = @layout([a b]), size=(800, 400)) display(intensity_plot) # compare the ratio of intensity on the focus point -println("I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = ", @sprintf("%d", round(abs.(Ex[m0_focus,opts.nz_low+l0_focus,2]).^2/abs.(Ex[m0_focus,opts.nz_low+l0_focus,1]).^2, digits=-2))) \ No newline at end of file +println("I_phase_congugation(y_0,z_0)/I_reg(y_0,z_0) = ", @sprintf("%d", round(abs.(Ex[m0_focus,opts.nz_low+l0_focus,2]).^2/abs.(Ex[m0_focus,opts.nz_low+l0_focus,1]).^2, digits=-2))) diff --git a/examples/2d_focusing_phase_conjugated_light_through_disorder/intensity_comparison.png b/examples/2d_focusing_inside_disorder_with_phase_conjugation/intensity_comparison.png similarity index 100% rename from examples/2d_focusing_phase_conjugated_light_through_disorder/intensity_comparison.png rename to examples/2d_focusing_inside_disorder_with_phase_conjugation/intensity_comparison.png diff --git a/examples/2d_focusing_phase_conjugated_light_through_disorder/phase_conjugated_focusing.gif.REMOVED.git-id b/examples/2d_focusing_inside_disorder_with_phase_conjugation/phase_conjugated_focusing.gif.REMOVED.git-id similarity index 100% rename from examples/2d_focusing_phase_conjugated_light_through_disorder/phase_conjugated_focusing.gif.REMOVED.git-id rename to examples/2d_focusing_inside_disorder_with_phase_conjugation/phase_conjugated_focusing.gif.REMOVED.git-id diff --git a/examples/2d_focusing_phase_conjugated_light_through_disorder/regular_focusing.gif.REMOVED.git-id b/examples/2d_focusing_inside_disorder_with_phase_conjugation/regular_focusing.gif.REMOVED.git-id similarity index 100% rename from examples/2d_focusing_phase_conjugated_light_through_disorder/regular_focusing.gif.REMOVED.git-id rename to examples/2d_focusing_inside_disorder_with_phase_conjugation/regular_focusing.gif.REMOVED.git-id diff --git a/examples/2d_open_channel_through_disorder/README.md b/examples/2d_open_channel_through_disorder/README.md index 1454350..220fce5 100644 --- a/examples/2d_open_channel_through_disorder/README.md +++ b/examples/2d_open_channel_through_disorder/README.md @@ -62,8 +62,9 @@ input.side = "low" output.side = "high" # put PML along z-direction -pml = get_optimal_PML(syst.wavelength/syst.dx) -pml.npixels = 15 +pml = mesti_optimal_pml_params(syst.wavelength/syst.dx) +pml_npixels = 15 +pml.npixels = pml_npixels syst.zPML = [pml] # transmission matrix: input from the low side, output to the high side @@ -89,7 +90,7 @@ Factorizing ... elapsed time: 122.678 secs ```julia # The most-open channels is the singular vector of the transmission matrix with # the largest singular value. -(_, sigma_max, v_max), _, _, _, _ = svds(t, nsv=1) +(_, sigma_max, v_open), _, _, _, _ = svds(t, nsv=1) N_prop_low = channels.low.N_prop # number of propagating channels on the low side ind_normal = Int(round((N_prop_low+1)/2)) # index of the normal-incident plane-wave @@ -114,13 +115,15 @@ println(" T_avg = ", @sprintf("%.2f", T_avg), "\n T_PW = ", @sprintf("%.2f" input = wavefront() input.v_low = zeros(ComplexF64, N_prop_low, 2) input.v_low[ind_normal, 1] = 1 -input.v_low[:, 2] = v_max +input.v_low[:, 2] = v_open # we will also get the field profile in the free spaces on the two sides, for # plotting purpose. opts = Opts() -opts.nz_low = round((L_tot-L)/2/dx) -opts.nz_high = opts.nz_low +nz_low = round(Int,(L_tot-L)/2/dx) +nz_high = nz_low +opts.nz_low = nz_low +opts.nz_high = nz_high # for field-profile computations Ex, _, _ = mesti2s(syst, input, opts) @@ -141,6 +144,81 @@ Solving ... elapsed time: 7.809 secs Total elapsed time: 147.811 secs ``` +# Define the matrix B by users and compare the field results + +```julia +# do the same field-profile computations through defining matrix B by users and using mesti() +syst = Syst() +ny_Ex = size(epsilon_xx,1) # total number of pixel along y for Ex grids + +# in mesti(), syst.epsilon_low and syst.epsilon_high are not needed +# in mesti(), we provide the whole epsilon_xx including the scattering region, source/detection region, and PML region +syst.epsilon_xx = cat(epsilon_low*ones(ny_Ex,pml_npixels+1), epsilon_xx, epsilon_high*ones(ny_Ex,pml_npixels+1), dims=2) +# in previous mesti2s() calculation, +# syst.epsilon_xx = epsilon_xx +# syst.epsilon_low = epsilon_low +# syst.epsilon_high = epsilon_high +syst.length_unit = "lambda_0" +syst.wavelength = 1 +syst.dx = dx +syst.yBC = yBC + +# put PML along z-direction +pml = mesti_optimal_pml_params(syst.wavelength/syst.dx) +pml.npixels = pml_npixels +pml.direction = "z" # put +syst.PML = [pml] +# in previous mesti2s() calculation, +# syst.zPML = [pml] (and do not need to specify pml.direction = "z") + +Bx = Source_struct() +u_prop_low_Ex = channels.u_x_m(channels.low.kydx_prop) # the transverse profiles, u_Ex, for the propagating ones + +# If it is a single propagating channel, a line source of -2i*sqrt(nu)*u_Ex(m) at l=0 will generate an z-flux-normalized incident field of exp(i*kzdx*|l|)*u_Ex(l)/sqrt(nu), where nu = sin(kzdx) +# here, we want input wavefront fields on the low side, which is a superposition of propagating channels +# these input wavefront fields on the low side are specified by v_low, and we take superpositions of the propagating channels using the v coefficients, with the sqrt(nu)*exp(-i*kzdx*0.5) prefactors included. +# we will multiply the -2i prefactor at the end. +# the source B_Ex would generate the input wavefront fields +B_Ex_low = u_prop_low_Ex*(channels.low.sqrt_nu_prop.*exp.((-1im*0.5)*channels.low.kzdx_prop).*v_low) # 0.5 pixel backpropagation indicates that the source is half a pixel away from z = 0 +Bx.data = [B_Ex_low] # the value of the source + +# the position of the source specify by a rectangle. +# [m1, l1, w, h] specifies the location and the size of the +# rectangle. Here, (m1, l1) is the index of the (y,z) coordinate of +# the smaller-y, smaller-z corner of the rectangle, at the location +# of Ex(m1, l1); (w, h) is the width and height of the rectangle +# line source is put on the lower source region, where is one pixel outside PML region. +Bx.pos = [[1,pml_npixels+1,ny_Ex,1]] + +opts = Opts() +# the -2i prefactor would be multiplied by. +opts.prefactor = -2im + +Ex_prime, _ = mesti(syst, [Bx], opts) +``` +```text:Output +===System size=== +ny_Ex = 5400; nz_Ex = 1381 for Ex(y,z) +UPML on -z +z sides; ; yBC = periodic; zBC = PEC +Building B,C... elapsed time: 0.001 secs +Building A ... elapsed time: 6.449 secs +< Method: factorize_and_solve using MUMPS in single precision with AMD ordering > +Analyzing ... elapsed time: 3.388 secs +Factorizing ... elapsed time: 102.293 secs +Solving ... elapsed time: 5.059 secs + Total elapsed time: 112.826 secs +``` + +```julia +# Excluding the extra padding and PML region, compare these field profiles +println("Maximum absolute value of field difference between constructing the source matrix B through mesti2s() and constructing by users = ", +maximum(abs.(Ex[:,nz_low+1:end-nz_high-1,:] - Ex_prime[:,pml_npixels+1+1:end-pml_npixels-1-1,:]))) +``` +```text:Output +Maximum absolute value of field difference between constructing the source matrix B through mesti2s() and constructing by users = 0.0 +``` + + # Animate the field profiles ```julia diff --git a/examples/2d_open_channel_through_disorder/open_channel_through_disorder.ipynb b/examples/2d_open_channel_through_disorder/open_channel_through_disorder.ipynb index d94284e..0f239d2 100644 --- a/examples/2d_open_channel_through_disorder/open_channel_through_disorder.ipynb +++ b/examples/2d_open_channel_through_disorder/open_channel_through_disorder.ipynb @@ -106,8 +106,9 @@ "output.side = \"high\"\n", "\n", "# put PML along z-direction\n", - "pml = get_optimal_PML(syst.wavelength/syst.dx)\n", - "pml.npixels = 15\n", + "pml = mesti_optimal_pml_params(syst.wavelength/syst.dx)\n", + "pml_npixels = 15\n", + "pml.npixels = pml_npixels\n", "syst.zPML = [pml]\n", "\n", "# transmission matrix: input from the low side, output to the high side\n", @@ -131,7 +132,7 @@ "source": [ "# The most-open channels is the singular vector of the transmission matrix with \n", "# the largest singular value.\n", - "(_, sigma_max, v_max), _, _, _, _ = svds(t, nsv=1)\n", + "(_, sigma_max, v_open), _, _, _, _ = svds(t, nsv=1)\n", "\n", "N_prop_low = channels.low.N_prop # number of propagating channels on the low side\n", "ind_normal = Int(round((N_prop_low+1)/2)) # index of the normal-incident plane-wave\n", @@ -157,18 +158,101 @@ "input = wavefront()\n", "input.v_low = zeros(ComplexF64, N_prop_low, 2)\n", "input.v_low[ind_normal, 1] = 1\n", - "input.v_low[:, 2] = v_max\n", + "input.v_low[:, 2] = v_open\n", "\n", "# we will also get the field profile in the free spaces on the two sides, for\n", "# plotting purpose.\n", "opts = Opts()\n", - "opts.nz_low = round((L_tot-L)/2/dx)\n", - "opts.nz_high = opts.nz_low\n", + "nz_low = round(Int,(L_tot-L)/2/dx)\n", + "nz_high = nz_low\n", + "opts.nz_low = nz_low\n", + "opts.nz_high = nz_high\n", "\n", "# for field-profile computations\n", "Ex, _, _ = mesti2s(syst, input, opts);" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Define the matrix B by users and compare the field results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# do the same field-profile computations through defining matrix B by users and using mesti()\n", + "syst = Syst()\n", + "ny_Ex = size(epsilon_xx,1) # total number of pixel along y for Ex grids\n", + "\n", + "# in mesti(), syst .epsilon_low and syst.epsilon_high are not needed\n", + "# in mesti(), we provide the whole epsilon_xx including the scattering region, source/detection region, and PML region\n", + "syst.epsilon_xx = cat(epsilon_low*ones(ny_Ex,pml_npixels+1), epsilon_xx, epsilon_high*ones(ny_Ex,pml_npixels+1), dims=2)\n", + "# in previous mesti2s() calculation, \n", + "# syst.epsilon_xx = epsilon_xx\n", + "# syst.epsilon_low = epsilon_low \n", + "# syst.epsilon_high = epsilon_high\n", + "syst.length_unit = \"lambda_0\"\n", + "syst.wavelength = 1\n", + "syst.dx = dx\n", + "syst.yBC = yBC \n", + "\n", + "# put PML along z-direction\n", + "pml = mesti_optimal_pml_params(syst.wavelength/syst.dx)\n", + "pml.npixels = pml_npixels\n", + "pml.direction = \"z\" # put\n", + "syst.PML = [pml]\n", + "# in previous mesti2s() calculation, \n", + "# syst.zPML = [pml] (and do not need to specify pml.direction = \"z\")\n", + "\n", + "Bx = Source_struct()\n", + "u_prop_low_Ex = channels.u_x_m(channels.low.kydx_prop) # the transverse profiles, u_Ex, for the propagating ones\n", + "\n", + "# If it is a single propagating channel, a line source of -2i*sqrt(nu)*u_Ex(m) at l=0 will generate an z-flux-normalized incident field of exp(i*kzdx*|l|)*u_Ex(l)/sqrt(nu), where nu = sin(kzdx)\n", + "# here, we want input wavefront fields on the low side, which is a superposition of propagating channels\n", + "# these input wavefront fields on the low side are specified by v_low, and we take superpositions of the propagating channels using the v coefficients, with the sqrt(nu)*exp(-i*kzdx*0.5) prefactors included.\n", + "# we will multiply the -2i prefactor at the end.\n", + "# the source B_Ex would generate the input wavefront fields\n", + "B_Ex_low = u_prop_low_Ex*(channels.low.sqrt_nu_prop.*exp.((-1im*0.5)*channels.low.kzdx_prop).*v_low) # 0.5 pixel backpropagation indicates that the source is half a pixel away from z = 0\n", + "Bx.data = [B_Ex_low] # the value of the source\n", + "\n", + "# the position of the source specify by a rectangle.\n", + "# [m1, l1, w, h] specifies the location and the size of the\n", + "# rectangle. Here, (m1, l1) is the index of the (y,z) coordinate of\n", + "# the smaller-y, smaller-z corner of the rectangle, at the location\n", + "# of Ex(m1, l1); (w, h) is the width and height of the rectangle\n", + "# line source is put on the lower source region, where is one pixel outside PML region.\n", + "Bx.pos = [[1,pml_npixels+1,ny_Ex,1]] \n", + "\n", + "opts = Opts()\n", + "# the -2i prefactor would be multiplied by.\n", + "opts.prefactor = -2im\n", + "\n", + "Ex_prime, _ = mesti(syst, [Bx], opts)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Excluding the extra padding and PML region, compare these field profiles\n", + "println(\"Maximum absolute value of field difference between constructing the source matrix B through mesti2s() and constructing by users = \", \n", + " maximum(abs.(Ex[:,nz_low+1:end-nz_high-1,:] - Ex_prime[:,pml_npixels+1+1:end-pml_npixels-1-1,:])))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Animate the field profiles" + ] + }, { "cell_type": "code", "execution_count": null, @@ -215,15 +299,15 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.9.3", + "display_name": "Julia 1.8.0", "language": "julia", - "name": "julia-1.9" + "name": "julia-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.9.3" + "version": "1.8.0" } }, "nbformat": 4, diff --git a/examples/2d_open_channel_through_disorder/open_channel_through_disorder.jl b/examples/2d_open_channel_through_disorder/open_channel_through_disorder.jl index e57cdab..147aebb 100644 --- a/examples/2d_open_channel_through_disorder/open_channel_through_disorder.jl +++ b/examples/2d_open_channel_through_disorder/open_channel_through_disorder.jl @@ -13,7 +13,6 @@ using MESTI, GeometryPrimitives, Arpack, Printf include("build_epsilon_disorder.jl") ## System parameters - # dimensions of the system, in units of the wavelength lambda_0 dx = 1/15 # discretization grid size W = 360 # width of the scattering region @@ -22,8 +21,8 @@ L_tot = 150 # full length of the system for plotting r_min = 0.2 # minimal radius of the cylindrical scatterers r_max = 0.4 # maximal radius of the cylindrical scatterers min_sep = 0.05 # minimal separation between cylinders -number_density = 1.3 # number density, in units of 1/lambda_0^2 -rng_seed = 0 # random number generator seed +number_density = 1.3 # number density, in units of 1/lambda_0^2 +rng_seed = 0 # random number generator seed # relative permittivity, unitless epsilon_scat = 1.2^2 # cylindrical scatterers @@ -41,7 +40,7 @@ build_epsilon_disorder(W, L, r_min, r_max, min_sep, number_density, rng_seed, dx, epsilon_scat, epsilon_bg, build_TM) -## Compute the transmission matrix for TM waves +## Compute the transmission matrix syst = Syst() syst.epsilon_xx = epsilon_xx syst.epsilon_low = epsilon_low @@ -60,8 +59,9 @@ input.side = "low" output.side = "high" # put PML along z-direction -pml = get_optimal_PML(syst.wavelength/syst.dx) -pml.npixels = 15 +pml = mesti_optimal_pml_params(syst.wavelength/syst.dx) +pml_npixels = 15 +pml.npixels = pml_npixels syst.zPML = [pml] # transmission matrix: input from the low side, output to the high side @@ -70,7 +70,7 @@ t, channels, _ = mesti2s(syst, input, output) ## Compare an open channel and a plane-wave input # The most-open channels is the singular vector of the transmission matrix with # the largest singular value. -(_, sigma_max, v_max), _, _, _, _ = svds(t, nsv=1) +(_, sigma_max, v_open), _, _, _, _ = svds(t, nsv=1) N_prop_low = channels.low.N_prop # number of propagating channels on the low side ind_normal = Int(round((N_prop_low+1)/2)) # index of the normal-incident plane-wave @@ -85,20 +85,77 @@ println(" T_avg = ", @sprintf("%.2f", T_avg), "\n T_PW = ", @sprintf("%.2f" # specify two input incident wavefronts: # (1) normal-incident plane-wave # (2) open channel +v_low = zeros(ComplexF64, N_prop_low, 2) +v_low[ind_normal, 1] = 1 +v_low[:, 2] = v_open input = wavefront() -input.v_low = zeros(ComplexF64, N_prop_low, 2) -input.v_low[ind_normal, 1] = 1 -input.v_low[:, 2] = v_max +input.v_low = v_low # we will also get the field profile in the free spaces on the two sides, for # plotting purpose. opts = Opts() -opts.nz_low = round((L_tot-L)/2/dx) -opts.nz_high = opts.nz_low +nz_low = round(Int,(L_tot-L)/2/dx) +nz_high = nz_low +opts.nz_low = nz_low +opts.nz_high = nz_high -# for field-profile computations +# compute the field-profiles through mesti() Ex, _, _ = mesti2s(syst, input, opts) +## Define the matrix B by users and compare the field results +# do the same field-profile computations through defining matrix B by users and using mesti() +syst = Syst() +ny_Ex = size(epsilon_xx,1) # total number of pixel along y for Ex grids + +# in mesti(), syst.epsilon_low and syst.epsilon_high are not needed +# in mesti(), we provide the whole epsilon_xx including the scattering region, source/detection region, and PML region +syst.epsilon_xx = cat(epsilon_low*ones(ny_Ex,pml_npixels+1), epsilon_xx, epsilon_high*ones(ny_Ex,pml_npixels+1), dims=2) +# in previous mesti2s() calculation, +# syst.epsilon_xx = epsilon_xx +# syst.epsilon_low = epsilon_low +# syst.epsilon_high = epsilon_high +syst.length_unit = "lambda_0" +syst.wavelength = 1 +syst.dx = dx +syst.yBC = yBC + +# put PML along z-direction +pml = mesti_optimal_pml_params(syst.wavelength/syst.dx) +pml.npixels = pml_npixels +pml.direction = "z" # put +syst.PML = [pml] +# in previous mesti2s() calculation, +# syst.zPML = [pml] (and do not need to specify pml.direction = "z") + +Bx = Source_struct() +u_prop_low_Ex = channels.u_x_m(channels.low.kydx_prop) # the transverse profiles, u_Ex, for the propagating ones + +# If it is a single propagating channel, a line source of -2i*sqrt(nu)*u_Ex(m) at l=0 will generate an z-flux-normalized incident field of exp(i*kzdx*|l|)*u_Ex(l)/sqrt(nu), where nu = sin(kzdx) +# here, we want input wavefront fields on the low side, which is a superposition of propagating channels +# these input wavefront fields on the low side are specified by v_low, and we take superpositions of the propagating channels using the v coefficients, with the sqrt(nu)*exp(-i*kzdx*0.5) prefactors included. +# we will multiply the -2i prefactor at the end. +# the source B_Ex would generate the input wavefront fields +B_Ex_low = u_prop_low_Ex*(channels.low.sqrt_nu_prop.*exp.((-1im*0.5)*channels.low.kzdx_prop).*v_low) # 0.5 pixel backpropagation indicates that the source is half a pixel away from z = 0 +Bx.data = [B_Ex_low] # the value of the source + +# the position of the source specify by a rectangle. +# [m1, l1, w, h] specifies the location and the size of the +# rectangle. Here, (m1, l1) is the index of the (y,z) coordinate of +# the smaller-y, smaller-z corner of the rectangle, at the location +# of Ex(m1, l1); (w, h) is the width and height of the rectangle +# line source is put on the lower source region, where is one pixel outside PML region. +Bx.pos = [[1,pml_npixels+1,ny_Ex,1]] + +opts = Opts() +# the -2i prefactor would be multiplied by. +opts.prefactor = -2im + +Ex_prime, _ = mesti(syst, [Bx], opts) + +# Excluding the extra padding and PML region, compare these field profiles +println("Maximum absolute value of field difference between constructing the source matrix B through mesti2s() and constructing by users = ", + maximum(abs.(Ex[:,nz_low+1:end-nz_high-1,:] - Ex_prime[:,pml_npixels+1+1:end-pml_npixels-1-1,:]))) + ## Animate the field profiles using Plots # normalize the field amplitude with respect to the plane-wave-input profile diff --git a/examples/2d_reflection_matrix_Gaussian_beams/README.md b/examples/2d_reflection_matrix_Gaussian_beams/README.md index 49dda2a..05f5bc8 100644 --- a/examples/2d_reflection_matrix_Gaussian_beams/README.md +++ b/examples/2d_reflection_matrix_Gaussian_beams/README.md @@ -51,7 +51,7 @@ heatmap(z, y, epsilon_xx, We consider inputs being Gaussian beams focused at (*y*f, *z*f). In this example, we fix the focal depth at *z*f = *z*0 (i.e., the depth of the scatterer), and scan the transverse coordinate *y*f of the focus. -Perfect Gaussian beams can be generated with the total-field/scattered-field (TF/SF) method. But since the cross section of the beam decays exponentially in *y*, we can generate Gaussian beams to a high accuracy simply by placing line sources at a cross section on the left, which is what we do here. We place the line sources at *z* = *z*source, just in front of the PML. +Perfect Gaussian beams can be generated with the total-field/scattered-field (TF/SF) method. But since the cross section of the beam decays exponentially in *y*, we can generate Gaussian beams to a high accuracy simply by placing line sources at a cross section on the low side, which is what we do here. We place the line sources at *z* = *z*source, just in front of the PML. To determine the required line sources, we (1) take the field profile of the desired incident Gaussian beam at the focal plane, *E*in(*y*,*z*f) = *E*0exp(-(*y* - *y*f)2/*w*2), (2) project it onto the propagating channels (i.e., ignoring evanescent contributions) of free space, (3) back propagate it to the source plane to determine *E*in(*y*,*z*source) in the propagating-channel basis, and (4) determine the line source necessary to generate such *E*in(*y*,*z*source). diff --git a/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.ipynb b/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.ipynb index 1610c4c..0c486aa 100644 --- a/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.ipynb +++ b/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.ipynb @@ -98,7 +98,7 @@ "source": [ "We consider inputs being Gaussian beams focused at (*y*f, *z*f). In this example, we fix the focal depth at *z*f = *z*0 (i.e., the depth of the scatterer), and scan the transverse coordinate *y*f of the focus.\n", "\n", - "Perfect Gaussian beams can be generated with the total-field/scattered-field (TF/SF) method. But since the cross section of the beam decays exponentially in *y*, we can generate Gaussian beams to a high accuracy simply by placing line sources at a cross section on the left, which is what we do here. We place the line sources at *z* = *z*source, just in front of the PML.\n", + "Perfect Gaussian beams can be generated with the total-field/scattered-field (TF/SF) method. But since the cross section of the beam decays exponentially in *y*, we can generate Gaussian beams to a high accuracy simply by placing line sources at a cross section on the low side, which is what we do here. We place the line sources at *z* = *z*source, just in front of the PML.\n", "\n", "To determine the required line sources, we (1) take the field profile of the desired incident Gaussian beam at the focal plane, *E*in(*y*,*z*f) = *E*0exp(-(*y* - *y*f)2/*w*2), (2) project it onto the propagating channels (i.e., ignoring evanescent contributions) of free space, (3) back propagate it to the source plane to determine *E*in(*y*,*z*source) in the propagating-channel basis, and (4) determine the line source necessary to generate such *E*in(*y*,*z*source)." ] diff --git a/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.jl b/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.jl index 77b3500..ac312f8 100644 --- a/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.jl +++ b/examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.jl @@ -56,7 +56,7 @@ heatmap(z, y, epsilon_xx, # total-field/scattered-field (TF/SF) method. But since the cross section # of the beam decays exponentially in y, we can generate Gaussian beams to # a high accuracy simply by placing line sources at a cross section on the -# left, which is what we do here. We place the line sources at z = +# low side, which is what we do here. We place the line sources at z = # z_source, just in front of the PML. # # To determine the required line sources, we (1) take the field profile diff --git a/src/MESTI.jl b/src/MESTI.jl index 82ff118..0f92973 100644 --- a/src/MESTI.jl +++ b/src/MESTI.jl @@ -18,12 +18,12 @@ include("mumps3_convenience_wrappers.jl") include("mumps3_icntl_alibis.jl") include("mumps3_printing.jl") -include("get_optimal_PML.jl") -include("build_transverse_function_1d.jl") +include("mesti_optimal_pml_params.jl") +include("mesti_build_transverse_function_1d.jl") include("mesti_subpixel_smoothing.jl") include("mesti_build_fdfd_matrix.jl") include("mesti_matrix_solver.jl") -include("setup_longitudinal.jl") +include("mesti_setup_longitudinal.jl") include("mesti_main.jl") include("mesti_build_channels.jl") include("mesti2s.jl") diff --git a/src/mesti2s.jl b/src/mesti2s.jl index 50185af..f64e60c 100644 --- a/src/mesti2s.jl +++ b/src/mesti2s.jl @@ -529,7 +529,7 @@ end opts.nthreads_OMP (positive integer scalar; optional): Number of OpenMP threads used in MUMPS; overwrites the OMP_NUM_THREADS environment variable. - opts.parallel_dependency_graph (logical scalar; optional): + opts.parallel_dependency_graph (logical scalar; optional, defaults to false): If MUMPS is multithread, whether to use parallel dependency graph in MUMPS. This typically improve the time performance, but marginally increase the memory usage. @@ -610,7 +610,7 @@ end numbers given by vectors channels.low.kzdx_prop and channels.high.kzdx_prop. The phases of the elements of the scattering matrix depend on the reference plane. For channels on the low, the reference plane is at - z = 0. For channels on the high, the reference plane is at z = H. + z = 0. For channels on the high, the reference plane is at z = L. When a subset of the propagating channels are requested, in 3D ''input = channel_index()'' and ''output = channel_index()'' with in.ind_low_s, in.ind_low_p, in.ind_high_s, in.ind_high_p, @@ -716,7 +716,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront (ny_Ex, nz_Ex) = size(syst.epsilon_xx) end - dn = 0.5; # the source/detection is half a pixel away from z=0 and z=d + dn = 0.5 # the source/detection is half a pixel away from z = 0 and z = L if ~use_2D_TM # Check boundary condition in x @@ -918,7 +918,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront throw(ArgumentError("opts.nz_low must be a non-negative integer scalar, if given.")) end if ~isdefined(opts, :nz_high) || isa(opts.nz_high, Nothing) - opts.nz_high = 0; + opts.nz_high = 0 elseif ~(opts.nz_high >= 0 ) throw(ArgumentError("opts.nz_high must be a non-negative integer scalar, if given.")) end @@ -973,9 +973,9 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront end end - method_specified = true; + method_specified = true if ~isdefined(opts, :method) || isa(opts.method, Nothing) - method_specified = false; + method_specified = false # By default, if the argument ''output'' is not given or if # opts.iterative_refinement = true, then "factorize_and_solve" is used. # Otherwise, then "APF" is used when opts.solver = "MUMPS", and @@ -1101,7 +1101,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront @printf("\n") end - t1 = time(); timing_init = t1-t0; # Initialization time + t1 = time(); timing_init = t1-t0 # Initialization time ## Part 2.1: Parse the argument ''input'' if opts.verbal; @printf("Building B,C... "); end @@ -1404,7 +1404,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront if isdefined(output, :polarization) && ~(output.polarization in ["s", "p", "both"]) throw(ArgumentError("output.polarization, if given, must be \"s\", \"p\", or \"both\".")) elseif ~isdefined(output, :polarization) - # Pick the polarization to use when uers specify the side; + # Pick the polarization to use when uers specify the side output.polarization = "both" end end @@ -1642,12 +1642,12 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront # No need to build C if we symmetrize K if opts.return_field_profile || use_transpose_B - build_C = false; + build_C = false else - build_C = true; + build_C = true end - # We only need the transverse functions for Ex, Ey, and derivative of the transverse functions for Ez; see Ref XXX ...... + # We only need the transverse functions for Ex, Ey, and derivative of the transverse functions for Ez; see Ref which would be published later. if ~use_2D_TM @cast u_prop_low_Ex[(n,m),N_prop_low] := channels.u_x_m(channels.low.kydx_prop)[m,N_prop_low] * channels.u_x_n(channels.low.kxdx_prop)[n,N_prop_low] @cast u_prop_low_Ey[(n,m),N_prop_low] := channels.u_y_m(channels.low.kydx_prop)[m,N_prop_low] * channels.u_y_n(channels.low.kxdx_prop)[n,N_prop_low] @@ -1667,12 +1667,12 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront end end else - @cast u_prop_low_Ex[m,N_prop_low] := channels.u_x_m(channels.low.kydx_prop)[m,N_prop_low] + u_prop_low_Ex = channels.u_x_m(channels.low.kydx_prop) if two_sided if (syst.epsilon_high == syst.epsilon_low) u_prop_high_Ex = u_prop_low_Ex else - @cast u_prop_high_Ex[m,N_prop_high] := channels.u_x_m(channels.high.kydx_prop)[m,N_prop_high] + u_prop_high_Ex = channels.u_x_m(channels.high.kydx_prop) end end end @@ -1682,7 +1682,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront # B_Ej_low, C_Ej_low, B_Ej_high, C_Ej_high are all dense matrices with size(..., 1) = nx_Ej*ny_Ej. They are inputs/outputs on a slice surface in xy-plane. # Inputs/outputs are placed at one pixel outside syst.epsilon_xx/epsilon_yy/epsilon_zz. # A line source including back propagation of B_Ex_low = -2im*sqrt(nu)*alpha_x(a,b,sigma,+)*u_x(n,m,a,b)*exp(1im*kz(a,b)*dx*(-1/2))+2/sqrt(nu)*alpha_z(a,b,sigma,+)*(u_z(n+1,m,a,b)-u_z(n,m,a,b))*cos(kz(a,b)*dx/2)*exp(1im*kz(a,b)*dx*(-1/2)) and B_Ey_low = -2im*sqrt(nu)*alpha_y(a,b,sigma,+)*u_y(n,m,a,b)*exp(1im*kz(a,b)*dx*(-1/2))+2/sqrt(nu)*alpha_z(a,b,sigma,+)*(u_z(n,m+1,a,b)-u_z(n,m,a,b))*cos(kz(a,b)*dx/2)*exp(1im*kz(a,b)*dx*(-1/2)) at l=0 will generate an z-flux-normalized incident field of Ex, Ey, and Ez. Ex = alpha_x(a,b,sigma,+/-)/sqrt(nu)*alpha_x(a,b,sigma,+/-)*u_x(n,m,a,b)*exp(1im*kz(a,b)*dx*(|l|-1/2)), where + for l>= 0 and - for l< 0. Ey = alpha_y(a,b,sigma,+/-)/sqrt(nu)*alpha_y(a,b,sigma,+/-)*u_x(n,m,a,b)*exp(1im*kz(a,b)*dx*(|l|-1/2)), where + for l>= 0 and - for l< 0. Ez = alpha_y(a,b,sigma,+/-)/sqrt(nu)*alpha_y(a,b,sigma,+/-)*u_x(n,m,a,b)*exp(1im*kz(a,b)*dx*|l|), where + for l>= 0 and - for l< 0. Note nu = sin(kzdx). - # We will multiple the -2i prefactor at the end. + # We will multiply the -2i prefactor at the end. # The flux-normalized output projection is \sum\limits_{j = x,y}{sqrt(nu)*conj(alpha_j(c,d,sigma,+/-)*u_x(n,m,c,d)*exp((+/-)1im*kz(c,d)*dx*(1/2)))+1/sqrt(nu)*conj(1im*alpha_z(c,d,sigma,+/-)*(u_z(n+delta_ix,m+delta_iy,c,d)-u_z(n,m,c,d))*cos(kz(c,d)*dx/2)))*exp(1im*kz(c,d)*dx*(-1/2)}; it will be transposed in mesti(). We also need to shift the phase by back propagating. # Therefore, we want C_Ex_low = sqrt(nu)*conj(alpha_x(c,d,sigma,-)*u_x(n,m,c,d)*exp(-1im*kz(c,d)*dx*(1/2)))+1/sqrt(nu)*conj(1im*alpha_z(c,d,sigma,-)*(u_z(n+1,m,c,d)-u_z(n,m,c,d))*cos(kz(c,d)*dx/2)))*exp(1im*kz(c,d)*dx*(-1/2) and C_Ey_low = sqrt(nu)*conj(alpha_y(c,d,sigma,-)*u_x(n,m,c,d)*exp(-1im*kz(c,d)*dx*(1/2)))+1/sqrt(nu)*conj(1im*alpha_z(c,d,sigma,-)*(u_z(n+1,m,c,d)-u_z(n,m,c,d))*cos(kz(c,d)*dx/2)))*exp(1im*kz(c,d)*dx*(-1/2). # Note that the complex conjugation only applies to u; the sqrt(nu) prefactor is not conjugated. (At real-valued frequency, nu is real-valued, so this doesn't matter. But at complex-valued frequency, nu is complex-valued, and we should not conjugate it. Note that the transverse basis is complete and orthonormal even when the frequency is complex, so the output projection doesn't need to be modified when the frequency is complex.) @@ -1917,7 +1917,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront C_high_s_Ex = nothing; C_high_s_Ey = nothing; C_high_p_Ex = nothing; C_high_p_Ey = nothing C_high_p_dEz_over_dx = nothing; C_high_p_dEz_over_dy = nothing else - u_prop_low_Ex = nothing; u_prop_high_Ex = nothing; + u_prop_low_Ex = nothing; u_prop_high_Ex = nothing end GC.gc() end @@ -1997,7 +1997,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront B_Ex = Source_struct() B_Ey = Source_struct() B_Ez = Source_struct() - B_Ez.isempty = true # The source only has x and y components; see Ref XXX ...... + B_Ez.isempty = true # The source only has x and y components; see Ref which would be published later. if ~two_sided B_Ex.pos = [[1, 1, l_low, nx_Ex, ny_Ex, 1]] B_Ey.pos = [[1, 1, l_low, nx_Ey, ny_Ey, 1]] @@ -2028,7 +2028,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront C_Ex = Source_struct() C_Ey = Source_struct() C_Ez = Source_struct() - C_Ez.isempty = true # The output projection only has x and y components; see Ref XXX ...... + C_Ez.isempty = true # The output projection only has x and y components; see Ref which would be published later. if ~two_sided C_Ex.pos = [[1, 1, l_low, nx_Ex, ny_Ex, 1]] C_Ey.pos = [[1, 1, l_low, nx_Ey, ny_Ey, 1]] @@ -2234,11 +2234,11 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront prefactor = [prefactor; channels.high.sqrt_nu_prop[ind_out_high].*exp.((-1im*dn)*channels.high.kzdx_prop[ind_out_high])] end end - S = prefactor.*S; # use implicit expansion + S = prefactor.*S # use implicit expansion end # Subtract D = C*inv(A_0)*B - S_0 where A_0 is a reference system and S_0 is its scattering matrix - # The form of D is summation phase shift by dn over x, y, and z components; see Ref XXX ...... + # The form of D is summation phase shift by dn over x, y, and z components; see Ref which would be published later. # When user-specified input and output wavefronts are used, we have D_low = (v_out_low')*diag(D)*v_in_low if ~use_2D_TM phase_factor = repeat(exp.((-1im*2*dn)*channels.low.kzdx_prop), 2) @@ -2390,7 +2390,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront # u where the a-th column is the a-th channels; it includes all the propagating and evanescent channels. # u_low = [[u_Ex.*reshape(alpha_x_all_low_s,1,:); u_Ey.*reshape(alpha_y_all_low_s,1,:)] - # [u_Ex.*reshape(alpha_x_all_low_p,1,:)+1im*u_dEz_over_dx*reshape(cos.(channels.low.kzdx_all/2).*alpha_z_all_low_p./sin.(channels.low.kzdx_all),1,:); + # [u_Ex.*reshape(alpha_x_all_low_p,1,:)+1im*u_dEz_over_dx*reshape(cos.(channels.low.kzdx_all/2).*alpha_z_all_low_p./sin.(channels.low.kzdx_all),1,:) # u_Ey.*reshape(alpha_y_all_low_p,1,:)+1im*u_dEz_over_dy*reshape(cos.(channels.low.kzdx_all/2).*alpha_z_all_low_p./sin.(channels.low.kzdx_all),1,:)]] # u_low = [[u_Ex.*reshape(alpha_x_all_low_s,1,:); u_Ey.*reshape(alpha_y_all_low_s,1,:); u_Ez.*reshape(alpha_z_all_low_s,1,:)] # [u_Ex.*reshape(alpha_x_all_low_p,1,:); u_Ey.*reshape(alpha_y_all_low_p,1,:); u_Ez.*reshape(alpha_z_all_low_p,1,:)]] @@ -2542,7 +2542,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront kz_z = reshape(channels.low.kzdx_all, ny_Ex, 1).*reshape(l, 1, :) # kz*z; ny_Ex-by-nz_low_extra matrix through implicit expansion exp_mikz = exp.(-1im*kz_z) # exp(-i*kz*z) kz_z_prop = reshape(channels.low.kzdx_prop, channels.low.N_prop, 1).*reshape(l, 1, :) # kz*z; channels.low.N_prop-by-nz_low_extra matrix through implicit expansion - exp_pikz_prop = exp.( 1im*kz_z_prop); # exp(+i*kz*z) + exp_pikz_prop = exp.( 1im*kz_z_prop) # exp(+i*kz*z) c_in = zeros(ComplexF64, ny_Ex, 1) c_in_prop = zeros(ComplexF64, channels.low.N_prop, 1) l_low = 1 # index for the inputs/outputs on the low surface @@ -2687,7 +2687,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront end for ii = 1:M_in_high # input from high c = u_prime*Ex[:, l_high, M_in_low+ii] # c is a ny_Ex-by-1 column vector of transverse mode coefficients - # c_in is the incident wavefront at n_R; note we need to back propagate dn pixel from z=d + # c_in is the incident wavefront at n_R; note we need to back propagate dn pixel from z = L c_in[:] .= 0 if use_ind_in c_in_prop[ind_in_high[ii]] = prefactor[ind_in_high[ii]] diff --git a/src/mesti_build_channels.jl b/src/mesti_build_channels.jl index b4b7085..2b57227 100644 --- a/src/mesti_build_channels.jl +++ b/src/mesti_build_channels.jl @@ -270,7 +270,7 @@ function mesti_build_channels(nx_Ex::Union{Int64,Nothing}, nx_Ey::Union{Int64,No yBC = "Bloch" end - # f = [f(1), ..., f(nx)].'; + # f = [f(1), ..., f(nx)].' # For periodic and Bloch periodic boundary, we order channels.kxdx_all (channels.kydx_all) such that it increases monotonically from negative to positive # For other boundary conditions, kx >= 0 (ky >= 0), and we order channels.kxdx_all (channels.kydx_all) such that it increases monotonically from smallest to largest if two_sided @@ -280,18 +280,18 @@ function mesti_build_channels(nx_Ex::Union{Int64,Nothing}, nx_Ey::Union{Int64,No end if ~use_2D_TM - (channels.u_x_n, channels.kxdx_all) = build_transverse_function_1d(nx_Ex, BC_x_x, n0, true) - (channels.u_x_m, _) = build_transverse_function_1d(ny_Ex, BC_x_y, n0) - (channels.u_y_n, _) = build_transverse_function_1d(nx_Ey, BC_y_x, n0) - (channels.u_y_m, channels.kydx_all) = build_transverse_function_1d(ny_Ey, BC_y_y, m0, true) + (channels.u_x_n, channels.kxdx_all) = mesti_build_transverse_function_1d(nx_Ex, BC_x_x, n0, true) + (channels.u_x_m, _) = mesti_build_transverse_function_1d(ny_Ex, BC_x_y, n0) + (channels.u_y_n, _) = mesti_build_transverse_function_1d(nx_Ey, BC_y_x, n0) + (channels.u_y_m, channels.kydx_all) = mesti_build_transverse_function_1d(ny_Ey, BC_y_y, m0, true) channels.u_z_n = channels.u_y_n # u_z_n and u_y_n are same transverse function channels.u_z_m = channels.u_x_m # u_z_m and u_x_m are same transverse function - channels.du_z_n = build_transverse_function_1d_derivative(nx_Ey, BC_z_x, n0, 1) - channels.du_z_m = build_transverse_function_1d_derivative(ny_Ex, BC_z_y, m0, 1) + channels.du_z_n = mesti_build_transverse_function_1d_derivative(nx_Ey, BC_z_x, n0, 1) + channels.du_z_m = mesti_build_transverse_function_1d_derivative(ny_Ex, BC_z_y, m0, 1) else - (channels.u_x_m, channels.kydx_all) = build_transverse_function_1d(ny_Ex, BC_x_y, m0) + (channels.u_x_m, channels.kydx_all) = mesti_build_transverse_function_1d(ny_Ex, BC_x_y, m0) channels.kxdx_all = nothing end @@ -312,7 +312,7 @@ function mesti_build_channels(nx_Ex::Union{Int64,Nothing}, nx_Ey::Union{Int64,No end # Properties for the homogeneous space on the low (kzdx, sqrt_nu_prop, number of propagating channels, etc; depends on epsilon_low/high) - side = setup_longitudinal(k0dx, epsilon_low, channels.kxdx_all, channels.kydx_all, kLambda_x, kLambda_y, ind_zero_kx, ind_zero_ky, use_continuous_dispersion) + side = mesti_setup_longitudinal(k0dx, epsilon_low, channels.kxdx_all, channels.kydx_all, kLambda_x, kLambda_y, ind_zero_kx, ind_zero_ky, use_continuous_dispersion) if two_sided channels.low = side @@ -320,7 +320,7 @@ function mesti_build_channels(nx_Ex::Union{Int64,Nothing}, nx_Ey::Union{Int64,No if epsilon_high == epsilon_low channels.high = side elseif ~isnan(epsilon_high) - channels.high = setup_longitudinal(k0dx, epsilon_high, channels.kxdx_all, channels.kydx_all, kLambda_x, kLambda_y, ind_zero_kx, ind_zero_ky, use_continuous_dispersion) + channels.high = mesti_setup_longitudinal(k0dx, epsilon_high, channels.kxdx_all, channels.kydx_all, kLambda_x, kLambda_y, ind_zero_kx, ind_zero_ky, use_continuous_dispersion) end else # Add the fields of "side" to "channels" @@ -357,12 +357,12 @@ function mesti_build_channels(syst::Syst) end if ~isdefined(syst, :epsilon_low); throw(ArgumentError("Input argument syst must have field \"epsilon_low\".")); end - epsilon_low = syst.epsilon_low; + epsilon_low = syst.epsilon_low if ~isdefined(syst, :wavelength); throw(ArgumentError("Input argument syst must have field \"wavelength\".")); end if ~isdefined(syst, :dx); throw(ArgumentError("Input argument syst must have field \"dx\".")); end if ~(syst.dx > 0); throw(ArgumentError("syst.dx must be a positive scalar.")); end - k0dx = (2*pi/syst.wavelength)*(syst.dx); + k0dx = (2*pi/syst.wavelength)*(syst.dx) # Check boundary condition in x if ~use_2D_TM diff --git a/src/mesti_build_fdfd_matrix.jl b/src/mesti_build_fdfd_matrix.jl index eb3dbc7..13f10b5 100644 --- a/src/mesti_build_fdfd_matrix.jl +++ b/src/mesti_build_fdfd_matrix.jl @@ -159,7 +159,7 @@ end """ function mesti_build_fdfd_matrix(epsilon_xx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Matrix{Int64},Matrix{Float64},Matrix{ComplexF64}}, epsilon_xy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_xz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_yx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_yy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_yz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_zx::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_zy::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, epsilon_zz::Union{Array{Int64,3},Array{Float64,3},Array{ComplexF64,3},Nothing}, k0dx::Union{Float64,ComplexF64}, xBC::Union{String,Int64,Float64,ComplexF64,Nothing}, yBC::Union{String,Int64,Float64,ComplexF64}, zBC::Union{String,Int64,Float64,ComplexF64}, xPML::Union{Vector{PML},Nothing} = [PML(0), PML(0)], yPML::Vector{PML} = [PML(0), PML(0)], zPML::Vector{PML} = [PML(0), PML(0)], use_UPML::Bool=true) # Make deepcopy of them to avoid mutating input argument - xPML = deepcopy(xPML); yPML = deepcopy(yPML); zPML = deepcopy(zPML); + xPML = deepcopy(xPML); yPML = deepcopy(yPML); zPML = deepcopy(zPML) # Take care of the 2D TM case if ndims(epsilon_xx) == 2 diff --git a/src/build_transverse_function_1d.jl b/src/mesti_build_transverse_function_1d.jl similarity index 95% rename from src/build_transverse_function_1d.jl rename to src/mesti_build_transverse_function_1d.jl index 4727e49..ed6bfa4 100644 --- a/src/build_transverse_function_1d.jl +++ b/src/mesti_build_transverse_function_1d.jl @@ -1,5 +1,5 @@ """ - BUILD_TRANSVERSE_FUNCTION_1D sets up 1D transverse function and wave number in the homogeneous space. + MESTI_BUILD_TRANSVERSE_FUNCTION_1D sets up 1D transverse function and wave number. === Input Arguments === nx (positive integer scalar; required): @@ -44,7 +44,7 @@ and are ordered from small to large. """ -function build_transverse_function_1d(nx::Int64, xBC::Union{String,Int64,Float64}, n0::Union{Int64,Float64}=0, offset::Bool=false) +function mesti_build_transverse_function_1d(nx::Int64, xBC::Union{String,Int64,Float64}, n0::Union{Int64,Float64}=0, offset::Bool=false) # Check input parameters if ~(nx>=0) throw(ArgumentError("Input argument nx must be a natural number.")) @@ -66,7 +66,7 @@ function build_transverse_function_1d(nx::Int64, xBC::Union{String,Int64,Float64 xBC = "Bloch" end - # f = [f(1), ..., f(nx)].'; + # f = [f(1), ..., f(nx)].' # For periodic and Bloch periodic boundary, we order kxdx_all such that it increases monotonically from negative to positive # For other boundary conditions, kx >= 0, and we order kxdx_all such that it increases monotonically from smallest to largest # Transverse modes in x (form a complete basis in x) @@ -121,7 +121,7 @@ function build_transverse_function_1d(nx::Int64, xBC::Union{String,Int64,Float64 return fun_u_1d, kxdx_all end """ - BUILD_HOMOGENEOUS_TRANSVERSE_FUNCTION_1D_DERIVATIVE sets up derivative of a transverse function. + MESTI_BUILD_TRANSVERSE_FUNCTION_1D_DERIVATIVE sets up derivative of a transverse function. === Input Arguments === nx (positive integer scalar; required): @@ -154,7 +154,7 @@ end when the input is a vector, it returns a matrix where each column is the respective derivative transverse profile. """ -function build_transverse_function_1d_derivative(nx::Int64, xBC::Union{String,Int64,Float64}, n0::Union{Int64,Float64}=0, changegrid::Union{Int64,Float64}=0) +function mesti_build_transverse_function_1d_derivative(nx::Int64, xBC::Union{String,Int64,Float64}, n0::Union{Int64,Float64}=0, changegrid::Union{Int64,Float64}=0) # Check input parameters if ~(nx>=0) error("Input argument nx must be a natural number.") diff --git a/src/mesti_main.jl b/src/mesti_main.jl index de163c0..6be3bea 100644 --- a/src/mesti_main.jl +++ b/src/mesti_main.jl @@ -571,7 +571,7 @@ end opts.nthreads_OMP (positive integer scalar; optional): Number of OpenMP threads used in MUMPS; overwrites the OMP_NUM_THREADS environment variable. - opts.parallel_dependency_graph (logical scalar; optional): + opts.parallel_dependency_graph (logical scalar; optional, defaults to false): If MUMPS is multithread, whether to use parallel dependency graph in MUMPS. This typically improve the time performance, but marginally increase the memory usage. @@ -1092,7 +1092,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC if opts.verbal # print basic system info if the calling function is not mesti2s() if stacktrace()[2].func == :mesti2s - called_from_mesti2s = true; + called_from_mesti2s = true @printf(" ... ") else called_from_mesti2s = false @@ -1249,7 +1249,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC n2 = n1 + pos[4] - 1 # last index in x m2 = m1 + pos[5] - 1 # last index in y l2 = l1 + pos[6] - 1 # last index in z - nxyz_data = pos[4]*pos[5]*pos[6]; # number of elements in this cuboid + nxyz_data = pos[4]*pos[5]*pos[6] # number of elements in this cuboid if n1 > nx_list[ii] if ii == 1; throw(ArgumentError("B[$ii].pos[$jj][1] = $(n1) exceeds nx_E$(component[ii]) = $(nx_Ex).")); end if ii == 2; throw(ArgumentError("B[$ii].pos[$jj][1] = $(n1) exceeds nx_E$(component[ii]) = $(nx_Ey).")); end @@ -1293,9 +1293,9 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC end if use_iv_pairs # convert to linear indices - n_list = repeat((n1:n2), 1, pos(5), pos(6)) - m_list = repeat(transpose(m1:m2), pos(4), 1, pos(6)) - l_list = repeat(reshape((l1:l2),1,1,:), pos(4), pos(5), 1) + n_list = reshape(repeat((n1:n2), 1, pos[5], pos[6]),:) + m_list = reshape(repeat(transpose(m1:m2), pos[4], 1, pos[6]),:) + l_list = reshape(repeat(reshape((l1:l2),1,1,:), pos[4], pos[5], 1),:) #ind = LinearIndices((nx_list[ii], ny_list[ii], nz_list[ii]))[CartesianIndex.(n_list, m_list, l_list)] ind = Base._sub2ind((nx_list[ii], ny_list[ii], nz_list[ii]), n_list, m_list, l_list) end @@ -1312,7 +1312,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC end if use_iv_pairs # Build index-value pairs: (ind_list, a_list, v_list) - N_ii = nxyz_data*M_ii; # number of nonzero elements in the jj-th part of matrix B_ii + N_ii = nxyz_data*M_ii # number of nonzero elements in the jj-th part of matrix B_ii ind_temp = N .+ (1:N_ii) ind_list[ind_temp] = repeat(ind, M_ii) # spatial index a_list[ind_temp] = reshape(repeat(M.+(1:M_ii), nxyz_data), N_ii) # input index @@ -1438,7 +1438,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC end if use_iv_pairs # Construct matrix C_ii from the complete set of index-value pairs - N_tot = 0; # total number of nonzero elements in C_ii + N_tot = 0 # total number of nonzero elements in C_ii for jj = 1:length(C_struct.data) N_tot = N_tot + length(C_struct.data[jj]) end @@ -1467,7 +1467,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC n2 = n1 + pos[4] - 1 # last index in x m2 = m1 + pos[5] - 1 # last index in y l2 = l1 + pos[6] - 1 # last index in z - nxyz_data = pos[4]*pos[5]*pos[6]; # number of elements in this cuboid + nxyz_data = pos[4]*pos[5]*pos[6] # number of elements in this cuboid if n1 > nx_list[ii] if ii == 1; throw(ArgumentError("C[$ii].pos[$jj][1] = $(n1) exceeds nx_E$(component[ii]) = $(nx_Ex).")); end if ii == 2; throw(ArgumentError("C[$ii].pos[$jj][1] = $(n1) exceeds nx_E$(component[ii]) = $(nx_Ey).")); end @@ -1511,9 +1511,9 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC end if use_iv_pairs # convert to linear indices - n_list = repeat((n1:n2), 1, pos[5], pos[6]) - m_list = repeat(transpose(m1:m2), pos[4], 1, pos[6]) - l_list = repeat(reshape((l1:l2),1,1,:), pos[4], pos[5], 1) + n_list = reshape(repeat((n1:n2), 1, pos[5], pos[6]),:) + m_list = reshape(repeat(transpose(m1:m2), pos[4], 1, pos[6]),:) + l_list = reshape(repeat(reshape((l1:l2),1,1,:), pos[4], pos[5], 1),:) #ind = LinearIndices((nx_list[ii], ny_list[ii], nz_list[ii]))[CartesianIndex.(n_list, m_list, l_list)] ind = Base._sub2ind((nx_list[ii], ny_list[ii], nz_list[ii]), n_list, m_list, l_list) end @@ -1530,7 +1530,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC end if use_iv_pairs # Build index-value pairs: (a_list, ind_list, v_list) - N_ii = nxyz_data*M_ii; # number of nonzero elements in the jj-th part of matrix C_ii + N_ii = nxyz_data*M_ii # number of nonzero elements in the jj-th part of matrix C_ii ind_temp = N .+ (1:N_ii) ind_list[ind_temp] = repeat(ind, M_ii) # spatial index a_list[ind_temp] = reshape(repeat(M.+(1:M_ii), nxyz_data), N_ii) # input index @@ -1771,7 +1771,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC if info.opts.exclude_PML_in_field_profiles # Exclude the PML pixels from the returned field profiles if ~use_2D_TM - n_start = 1; n_Ex_end = nx; n_Ey_end = nx; + n_start = 1; n_Ex_end = nx; n_Ey_end = nx Ex = reshape(S[1:nt_Ex, :], nx_Ex, ny_Ex, nz_Ex, M)[(info.xPML[1].npixels+1):(nx_Ex-info.xPML[2].npixels),(info.yPML[1].npixels+1):(ny_Ex-info.yPML[2].npixels),(info.zPML[1].npixels+1):(nz_Ex-info.zPML[2].npixels),:] Ey = reshape(S[nt_Ex+1:nt_Ex+nt_Ey, :], nx_Ey, ny_Ey, nz_Ey, M)[(info.xPML[1].npixels+1):(nx_Ey-info.xPML[2].npixels),(info.yPML[1].npixels+1):(ny_Ey-info.yPML[2].npixels),(info.zPML[1].npixels+1):(nz_Ey-info.zPML[2].npixels),:] Ez = reshape(S[nt_Ex+nt_Ey+1:nt_Ex+nt_Ey+nt_Ez, :], nx_Ez, ny_Ez, nz_Ez, M)[(info.xPML[1].npixels+1):(nx_Ez-info.xPML[2].npixels),(info.yPML[1].npixels+1):(ny_Ez-info.yPML[2].npixels),(info.zPML[1].npixels+1):(nz_Ez-info.zPML[2].npixels),:] diff --git a/src/mesti_matrix_solver.jl b/src/mesti_matrix_solver.jl index 815dc2b..1c5b66b 100644 --- a/src/mesti_matrix_solver.jl +++ b/src/mesti_matrix_solver.jl @@ -210,7 +210,7 @@ end opts.nthreads_OMP (positive integer scalar; optional): Number of OpenMP threads used in MUMPS; overwrites the OMP_NUM_THREADS environment variable. - opts.parallel_dependency_graph (logical scalar; optional): + opts.parallel_dependency_graph (logical scalar; optional, defaults to false): If MUMPS is multithread, whether to use parallel dependency graph in MUMPS. This typically improve the time performance, but marginally increase the memory usage. @@ -650,7 +650,7 @@ function mesti_matrix_solver!(matrices::Matrices, opts::Union{Opts,Nothing}=noth is_symmetric_K = false # even if A is symmetric, generally C won't equal transpose(B); we will not check whether C equals B.' or not; the user should set C = "transpose(B)" if C=B.' end if opts.clear_memory - matrices.A = nothing; matrices.B = nothing; matrices.C = nothing; D = nothing; + matrices.A = nothing; matrices.B = nothing; matrices.C = nothing; D = nothing GC.gc() end ind_schur = N .+ (1:M_tot) # indices for the Schur variables; must be a row vector @@ -743,7 +743,7 @@ function mesti_matrix_solver!(matrices::Matrices, opts::Union{Opts,Nothing}=noth S = Matrix(C_inv_U/L)*(P*(R\matrices.B)) # [[C*inv(U)]*inv(L)]*B end end - t2 = time(); info.timing_solve = t2-t1; + t2 = time(); info.timing_solve = t2-t1 if opts.verbal; @printf("elapsed time: %7.3f secs\n", info.timing_solve); end elseif opts.method == "factorize_and_solve" ## Compute S=C*inv(A)*B or X=inv(A)*B by factorizing A and solving for X column by column @@ -959,7 +959,6 @@ function MUMPS_analyze_and_factorize(A::Union{SparseMatrixCSC{Int64, Int64},Spar MPI.Initialized() ? nothing : MPI.Init() id = Mumps(A, sym=sym, par=par) # get the default parameters - set_icntl!(id,4,0;displaylevel=0); # Turn off diagnostic messages if opts.verbal_solver # Output to standard output stream, which is labeled by 6 in fortran # Note that the output behavior depends on the compiler used to compile MUMPS: diff --git a/src/get_optimal_PML.jl b/src/mesti_optimal_pml_params.jl similarity index 89% rename from src/get_optimal_PML.jl rename to src/mesti_optimal_pml_params.jl index f9cc772..0ce6f27 100644 --- a/src/get_optimal_PML.jl +++ b/src/mesti_optimal_pml_params.jl @@ -1,7 +1,7 @@ -# Export a function get_optimal_PML() -export get_optimal_PML +# Export a function mesti_optimal_pml_params() +export mesti_optimal_pml_params -function get_optimal_PML(wavelength_over_dx) +function mesti_optimal_pml_params(wavelength_over_dx) # Return an optimal PML parameters from double-log FOM given the mesh resolution wavelength_over_dx # wavelength_over_dx should be in the range of [5, 500]. diff --git a/src/setup_longitudinal.jl b/src/mesti_setup_longitudinal.jl similarity index 94% rename from src/setup_longitudinal.jl rename to src/mesti_setup_longitudinal.jl index c8ed8c9..e111d97 100644 --- a/src/setup_longitudinal.jl +++ b/src/mesti_setup_longitudinal.jl @@ -3,7 +3,7 @@ export Side mutable struct Side # A composite data type to store the items on a side - # See also: setup_longitudinal + # See also: mesti_setup_longitudinal N_prop::Integer kzdx_all::Vector{ComplexF64} ind_prop::Vector{Int64} @@ -16,7 +16,7 @@ mutable struct Side end """ - SETUP_LONGITUDINAL sets up a structure "Side" for one component in the homogeneous space. + MESTI_SETUP_LONGITUDINAL sets up a structure "Side" for one component in the homogeneous space. === Input Arguments === k0dx (numeric scalar, real or complex; required): Dimensionless frequency, k0*dx = (2*pi/vacuum_wavelength)*dx. @@ -83,7 +83,7 @@ end A permutation vector that switches one propagating channel with one having a complex-conjugated transverse profile. Now it only uses in 2D systems. """ -function setup_longitudinal(k0dx::Union{Float64,ComplexF64}, epsilon_bg::Union{Int64,Float64,ComplexF64}, kxdx_all::Union{StepRangeLen{Float64}, Vector{Float64}, Nothing}, kydx_all::Union{StepRangeLen{Float64}, Vector{Float64}}, kLambda_x::Union{Int64,Float64,ComplexF64,Nothing}=nothing, kLambda_y::Union{Int64,Float64,ComplexF64,Nothing}=nothing, ind_zero_kx::Union{Int64,Nothing}=nothing, ind_zero_ky::Union{Int64,Nothing}=nothing, use_continuous_dispersion::Bool=false) +function mesti_setup_longitudinal(k0dx::Union{Float64,ComplexF64}, epsilon_bg::Union{Int64,Float64,ComplexF64}, kxdx_all::Union{StepRangeLen{Float64}, Vector{Float64}, Nothing}, kydx_all::Union{StepRangeLen{Float64}, Vector{Float64}}, kLambda_x::Union{Int64,Float64,ComplexF64,Nothing}=nothing, kLambda_y::Union{Int64,Float64,ComplexF64,Nothing}=nothing, ind_zero_kx::Union{Int64,Nothing}=nothing, ind_zero_ky::Union{Int64,Nothing}=nothing, use_continuous_dispersion::Bool=false) if kxdx_all == nothing && kLambda_x == nothing && ind_zero_kx == nothing use_2D_TM = true diff --git a/src/mesti_subpixel_smoothing.jl b/src/mesti_subpixel_smoothing.jl index a55a5cb..79eca55 100644 --- a/src/mesti_subpixel_smoothing.jl +++ b/src/mesti_subpixel_smoothing.jl @@ -416,7 +416,7 @@ export mesti_subpixel_smoothing end if (~use_2D_TM && ~use_2D_TE) - throw(ArgumentError("In 2D case, use_2D_TM and/or use_2D_TE should be true")); + throw(ArgumentError("In 2D case, use_2D_TM and/or use_2D_TE should be true")) end # If the domain does not start from (0,0,0), translate the coordinates of the domain and obejcts. @@ -553,7 +553,8 @@ export mesti_subpixel_smoothing object = object_list[obj_ind] # looping over the range of the object in Eo grids. # Note that adding eps() can help avoid non-ideal situations where the last pixel might not be properly counted due to floating-point precision issues. - for kk = max((Int(div(bounds(object)[1][2]+delta_x/2-eps(), delta_x)) + 1),1):min((Int(div(bounds(object)[2][2]+delta_x/2+eps(), delta_x)) + 1),Nz), jj = max((Int(div(bounds(object)[1][1]+delta_x/2-eps(), delta_x)) + 1),1):min((Int(div(bounds(object)[2][1]+delta_x/2+eps(), delta_x)) + 1),Ny) iz = Ey_z_coord[kk]; iy = Ez_y_coord[jj] + for kk = max((Int(div(bounds(object)[1][2]+delta_x/2-eps(), delta_x)) + 1),1):min((Int(div(bounds(object)[2][2]+delta_x/2+eps(), delta_x)) + 1),Nz), jj = max((Int(div(bounds(object)[1][1]+delta_x/2-eps(), delta_x)) + 1),1):min((Int(div(bounds(object)[2][1]+delta_x/2+eps(), delta_x)) + 1),Ny) + iz = Ey_z_coord[kk]; iy = Ez_y_coord[jj] if without_sb if [iy,iz] ∈ object inv_epsilon_Eo_site[jj,kk,:,:] = ((object_epsilon_list[obj_ind]^(-1)) * Matrix(I,3,3)) @@ -672,7 +673,7 @@ export mesti_subpixel_smoothing end end - function pick_epsilon_3d(epsilon_Eo::Array{ComplexF64}, epsilon_Ex::Array{ComplexF64}, epsilon_Ey::Array{ComplexF64}, epsilon_Ez::Array{ComplexF64}, xBC::String, yBC::String, zBC::String) + function pick_epsilon_3d(epsilon_Eo::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, epsilon_Ex::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, epsilon_Ey::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, epsilon_Ez::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, xBC::String, yBC::String, zBC::String) # Pick the right site of components of epsilon based on Yee cell and boundary condition epsilon_xx = epsilon_Ex[:,:,:,1,1]; epsilon_xy = epsilon_Eo[:,:,:,1,2]; epsilon_xz = epsilon_Eo[:,:,:,1,3] @@ -687,7 +688,7 @@ export mesti_subpixel_smoothing #epsilon_xx = epsilon_xx[1:end-1,:,:]; epsilon_xy = epsilon_xy[:,:,:]; epsilon_xz = epsilon_xz[:,:,:] #epsilon_yx = epsilon_yx[:,:,:]; epsilon_yy = epsilon_yy[:,:,:]; epsilon_yz = epsilon_yz[:,:,:] #epsilon_zx = epsilon_zx[:,:,:]; epsilon_zy = epsilon_zy[:,:,:]; epsilon_zz = epsilon_zz[:,:,:] - epsilon_xx = epsilon_xx[1:end-1,:,:]; + epsilon_xx = epsilon_xx[1:end-1,:,:] elseif xBC == "PECPMC" epsilon_xx = epsilon_xx[1,end-1,:]; epsilon_xy = epsilon_xy[2:end,:,:]; epsilon_xz = epsilon_xz[2:end,:,:] epsilon_yx = epsilon_yx[2:end,:,:]; epsilon_yy = epsilon_yy[2:end,:,:]; epsilon_yz = epsilon_yz[2:end,:,:] @@ -706,7 +707,7 @@ export mesti_subpixel_smoothing #epsilon_xx = epsilon_xx[:,:,:]; epsilon_xy = epsilon_xy[:,:,:]; epsilon_xz = epsilon_xz[:,:,:] #epsilon_yx = epsilon_yx[:,:,:]; epsilon_yy = epsilon_yy[:,1:end-1,:]; epsilon_yz = epsilon_yz[:,:,:] #epsilon_zx = epsilon_zx[:,:,:]; epsilon_zy = epsilon_zy[:,:,:]; epsilon_zz = epsilon_zz[:,:,:] - epsilon_yy = epsilon_yy[:,1:end-1,:]; + epsilon_yy = epsilon_yy[:,1:end-1,:] elseif yBC == "PECPMC" epsilon_xx = epsilon_xx[:,2:end,:]; epsilon_xy = epsilon_xy[:,2:end,:]; epsilon_xz = epsilon_xz[:,2:end,:] epsilon_yx = epsilon_yx[:,2:end,:]; epsilon_yy = epsilon_yy[:,1:end-1,:]; epsilon_yz = epsilon_yz[:,2:end,:] @@ -725,7 +726,7 @@ export mesti_subpixel_smoothing #epsilon_xx = epsilon_xx[:,:,:]; epsilon_xy = epsilon_xy[:,:,:]; epsilon_xz = epsilon_xz[:,:,:] #epsilon_yx = epsilon_yx[:,:,:]; epsilon_yy = epsilon_yy[:,:,:]; epsilon_yz = epsilon_yz[:,:,:] #epsilon_zx = epsilon_zx[:,:,:]; epsilon_zy = epsilon_zy[:,:,:]; epsilon_zz = epsilon_zz[:,:,1:end-1] - epsilon_zz = epsilon_zz[:,:,1:end-1]; + epsilon_zz = epsilon_zz[:,:,1:end-1] elseif zBC == "PECPMC" epsilon_xx = epsilon_xx[:,:,2:end]; epsilon_xy = epsilon_xy[:,:,2:end]; epsilon_xz = epsilon_xz[:,:,2:end] epsilon_yx = epsilon_yx[:,:,2:end]; epsilon_yy = epsilon_yy[:,:,2:end]; epsilon_yz = epsilon_yz[:,:,2:end] @@ -739,7 +740,7 @@ export mesti_subpixel_smoothing end - function pick_epsilon_2d_TM(epsilon_xx::Union{Array{Int64},Array{Float64},Array{ComplexF64}}, yBC::String, zBC::String) + function pick_epsilon_2d_TM(epsilon_xx::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, yBC::String, zBC::String) # Pick the right site of components of epsilon based on Yee cell and boundary condition if yBC == "PEC" @@ -765,7 +766,7 @@ export mesti_subpixel_smoothing end - function pick_inv_epsilon_2d_TE(inv_epsilon_Ey_site::Union{Array{Int64},Array{Float64},Array{ComplexF64}}, inv_epsilon_Ez_site::Union{Array{Int64},Array{Float64},Array{ComplexF64}}, inv_epsilon_Eo_site::Union{Array{Int64},Array{Float64},Array{ComplexF64}}, yBC::String, zBC::String) + function pick_inv_epsilon_2d_TE(inv_epsilon_Ey_site::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, inv_epsilon_Ez_site::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, inv_epsilon_Eo_site::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, yBC::String, zBC::String) # Pick the right site of components of epsilon based on Yee cell and boundary condition inv_epsilon_yy = inv_epsilon_Ey_site[:,:,2,2] @@ -773,7 +774,7 @@ export mesti_subpixel_smoothing inv_epsilon_yz = inv_epsilon_Eo_site[:,:,2,3] if yBC == "PEC" - inv_epsilon_yy = inv_epsilon_yy[:,:]; + inv_epsilon_yy = inv_epsilon_yy[:,:] inv_epsilon_yz = inv_epsilon_yz[2:end,:] inv_epsilon_zz = inv_epsilon_zz[2:end,:] elseif yBC == "PMC" @@ -810,7 +811,7 @@ export mesti_subpixel_smoothing return inv_epsilon_yy, inv_epsilon_zz, inv_epsilon_yz end - function Kottke_smoothing(vol_frac::Float64, n0::SVector{3, Float64}, epsilon_object::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, epsilon_voxel::Union{Array{<:Real},Array{<:Complex}}) + function Kottke_smoothing(vol_frac::Float64, n0::SVector{3, Float64}, epsilon_object::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}, epsilon_voxel::Union{Array{<:Int},Array{<:Real},Array{<:Complex}}) Scomp = @SMatrix rand(3,3-1) # directions complementary to n12; works even for K = 1 Stemp = [n0 Scomp] # SMatrix{K,K} S = qr(Stemp).Q # nonallocating; 1st column is normalized n12 diff --git a/src/mumps3_interface.jl b/src/mumps3_interface.jl index 9e76701..296276c 100644 --- a/src/mumps3_interface.jl +++ b/src/mumps3_interface.jl @@ -4,14 +4,15 @@ # None of these functions change the JOB parameter. export invoke_mumps!, +set_keep!, set_cntl!, set_icntl!, set_job!, provide_matrix!, provide_rhs!, get_rhs!, get_rhs, get_sol!, get_sol, set_schur_centralized_by_column!, -get_schur_complement!, get_schur_complement - +get_schur_complement!, get_schur_complement, +finalize! """ invoke_mumps_unsafe!(mumps) @@ -81,31 +82,46 @@ end """ - set_icntl!(mumps,i,val; [displaylevel=1]) + set_keep!(mumps,i,val; [displaylevel=0]) + +Set the keep parameters according to KEEP[i]=val + +See also: [`display_keep`](@ref) +""" +function set_keep!(mumps::Mumps,i::Int,val::Int; displaylevel=0) + keep = mumps.keep + mumps.keep = (keep[1:i-1]...,convert(MUMPS_INT,val),keep[i+1:end]...) + #displaylevel>0 ? display_keep(stdout,mumps.keep,i,val) : nothing + return nothing +end + + +""" + set_icntl!(mumps,i,val; [displaylevel=0]) Set the integer control parameters according to ICNTL[i]=val See also: [`display_icntl`](@ref) """ -function set_icntl!(mumps::Mumps,i::Int,val::Int; displaylevel=mumps.icntl[4]-1) +function set_icntl!(mumps::Mumps,i::Int,val::Int; displaylevel=0) icntl = mumps.icntl mumps.icntl = (icntl[1:i-1]...,convert(MUMPS_INT,val),icntl[i+1:end]...) - displaylevel>0 ? display_icntl(stdout,mumps.icntl,i,val) : nothing + #displaylevel>0 ? display_icntl(stdout,mumps.icntl,i,val) : nothing return nothing end """ - set_cntl!(mumps,i,val; [displaylevel=1]) + set_cntl!(mumps,i,val; [displaylevel=0]) Set the real/complex control parameters according to CNTL[i]=val See also: [`display_cntl`](@ref) """ -function set_cntl!(mumps::Mumps{TC,TR},i::Int,val::Float64; displaylevel=mumps.icntl[4]-1) where {TC,TR} +function set_cntl!(mumps::Mumps{TC,TR},i::Int,val::Float64; displaylevel=0) where {TC,TR} cntl = mumps.cntl mumps.cntl = (cntl[1:i-1]...,convert(TR,val),cntl[i+1:end]...) - displaylevel>0 ? display_cntl(stdout,mumps.cntl,i,val) : nothing + #displaylevel>0 ? display_cntl(stdout,mumps.cntl,i,val) : nothing return nothing end diff --git a/src/mumps3_printing.jl b/src/mumps3_printing.jl index ecd7b1f..50ec660 100644 --- a/src/mumps3_printing.jl +++ b/src/mumps3_printing.jl @@ -1,4 +1,4 @@ -export display_icntl +export display_icntl, display_cntl, display_keep function Base.show(io::IO,mumps::Mumps{TC,TR}) where {TC,TR} print(io,"Mumps{$TC,$TR}: ") @@ -84,6 +84,7 @@ function Base.show(io::IO,mumps::Mumps{TC,TR}) where {TC,TR} end end + """ display_icntl(mumps) @@ -415,7 +416,6 @@ function display_icntl(io::IO,icntl,i,val) end - """ display_cntl(mumps) @@ -449,3 +449,24 @@ function display_cntl(io::IO,cntl,i,val) print(io,"\n") return nothing end + + +""" + display_keep(mumps) + +Show the KEEP array of `mumps` + +See also: [`set_keep!`](@ref) +""" +display_keep(io::IO,mumps::Mumps) = display_keep(io,mumps.mumps.keep) +function display_keep(io::IO,keep) + for i ∈ eachindex(keep) + display_keep(io,keep,i,icntl[i]) + end +end +function display_keep(io::IO,keep,i,val) + if i == 401 + print(io,"keep(401): ") + print(io,"$val") + end +end \ No newline at end of file diff --git a/test/interface_t_r_test.jl b/test/interface_t_r_test.jl index b495924..75b195d 100644 --- a/test/interface_t_r_test.jl +++ b/test/interface_t_r_test.jl @@ -11,7 +11,7 @@ syst.wavelength = 1 # vacuum wavelength syst.dx = 1/resolution # grid size # Use optimized PML parameters for this resolution to reduce error -zpml = get_optimal_PML(syst.wavelength/syst.dx) +zpml = mesti_optimal_pml_params(syst.wavelength/syst.dx) zpml.npixels = 30 syst.zPML = [zpml] k0dx = 2*pi/syst.wavelength*syst.dx # Dimensionless frequency k0*dx diff --git a/test/unitary_test.jl b/test/unitary_test.jl index 0d8f58c..21cd638 100644 --- a/test/unitary_test.jl +++ b/test/unitary_test.jl @@ -19,7 +19,7 @@ nx_Ey = nx; ny_Ey = ny; nz_Ey = nz -1 nx_Ez = nx; ny_Ez = ny; nz_Ez = nz # Use optimized PML parameters for this resolution to reduce error -zpml = get_optimal_PML(syst.wavelength/syst.dx) +zpml = mesti_optimal_pml_params(syst.wavelength/syst.dx) zpml.npixels = 25 syst.zPML = [zpml]