Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

mode solver preconditioner improvement and turning off incidence matrix #1928

Merged
merged 2 commits into from
Sep 28, 2024
+67 −16
Merged

mode solver preconditioner improvement and turning off incidence matrix #1928

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7abff0e
mode solver preconditioner improvement and turning off incidence matrix
weiliangjin2021 Aug 29, 2024
eec46bd
Apply matrix trimming after preconditioning
weiliangjin2021 Aug 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
83 changes: 67 additions & 16 deletions tidy3d/plugins/mode/solver.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,8 @@
TOL_TENSORIAL = 1e-6
# shift target neff by this value, both rel and abs, whichever results in larger shift
TARGET_SHIFT = 10 * fp_eps
# Preconditioner: "Jacobi" or "Material" based
PRECONDITIONER = "Material"
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Now I also include the previous preconditioning method (Jacobi) along with the new one. For non-magnetic material (mu=1 everywhere), the two are actually very similar. The new one performs better when both epsilon and mu can be huge.



class EigSolver(Tidy3dBaseModel):
Expand Down Expand Up @@ -90,7 +92,7 @@ def compute_modes(
angle_phi = mode_spec.angle_phi
omega = 2 * np.pi * freq
k0 = omega / C_0
enable_incidence_matrices = split_curl_scaling is not None or mu_cross is not None
enable_incidence_matrices = False # Experimental feature, always off for now

eps_formated = cls.format_medium_data(eps_cross)
eps_xx, eps_xy, eps_xz, eps_yx, eps_yy, eps_yz, eps_zx, eps_zy, eps_zz = eps_formated
Expand Down Expand Up @@ -331,6 +333,7 @@ def conductivity_model_for_pec(eps, threshold=0.9 * pec_val):
return eps

eps_tensor = conductivity_model_for_pec(eps_tensor)
mu_tensor = conductivity_model_for_pec(mu_tensor)

# Determine if ``eps`` and ``mu`` are diagonal or tensorial
off_diagonals = (np.ones((3, 3)) - np.eye(3)).astype(bool)
Expand Down Expand Up @@ -464,7 +467,7 @@ def any_pec(eps_vec, threshold=0.9 * np.abs(pec_val)):
"""Check if there are any PEC values in the given permittivity array."""
return np.any(np.abs(eps_vec) >= threshold)

if any(any_pec(i) for i in [eps_xx, eps_yy, eps_zz]):
if any(any_pec(i) for i in [eps_xx, eps_yy, eps_zz, mu_xx, mu_yy, mu_zz]):
enable_preconditioner = True

# Compute the matrix for diagonalization
Expand Down Expand Up @@ -493,10 +496,6 @@ def any_pec(eps_vec, threshold=0.9 * np.abs(pec_val)):
mat_dtype = cls.matrix_data_type(eps, mu, der_mats, mat_precision, is_tensorial=False)
mat = cls.type_conversion(mat, mat_dtype)

# Trim small values in single precision case
if mat_precision == "single":
cls.trim_small_values(mat, tol=fp_eps)

# Casting starting vector to target data type
vec_init = cls.type_conversion(vec_init, mat_dtype)

Expand All @@ -507,18 +506,59 @@ def any_pec(eps_vec, threshold=0.9 * np.abs(pec_val)):
mat = dnz * mat * dnz.T
vec_init = dnz * vec_init

# Denote the original eigenvalue problem as Ax = lambda x,
# with left and right preconditioners, we solve for the following eigenvalue problem,
# L A R y = lambda LR y, where x = R y
precon_left = None
precon_right = None
generalized_M = None # matrix in the generalized eigenvalue problem
if enable_preconditioner:
precon = sp.diags(1 / mat.diagonal())
mat = mat * precon
else:
precon = None
if PRECONDITIONER == "Jacobi":
precon_right = sp.diags(1 / mat.diagonal())

elif PRECONDITIONER == "Material":

def conditional_inverted_vec(eps_vec, threshold=1):
"""Returns a diagonal sparse matrix whose i-th element in the diagonal
is |eps_i|^-1 if |eps_i|>threshold, and |eps_i| otherwise.
"""
abs_vec = np.abs(eps_vec)
return sp.spdiags(
np.where(abs_vec > threshold, 1.0 / abs_vec, abs_vec), [0], N, N
)

precon_left = sp.bmat(
[
[conditional_inverted_vec(mu_yy), None],
[None, conditional_inverted_vec(mu_xx)],
]
)
precon_right = sp.bmat(
[
[conditional_inverted_vec(eps_xx), None],
[None, conditional_inverted_vec(eps_yy)],
]
)
generalized_M = precon_right
mat = mat @ precon_right
if precon_left is not None:
generalized_M = precon_left @ generalized_M
mat = precon_left @ mat

# Trim small values in single precision case
if mat_precision == "single":
cls.trim_small_values(mat, tol=fp_eps)

if analyze_conditioning:
aca = mat.conjugate().T * mat
aac = mat * mat.conjugate().T
diff = aca - aac
print(spl.norm(diff, ord=np.inf), spl.norm(aca, ord=np.inf), spl.norm(aac, ord=np.inf))
print(spl.norm(diff, ord="fro"), spl.norm(aca, ord="fro"), spl.norm(aac, ord="fro"))
print(
f"inf-norm: A*A: {spl.norm(aca, ord=np.inf)}, AA*: {spl.norm(aac, ord=np.inf)}, nonnormality: {spl.norm(diff, ord=np.inf)}, relative nonnormality: {spl.norm(diff, ord=np.inf)/spl.norm(aca, ord=np.inf)}"
)
print(
f"fro-norm: A*A: {spl.norm(aca, ord='fro')}, AA*: {spl.norm(aac, ord='fro')}, nonnormality: {spl.norm(diff, ord='fro')}, relative nonnormality: {spl.norm(diff, ord='fro')/spl.norm(aca, ord='fro')}"
)
weiliangjin2021 marked this conversation as resolved.
Show resolved Hide resolved

# preprocess basis modes
basis_vecs = None
Expand All @@ -541,7 +581,7 @@ def any_pec(eps_vec, threshold=0.9 * np.abs(pec_val)):
vec_init,
guess_value=eig_guess,
mode_solver_type=mode_solver_type,
M=precon,
M=generalized_M,
)
else:
vals, vecs = cls.solver_eigs_relative(
Expand All @@ -550,7 +590,7 @@ def any_pec(eps_vec, threshold=0.9 * np.abs(pec_val)):
vec_init,
guess_value=eig_guess,
mode_solver_type=mode_solver_type,
M=precon,
M=generalized_M,
basis_vecs=basis_vecs,
)
neff, keff = cls.eigs_to_effective_index(vals, mode_solver_type)
Expand All @@ -562,8 +602,8 @@ def any_pec(eps_vec, threshold=0.9 * np.abs(pec_val)):

E, H = None, None
if basis_E is None:
if enable_preconditioner:
vecs = precon * vecs
if precon_right is not None:
vecs = precon_right * vecs

if enable_incidence_matrices:
vecs = dnz.T * vecs
Expand Down Expand Up @@ -744,6 +784,17 @@ def solver_eigs(
values, vectors = spl.eigs(
mat, k=num_modes, sigma=guess_value, tol=TOL_EIGS, v0=vec_init, M=M
)

# for i, eig_i in enumerate(values):
# vec = vectors[:, i]
# rhs = vec
# if M is not None:
# rhs = M @ rhs
# eig_from_vec = (vec.T @ (mat @ vec)) / (vec.T @ rhs)
# residue = np.linalg.norm(mat @ vec - eig_i * rhs) / np.linalg.norm(vec)
# print(
# f"{i}-th eigenvalue: {eig_i}, referred from eigenvectors: {eig_from_vec}, relative residual: {residue}."
# )
return values, vectors

@classmethod
Expand Down
Loading