diff --git a/examples/2d_metalens_focusing_via_angular_spectrum_propagation/README.md b/examples/2d_metalens_focusing_via_angular_spectrum_propagation/README.md
index 31f4c4a..e97b9b3 100644
--- a/examples/2d_metalens_focusing_via_angular_spectrum_propagation/README.md
+++ b/examples/2d_metalens_focusing_via_angular_spectrum_propagation/README.md
@@ -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,
diff --git a/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.ipynb b/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.ipynb
index bb99700..034f89a 100644
--- a/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.ipynb
+++ b/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.ipynb
@@ -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",
@@ -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"
]
},
{
@@ -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",
diff --git a/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.jl b/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.jl
index 4d3e53f..d5e0395 100644
--- a/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.jl
+++ b/examples/2d_metalens_focusing_via_angular_spectrum_propagation/metalens_focusing_via_angular_spectrum_propagation.jl
@@ -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,
diff --git a/examples/2d_reflection_matrix_Gaussian_beams/README.md b/examples/2d_reflection_matrix_Gaussian_beams/README.md
index fb28e25..967e7e8 100644
--- a/examples/2d_reflection_matrix_Gaussian_beams/README.md
+++ b/examples/2d_reflection_matrix_Gaussian_beams/README.md
@@ -195,7 +195,9 @@ C = "transpose(B)";
The scattering matrix is given by S = C\*inv(A)\*B - D, with D = C\*inv(A0)\*B - S0 where A0 is a reference system for which its scattering matrix S0 is known. We consider A0 to be a homogeneous space with no scatterers, for which the reflection matrix S0 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
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 4ae16d8..9f505cf 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
@@ -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",
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 3259c83..8db8501 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
@@ -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
diff --git a/src/mesti2s.jl b/src/mesti2s.jl
index bd34065..042dfb5 100644
--- a/src/mesti2s.jl
+++ b/src/mesti2s.jl
@@ -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])]
diff --git a/src/mesti_build_fdfd_matrix.jl b/src/mesti_build_fdfd_matrix.jl
index fccb6d3..aa3202d 100644
--- a/src/mesti_build_fdfd_matrix.jl
+++ b/src/mesti_build_fdfd_matrix.jl
@@ -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:
diff --git a/src/mesti_main.jl b/src/mesti_main.jl
index 0cfe086..e21a9b7 100644
--- a/src/mesti_main.jl
+++ b/src/mesti_main.jl
@@ -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
@@ -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]
@@ -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
@@ -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
@@ -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 == "-"
@@ -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]
diff --git a/src/mesti_matrix_solver.jl b/src/mesti_matrix_solver.jl
index b629597..9f38471 100644
--- a/src/mesti_matrix_solver.jl
+++ b/src/mesti_matrix_solver.jl
@@ -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