Add build_EVP method to IVP solver

dedalus/core/field.py

@@ -255,7 +255,7 @@ def expression_matrices(self, subproblem, vars, **kw):
         """Build expression matrices for a specific subproblem and variables."""
         raise NotImplementedError()
 
-    def frechet_differential(self, variables, perturbations):
+    def frechet_differential(self, variables, perturbations, backgrounds=None):
         """
         Compute Frechet differential with respect to specified variables/perturbations.
 
@@ -265,29 +265,32 @@ def frechet_differential(self, variables, perturbations):
             Variables to differentiate around.
         perturbations : list of Field objects
             Perturbation directions for each variable.
+        backgrounds : list of Field objects, optional
+            Backgrounds for each variable. Default: variables.
 
         Notes
         -----
-        This method symbolically computes the functional directional derivative in the
-        direction of the specified perturbations:
-            sum_{vars, perts} lim_{ε -> 0} d/dε X(var + ε * pert)
+        This method symbolically computes the functional directional derivative around the
+        specified backgrounds in the direction of the specified perturbations:
+            F'(X0).X1 = lim_{ε -> 0} d/dε F(X0 + ε*X1)
         The result is a linear operator acting on the perturbations with NCCs that
-        depend on the original variables.
+        depend on the backgrounds.
         """
         dist = self.dist
         tensorsig = self.tensorsig
         dtype = self.dtype
-        # Perturbation variable
-        ep = Field(dist=dist, name='ep', dtype=dtype)
-        # Add differentials for each variable
-        diff = 0
+        # Compute differential
+        epsilon = Field(dist=dist, dtype=dtype)
+        diff = self
         for var, pert in zip(variables, perturbations):
-            diff_var = self.replace(var, var + ep*pert)
-            diff_var = diff_var.sym_diff(ep)
-            if diff_var:
-                diff_var = Operand.cast(diff_var, self.dist, tensorsig=tensorsig, dtype=dtype)
-                diff_var = diff_var.replace(ep, 0)
-                diff += diff_var
+            diff = diff.replace(var, var + epsilon*pert)
+        diff = diff.sym_diff(epsilon)
+        diff = Operand.cast(diff, self.dist, tensorsig=tensorsig, dtype=dtype)
+        diff = diff.replace(epsilon, 0)
+        # Replace backgrounds
+        if backgrounds:
+            for var, bg in zip(variables, backgrounds):
+                diff = diff.replace(var, bg)
         return diff
 
     @property

dedalus/core/problems.py

@@ -354,6 +354,56 @@ def _build_matrix_expressions(self, eqn):
         logger.debug(f"  L: {L}")
         logger.debug(f"  F: {F}")
 
+    def build_EVP(self, eigenvalue=None, backgrounds=None, perturbations=None, **kw):
+        """
+        Create an eigenvalue problem from an initial value problem.
+
+        Parameters
+        ----------
+        eigenvalue : Field, optional
+            Eigenvalue field.
+        backgrounds : list of Fields, optional
+            Background fields for linearizing the RHS. Default: the IVP variables.
+        perturbations : list of Fields, optional
+            Perturbation fields for the EVP. Default: copies of IVP variables.
+
+        Notes
+        -----
+        This method converts time-independent IVP equations of the form
+            M.dt(X) + L.X = F(X)
+        to EVP equations as
+            λ*M.X1 + L.X1 - F'(X0).X1 = 0.
+        If backgrounds (X0) are not specified, the IVP variables (X) are used.
+        """
+        # Create eigenvalue problem for perturbations
+        variables = self.variables
+        if eigenvalue is None:
+            eigenvalue = self.dist.Field(name='λ')
+        if perturbations is None:
+            perturbations = [var.copy() for var in variables]
+        EVP = EigenvalueProblem(perturbations, eigenvalue, **kw)
+        # Convert equations from IVP
+        for eqn in self.equations:
+            # Extract IVP expressions
+            M, L = eqn['LHS'].split(operators.TimeDerivative)
+            F = eqn['RHS']
+            # Convert M@dt(X) to λ*M@Y
+            if M:
+                M = M.replace(operators.TimeDerivative, lambda x: eigenvalue*x)
+                for var, pert in zip(variables, perturbations):
+                    M = M.replace(var, pert)
+            # Convert L@X to L@Y
+            if L:
+                for var, pert in zip(variables, perturbations):
+                    L = L.replace(var, pert)
+            # Take Frechet differential of F(X)
+            if F:
+                if F.has(self.time):
+                    raise UnsupportedEquationError("Cannot convert time-dependent IVP to EVP.")
+                F = F.frechet_differential(variables=variables, perturbations=perturbations, backgrounds=backgrounds)
+            EVP.add_equation((M + L - F, 0))
+        return EVP
+
 
 class EigenvalueProblem(ProblemBase):
     """
@@ -372,7 +422,7 @@ class EigenvalueProblem(ProblemBase):
     Notes
     -----
     This class supports linear eigenvalue problems of the form:
-        s*M.X + L.X = 0
+        λ*M.X + L.X = 0
     The LHS terms must be linear in the specified variables and affine in the eigenvalue.
     The RHS must be zero.
     """