mode solver preconditioner improvement and turning off incidence matrix

tidy3d/plugins/mode/solver.py

@@ -90,7 +90,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
@@ -331,6 +331,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)
@@ -464,7 +465,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
@@ -508,10 +509,26 @@ def any_pec(eps_vec, threshold=0.9 * np.abs(pec_val)):
             vec_init = dnz * vec_init
 
         if enable_preconditioner:
-            precon = sp.diags(1 / mat.diagonal())
-            mat = mat * precon
+
+            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)]]
+            )
+            # precon = sp.diags(1 / mat.diagonal())
+            precon = precon_left @ precon_right
+            mat = precon_left @ mat @ precon_right
         else:
             precon = None
+            precon_right = None
 
         if analyze_conditioning:
             aca = mat.conjugate().T * mat
@@ -563,7 +580,7 @@ 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
+                vecs = precon_right * vecs
 
             if enable_incidence_matrices:
                 vecs = dnz.T * vecs
@@ -744,6 +761,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 residue: {residue}."
+        #     )
         return values, vectors
 
     @classmethod