From d9d19dc2e77c18fe0a8ec0299d830b324abfdd2d Mon Sep 17 00:00:00 2001 From: "Keaton J. Burns" Date: Thu, 14 Dec 2023 22:46:55 -0500 Subject: [PATCH] Change left evec normalization and add problem formulation docs describing this --- dedalus/core/solvers.py | 45 ++++++++-------- docs/pages/problem_formulations.rst | 80 +++++++++++++++++++++++++++++ docs/pages/user_guide.rst | 5 +- 3 files changed, 105 insertions(+), 25 deletions(-) create mode 100644 docs/pages/problem_formulations.rst diff --git a/dedalus/core/solvers.py b/dedalus/core/solvers.py index 10553e6a..9e3eb84f 100644 --- a/dedalus/core/solvers.py +++ b/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): """ diff --git a/docs/pages/problem_formulations.rst b/docs/pages/problem_formulations.rst new file mode 100644 index 00000000..864f960b --- /dev/null +++ b/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. + diff --git a/docs/pages/user_guide.rst b/docs/pages/user_guide.rst index 118d9d45..ee43151c 100644 --- a/docs/pages/user_guide.rst +++ b/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: