Automatically enforce Hermitian symmetry by transforming to grid spac…

…e during IVP timestepping

dedalus/core/domain.py

@@ -43,6 +43,7 @@ def __init__(self, bases, grid_dtype=np.complex128, comm=None, mesh=None):
 
         # Iteratively set basis data types
         # (Grid-to-coefficient transforms proceed in the listed order)
+        self.grid_dtype = grid_dtype
         for basis in bases:
             grid_dtype = basis.set_dtype(grid_dtype)
 

dedalus/core/solvers.py

@@ -320,6 +320,8 @@ class InitialValueSolver:
         Timestepper to use in evolving initial conditions
     matsolver : matsolver class or name, optional
         Matrix solver routine (default set by config file).
+    enforce_real_cadence : int, optional
+        Iteration cadence for enforcing Hermitian symmetry on real variables (default: 100).
 
     Attributes
     ----------
@@ -340,7 +342,7 @@ class InitialValueSolver:
 
     """
 
-    def __init__(self, problem, timestepper, matsolver=None):
+    def __init__(self, problem, timestepper, matsolver=None, enforce_real_cadence=100):
 
         logger.debug('Beginning IVP instantiation')
 
@@ -350,6 +352,7 @@ def __init__(self, problem, timestepper, matsolver=None):
         self.problem = problem
         self.domain = domain = problem.domain
         self.matsolver = matsolver
+        self.enforce_real_cadence = enforce_real_cadence
         self._float_array = np.zeros(1, dtype=float)
         self.start_time = self.get_world_time()
 
@@ -506,17 +509,24 @@ def step(self, dt, trim=False):
         self.state.scatter()
         # Update iteration
         self.iteration += 1
+        # Enforce Hermitian symmetry for real variables
+        if self.domain.grid_dtype == np.float64:
+            # Enforce for as many iterations as timestepper uses internally
+            if self.iteration % self.enforce_real_cadence <= self.timestepper._history:
+                # Transform state variables to grid and back
+                for path in self.domain.dist.paths:
+                    path.increment(self.state.fields)
+                for path in self.domain.dist.paths[::-1]:
+                    path.decrement(self.state.fields)
         return dt
 
     def evolve(self, timestep_function):
         """Advance system until stopping criterion is reached."""
-
         # Check for a stopping criterion
         if np.isinf(self.stop_sim_time):
             if np.isinf(self.stop_wall_time):
                 if np.isinf(self.stop_iteration):
                     raise ValueError("No stopping criterion specified.")
-
         # Evolve
         while self.ok:
             dt = timestep_function()
@@ -535,7 +545,5 @@ def evaluate_handlers_now(self, dt, handlers=None):
         end_wall_time = end_world_time - self.start_time
         if handlers is None:
             handlers = self.evaluator.handlers
-
         self.evaluator.evaluate_handlers(handlers, timestep=dt, sim_time=self.sim_time, world_time=end_world_time, wall_time=end_wall_time, iteration=self.iteration)
 
-

dedalus/core/timesteppers.py

@@ -59,6 +59,7 @@ def __init__(self, pencil_length, domain):
 
         # Create deque for storing recent timesteps
         N = max(self.amax, self.bmax, self.cmax)
+        self._history = N
         self.dt = deque([0.]*N)
 
         # Create coefficient systems for multistep history
@@ -503,6 +504,8 @@ class RungeKuttaIMEX:
 
     """
 
+    _history = 0
+
     def __init__(self, pencil_length, domain):
 
         self.RHS = CoeffSystem(pencil_length, domain)

docs/pages/troubleshooting.rst

@@ -18,6 +18,16 @@ This is why you'll see these types of boundary conditions for the wall-normal ve
     problem.add_bc("right(w) = 0", condition="(nx != 0)")
     problem.add_bc("right(p) = 0", condition="(nx == 0)")
 
+Out of memory errors
+====================
+
+Spectral simulations with implicit timestepping can require a large amount of memory to store the LHS matrices and their factorizations.
+The best way to minimize the required memory is to minimize the LHS matrix size by using as few variables as possible and to minimize the LHS matrix bandwidth (see the :ref:`performance_tips` page).
+Beyond this, several of the Dedalus configuration options can be changed the minimize the simulation's memory footprint, potentially at the cost of reduced performance (see the :ref:`configuration` page).
+
+Reducing memory consumption in Dedalus is an ongoing effort.
+Any assistance with memory profiling and contributions reducing the code's memory footprint would be greatly appreciated!
+
 Maintaining Hermitian symmetry with real variables
 ==================================================
 
@@ -27,24 +37,9 @@ Subsequent transforms are then complex-to-complex.
 
 This strategy maintains the proper Hermitian symmetry for the majority of the modes, since only the non-negative kx half of the spectrum is evolved and the negative kx modes are computed from the conjugates of the positive kx modes whenever they are needed (i.e. for interpolation along x).
 However, modes with kx = 0 should have coefficients that are real after the x-transform, and should have Hermitian symmetry in e.g. ky if there's a subsequent Fourier transform in y.
-This is not enforced by Dedalus, so if the PDE has a linear instability for modes with kx = 0, then the imaginary part of the solution can grow tihout bound (the nonlinear terms are computed in grid space, so they only affect the real portion of the solution and can’t saturate a growing imaginary part).
+If these conditions are not enforced and the PDE has a linear instability for modes with kx = 0, then the imaginary part of the solution can potentially grow without bound (the nonlinear terms are computed in grid space, so they only affect the real portion of the solution and can’t saturate a growing imaginary part).
 
 A simple solution is to periodically transform all the state variables to grid space, which will project away any imaginary component that may have been building up during timestepping.
-This can be done with a simple statement in the main loop like:
-
-.. code-block:: python
-
-    if solver.iteration % 100 == 0:
-        for field in solver.state.fields:
-            field.require_grid_space()
-
-Out of memory errors
-====================
-
-Spectral simulations with implicit timestepping can require a large amount of memory to store the LHS matrices and their factorizations.
-The best way to minimize the required memory is to minimize the LHS matrix size by using as few variables as possible and to minimize the LHS matrix bandwidth (see the :ref:`performance_tips` page).
-Beyond this, several of the Dedalus configuration options can be changed the minimize the simulation's memory footprint, potentially at the cost of reduced performance (see the :ref:`configuration` page).
-
-Reducing memory consumption in Dedalus is an ongoing effort.
-Any assistance with memory profiling and contributions reducing the code's memory footprint would be greatly appreciated!
+This is done in Dedalus every 100 timesteps by default.
+This cadence can be modified via the ``enforce_real_cadence`` keyword when instantiating an IVP solver, and may need to be decreased in simulations with strong linear instabilities.