Change left evec normalization and add problem formulation docs descr…

…ibing this

dedalus/core/solvers.py

@@ -167,15 +167,6 @@ def print_subproblem_ranks(self, subproblems=None, target=0):
             A = L + target * M
             print(f"MPI rank: {self.dist.comm.rank}, subproblem: {i}, group: {sp.group}, matrix rank: {np.linalg.matrix_rank(A)}/{A.shape[0]}, cond: {np.linalg.cond(A):.1e}")
 
-    def _build_modified_left_eigenvectors(self):
-        sp = self.eigenvalue_subproblem
-        return - sp.pre_right@sp.M_min.T.conj()@self.left_eigenvectors
-
-    def _normalize_left_eigenvectors(self):
-        modified_left_eigenvectors = self._build_modified_left_eigenvectors()
-        norms = np.diag(modified_left_eigenvectors.T.conj() @ self.eigenvectors)
-        self.left_eigenvectors /= np.conj(norms)
-
     def solve_dense(self, subproblem, rebuild_matrices=False, left=False, normalize_left=True, **kw):
         """
         Perform dense eigenvector search for selected subproblem.
@@ -209,11 +200,14 @@ def solve_dense(self, subproblem, rebuild_matrices=False, left=False, normalize_
         eig_output = scipy.linalg.eig(A, b=B, left=left, **kw)
         # Unpack output
         if left:
-            self.eigenvalues, self.left_eigenvectors, pre_eigenvectors = eig_output
-            self.eigenvectors = sp.pre_right @ pre_eigenvectors
+            self.eigenvalues, pre_left_evecs, pre_right_evecs = eig_output
+            self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
+            self.left_eigenvectors = sp.pre_left.H @ pre_left_evecs
+            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).H @ pre_left_evecs
             if normalize_left:
-                self._normalize_left_eigenvectors()
-            self.modified_left_eigenvectors = self._build_modified_left_eigenvectors()
+                norms = np.diag(pre_left_evecs.T.conj() @ sp.M_min @ pre_right_evecs)
+                self.left_eigenvectors /= np.conj(norms)
+                self.modified_left_eigenvectors /= np.conj(norms)
         else:
             self.eigenvalues, pre_eigenvectors = eig_output
             self.eigenvectors = sp.pre_right @ pre_eigenvectors
@@ -256,15 +250,17 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
         A = sp.L_min
         B = - sp.M_min
         # Solve for the right eigenvectors
-        self.eigenvalues, pre_eigenvectors = scipy_sparse_eigs(A=A, B=B, N=N, target=target, matsolver=self.matsolver, **kw)
-        self.eigenvectors = sp.pre_right @ pre_eigenvectors
+        self.eigenvalues, pre_right_evecs = scipy_sparse_eigs(A=A, B=B, N=N, target=target, matsolver=self.matsolver, **kw)
+        self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
         if left:
             # Solve for the left eigenvectors
             # Note: this definition of "left eigenvectors" is consistent with the documentation for scipy.linalg.eig
-            self.left_eigenvalues, self.left_eigenvectors = scipy_sparse_eigs(A=A.getH(),
-                                                                              B=B.getH(),
-                                                                              N=N, target=np.conjugate(target),
-                                                                              matsolver=self.matsolver, **kw)
+            self.left_eigenvalues, pre_left_evecs = scipy_sparse_eigs(A=A.getH(), B=B.getH(),
+                                                                      N=N, target=np.conjugate(target),
+                                                                      matsolver=self.matsolver, **kw)
+            self.left_eigenvectors = sp.pre_left.H @ pre_left_evecs
+            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).H @ pre_left_evecs
+            # Check that eigenvalues match
             if not np.allclose(self.eigenvalues, np.conjugate(self.left_eigenvalues)):
                 if raise_on_mismatch:
                     raise RuntimeError("Conjugate of left eigenvalues does not match right eigenvalues. "
@@ -273,11 +269,14 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
                                        "solve_sparse().")
                 else:
                     logger.warning("Conjugate of left eigenvalues does not match right eigenvalues.")
-            # In absence of above warning, modified_left_eigenvectors forms a biorthogonal set with the right
-            # eigenvectors.
+                    if normalize_left:
+                        logger.warning("Cannot normalize left eigenvectors.")
+                        normalize_left = False
+            # Normalize
             if normalize_left:
-                self._normalize_left_eigenvectors()
-            self.modified_left_eigenvectors = self._build_modified_left_eigenvectors()
+                norms = np.diag(pre_left_evecs.T.conj() @ sp.M_min @ pre_right_evecs)
+                self.left_eigenvectors /= np.conj(norms)
+                self.modified_left_eigenvectors /= np.conj(norms)
 
     def set_state(self, index, subsystem):
         """

docs/pages/problem_formulations.rst

@@ -0,0 +1,80 @@
+Problem Formulations
+********************
+
+Dedalus parses all equations, including constraints and boundary conditions, into common forms based on the problem type.
+
+Linear Boundary-Value Problems (LBVPs)
+--------------------------------------
+
+Equations in linear boundary-value problems must all take the form:
+
+.. math::
+    L \cdot X = F,
+
+where :math:`X` is the state-vector of problem variables, :math:`L` are linear operators, and :math:`F` are inhomogeneous terms that are independent of :math:`X`.
+LBVPs are solved by explicitly evaluating the RHS and solving the sparse-matrix representation of the LHS to find :math:`X`.
+
+Initial-Value Problems (IVPs)
+-----------------------------
+
+Equations in initial value problems must all take the form:
+
+.. math::
+    M \cdot \partial_t X + L \cdot X = F(X,t),
+
+where :math:`X(t)` is the state-vector of problem variables, :math:`M` and :math:`L` are time-independent linear operators, and :math:`F` are inhomogeneous terms or nonlinear terms.
+Initial conditions :math:`X(t=0)` are set for the state, and the state is then evolved forward in time using mixed implicit-explicit timesteppers.
+During this process, the RHS is explicitly evaluated using :math:`X(t)` and the LHS is implicitly solved using the sparse-matrix representations of :math:`M` and :math:`L` to produce :math:`X(t+ \Delta t)`.
+
+Eigenvalue Problems (EVPs)
+--------------------------
+
+Equations in eigenvalue problems must all take the generalized form:
+
+.. math::
+    \lambda M \cdot X + L \cdot X = 0,
+
+where :math:`\lambda` is the eigenvalue, :math:`X` is the state-vector of problem variables, and :math:`M` and :math:`L` are linear operators. The standard *right eigenmodes* :math:`(\lambda_i, X_i)` are solved using the sparse-matrix representations of :math:`M` and :math:`L`, and satisfy:
+
+.. math::
+    \lambda_i M \cdot X_i + L \cdot X_i = 0.
+
+The *left eigenmodes* :math:`(\lambda_i, Y_i)` are solved (if requested) using the sparse-matrix representations of :math:`M` and :math:`L`, and satisfy:
+
+.. math::
+    \lambda_i Y_i^* \cdot M + Y_i^* \cdot L = 0.
+
+The left and right eigenmodes satisfy the generalized :math:`M`-orthogonality condition:
+
+.. math::
+    Y_i^* \cdot M \cdot X_j = 0 \quad \mathrm{if} \quad \lambda_i \neq \lambda_j.
+
+For convenience, we also provide *modified left eigenmodes* :math:`Z_i = M^* \cdot Y_i`.
+When the eigenvalues are nondegenerate, the left and modified left eigenvectors are rescaled by :math:`(X_i^* \cdot M^* \cdot Y_i)^{-1}` so that
+
+.. math::
+    Z_i^* \cdot X_j = \delta_{ij}.
+
+Nonlinear Boundary-Value Problems (NLBVPs)
+--------------------------------------
+
+Equations in nonlinear boundary-value problems must all take the form:
+
+.. math::
+    G(X) = H(X),
+
+where :math:`X` is the state-vector of problem variables and :math:`G` and :math:`H` are generic operators.
+All equations are immediately reformulated into the root-finding problem
+
+.. math::
+    F(X) = G(X) - H(X) = 0.
+
+NLBVPs are solved iteratively via Newton's method.
+The problem is reduced to a LBVP for an update :math:`\delta X` to the state using the symbolically-computed Frechet differential as:
+
+.. math::
+    F(X_n + \delta X) \approx 0 \quad \implies \quad \partial F(X_n) \cdot \delta X = - F(X_n)
+
+Each iteration entails reforming the LHS matrix, explicitly evaluating the RHS, solving the LHS to find :math:`\delta X`, and updating the new solution as :math:`X_{n+1} = X_n + \delta X`.
+The iterations proceed until convergence criteria are satisfied.
+

docs/pages/user_guide.rst

@@ -6,10 +6,11 @@ General user guide:
 .. toctree::
     :maxdepth: 1
 
-    changes_from_v2
-    configuration
+    problem_formulations
     performance_tips
+    configuration
     troubleshooting
+    changes_from_v2
 
 Specific how-to's: