Skip to content

Commit

Permalink
Require the pml.direction in mesti() in default
Browse files Browse the repository at this point in the history
  • Loading branch information
hochunlin committed Jun 10, 2024
1 parent 94b394a commit 0d6c294
Show file tree
Hide file tree
Showing 10 changed files with 76 additions and 25 deletions.
4 changes: 3 additions & 1 deletion examples/2d_metalens_focusing_via_angular_spectrum_propagation/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,9 @@ syst = Syst()
syst.epsilon_xx = epsilon_syst
syst.wavelength = wavelength
syst.dx = dx
syst.PML = [PML(nPML)] # Number of PML pixels (on all four sides)
pml = PML(nPML)
pml.direction = "all"
syst.PML = [pml] # PML on all four sides

# In mesti(), B_struct.pos = [m1, n1, h, w] specifies the position of a
# block source, where (m1, n1) is the index of the smaller-(y,z) corner,
Expand Down
7 changes: 5 additions & 2 deletions ...via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,9 @@
"syst.epsilon_xx = epsilon_syst\n",
"syst.wavelength = wavelength\n",
"syst.dx = dx\n",
"syst.PML = [PML(nPML)] # Number of PML pixels (on all four sides)\n",
"pml = PML(nPML)\n",
"pml.direction = \"all\"\n",
"syst.PML = [pml] # PML on all four sides\n",
"\n",
"# In mesti(), B_struct.pos = [m1, n1, h, w] specifies the position of a\n",
"# block source, where (m1, n1) is the index of the smaller-(y,z) corner,\n",
Expand Down Expand Up @@ -306,7 +308,7 @@
"id": "5443bd6a",
"metadata": {},
"source": [
"## Calculate the field on the left side, inside and after the metalens"
"# Calculate the field on the left side, inside and after the metalens"
]
},
{
Expand All @@ -316,6 +318,7 @@
"metadata": {},
"outputs": [],
"source": [
"# Calculate the field on the left side, inside and after the metalens\n",
"else\n",
" # Use mesti() to compute the field inside the metalens.\n",
" field_inside_metalens, _ = mesti(syst, [B_struct], opts)\n",
Expand Down
4 changes: 3 additions & 1 deletion ...ng_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,9 @@ syst = Syst()
syst.epsilon_xx = epsilon_syst
syst.wavelength = wavelength
syst.dx = dx
syst.PML = [PML(nPML)] # Number of PML pixels (on all four sides)
pml = PML(nPML)
pml.direction = "all"
syst.PML = [pml] # PML on all four sides

# In mesti(), B_struct.pos = [m1, n1, h, w] specifies the position of a
# block source, where (m1, n1) is the index of the smaller-(y,z) corner,
Expand Down
4 changes: 3 additions & 1 deletion examples/2d_reflection_matrix_Gaussian_beams/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,9 @@ C = "transpose(B)";
The scattering matrix is given by S = C\*inv(A)\*B - D, with D = C\*inv(A<sub>0</sub>)\*B - S<sub>0</sub> where A<sub>0</sub> is a reference system for which its scattering matrix S<sub>0</sub> is known. We consider A<sub>0</sub> to be a homogeneous space with no scatterers, for which the reflection matrix S<sub>0</sub> is zero.

```julia
syst.PML = [PML(pml_npixels)] # Put PML on all four sides
pml = PML(pml_npixels)
pml.direction = "all"
syst.PML = [pml] # PML on all four sides

# For a homogeneous space, the length of the simulation domain doesn't
# matter, so we choose a minimal thickness of nz_Ex_temp = n_source + pml_npixels
Expand Down
4 changes: 3 additions & 1 deletion examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -285,7 +285,9 @@
"metadata": {},
"outputs": [],
"source": [
"syst.PML = [PML(pml_npixels)] # Put PML on all four sides\n",
"pml = PML(pml_npixels)\n",
"pml.direction = \"all\"\n",
"syst.PML = [pml] # PML on all four sides\n",
"\n",
"# For a homogeneous space, the length of the simulation domain doesn't\n",
"# matter, so we choose a minimal thickness of nz_Ex_temp = n_source + pml_npixels\n",
Expand Down
4 changes: 3 additions & 1 deletion examples/2d_reflection_matrix_Gaussian_beams/reflection_matrix_Gaussian_beams.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,9 @@ C = "transpose(B)"
# scattering matrix S_0 is known. We consider A_0 to be a homogeneous space
# with no scatterers, for which the reflection matrix S_0 is zero.

syst.PML = [PML(pml_npixels)] # Put PML on all four sides
pml = PML(pml_npixels)
pml.direction = "all"
syst.PML = [pml] # PML on all four sides

# For a homogeneous space, the length of the simulation domain doesn't
# matter, so we choose a minimal thickness of nz_Ex_temp = n_source + pml_npixels
Expand Down
2 changes: 1 addition & 1 deletion src/mesti2s.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -787,7 +787,7 @@ function mesti2s(syst::Syst, input::Union{channel_type, channel_index, wavefront

if ~(length(syst.zPML)==2) && two_sided
if isdefined(syst.zPML, :side) && syst.zPML.side != "both"
throw(ArgumentError("syst.zPML.direction should be assigned as \"both\", if given, when the user just provided only one PML object is provided for a two-sided geometry."))
throw(ArgumentError("syst.zPML.side should be assigned as \"both\", if given, when the user just provided only one PML object is provided for a two-sided geometry."))
end
# Apply the same zPML on both sides
syst.zPML = [deepcopy(syst.zPML[1]), deepcopy(syst.zPML[1])]
Expand Down
2 changes: 1 addition & 1 deletion src/mesti_build_fdfd_matrix.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -560,7 +560,7 @@ end
Boundary condition in z direction, analogous to yBC.
yPML (a two-element PML vector; optional):
Parameters for perfectly matched layer (PML) in y direction.
xPML = [PML_left, PML_right]
yPML = [PML_left, PML_right]
If users do not want to put PML on either side, just set the field "npixels" = 0.
For example, PML_left.npixels = 0: no PML on the left.
In each case, PML is a PML structure with the following fields:
Expand Down
66 changes: 53 additions & 13 deletions src/mesti_main.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,12 +191,18 @@ end
Number of PML pixels.
Note this is within syst.epsilon_ij (i = x,y,z and j = x,y,z),
not in addition to.
direction (string; optional):
direction (string; required):
Direction(s) where PML is placed. Available choices are (case-insensitive):
"all" - (default) PML in x, y, and z directions for 3D and PML in y and z directions for 2D TM
"x" - PML in x direction
"y" - PML in y direction
"z" - PML in z direction
"z" - PML in z direction
"xy" - PML in x and y directions
"xz" - PML in x and z directions
"yz" - PML in y and z directions
"yx" - PML in x and y directions
"zx" - PML in x and z directions
"zy" - PML in y and z directions
side (string; optional):
Side(s) where PML is placed.Available choices are (case-insensitive):
"both" - (default) PML on both sides
Expand Down Expand Up @@ -845,8 +851,10 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC
end

# Defaults to no PML anywhere
if ~isdefined(syst, :PML)
syst.PML = [PML(0)]
if ~isdefined(syst, :PML)
pml = PML(0)
pml.direction = "all"
syst.PML = [pml]
elseif isa(syst.PML, PML)
# If user specifies PML structure instead of vector of PML strcture, transform it to vector of PML strcture.
syst.PML = [syst.PML]
Expand All @@ -869,7 +877,7 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC
end

# Number of PML pixels must be given
# Other fields are optional and will be checked in mesti_build_fdfd_matrix()
# Other PML parameters are optional and will be checked in mesti_build_fdfd_matrix()
if ~isdefined(PML_ii, :npixels)
throw(ArgumentError("syst.PML[$(ii)] must contain field \"npixels\"."))
end
Expand All @@ -878,11 +886,11 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC
use_PML = true
end

# If PML is specified, we put it on both x, y, and z directions by default for 3D, and both y and z directions by default for 2D
# Direction(s) where PML is placed must be given
if ~isdefined(PML_ii, :direction)
PML_ii.direction = "all"
elseif ~(lowercase(PML_ii.direction) in ["all", "x", "y", "z"])
throw(ArgumentError("syst.PML[$(ii)].direction = \"$(PML_ii.direction)\" is not a supported option; use \"all\", \"x\", \"y\", or \"z\"."))
throw(ArgumentError("syst.PML[$(ii)] must contain field \"direction\"."))
elseif ~(lowercase(PML_ii.direction) in ["all", "x", "y", "z", "xy", "xz", "yz", "yx", "zx", "zy"])
throw(ArgumentError("syst.PML[$(ii)].direction = \"$(PML_ii.direction)\" is not a supported option; use \"all\", \"x\", \"y\", \"z\", \"xy\", \"xz\", \"yz\", \"yx\", \"zx\", or \"zy\"."))
end

# If PML is specified, we put it on both sides by default
Expand Down Expand Up @@ -918,6 +926,9 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC
end
end
elseif lowercase(PML_ii.direction) == "x"
if use_2D_TM
@warn "syst.PML[$(ii)].direction = \"x\" would not apply to 2D cases; will be ignored."
end
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [1,2]
elseif PML_ii.side == "-"
Expand All @@ -933,16 +944,45 @@ function mesti(syst::Syst, B::Union{SparseMatrixCSC{Int64,Int64},SparseMatrixCSC
else # PML_ii.side = "+"
ind_ii = 4
end
else # PML_ii.direction == "z"
elseif lowercase(PML_ii.direction) == "z"
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [5,6]
elseif lowercase(PML_ii.side) == "-"
ind_ii = 5
else # PML_ii.side = "+"
ind_ii = 6
end
end

end
elseif lowercase(PML_ii.direction) == "xy" || lowercase(PML_ii.direction) == "yx"
if use_2D_TM
@warn "syst.PML[$(ii)].direction = \"$(PML_ii.side)\"; in 2D cases, PML in x-direction will not be applied and ignored."
end
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [1,2,3,4]
elseif lowercase(PML_ii.side) == "-"
ind_ii = [1,3]
else # PML_ii.side = "+"
ind_ii = [2,4]
end
elseif lowercase(PML_ii.direction) == "xz" || lowercase(PML_ii.direction) == "zx"
if use_2D_TM
@warn "syst.PML[$(ii)].direction = \"$(PML_ii.side)\"; in 2D cases, PML in x-direction will not be applied and ignored."
end
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [1,2,5,6]
elseif lowercase(PML_ii.side) == "-"
ind_ii = [1,5]
else # PML_ii.side = "+"
ind_ii = [2,6]
end
else # lowercase(PML_ii.direction) == "yz" || lowercase(PML_ii.direction) == "zy"
if lowercase(PML_ii.side) == "both" # low & high
ind_ii = [3,4,5,6]
elseif lowercase(PML_ii.side) == "-"
ind_ii = [3,5]
else # PML_ii.side = "+"
ind_ii = [4,6]
end
end
# Specify PML at those locations
for jj = 1:length(ind_ii)
ind_side = ind_ii[jj]
Expand Down
4 changes: 1 addition & 3 deletions src/mesti_matrix_solver.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1030,9 +1030,7 @@ function MUMPS_analyze_and_factorize(A::Union{SparseMatrixCSC{Int64, Int64},Spar
if opts.use_L0_threads
# Utilize L0-threads feature.
# This typically improves the time performance, but marginally increases the memory usage in full multithread.
set_icntl!(id,48,1;displaylevel=0)
else
set_icntl!(id,48,0;displaylevel=0)
set_keep!(id,401,1)
end
set_job!(id,1) # what to do: analysis
if opts.use_given_ordering
Expand Down

2 comments on commit 0d6c294

@hochunlin
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/108620

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.4.4 -m "<description of version>" 0d6c2941a5699ec77956afce358768d80b39efcd
git push origin v0.4.4

Please sign in to comment.