The following is the Discrete Spectral Vorticity Evolution PDE for incompressible flow: $$ \frac{\partial \Omega_{pq}}{\partial t} = \nu \left( \frac{\partial^2 \Omega_{pq}}{\partial x^2} + \frac{\partial^2 \Omega_{pq}}{\partial y^2} \right) - \mathrm{fft2} \left( u_{ij} \frac{\partial \omega_{ij}}{\partial x} + v_{ij} \frac{\partial \omega_{ij}}{\partial y} \right) $$ After significant manipulation, we can obtain the following equivalent form, where $\Xi_{pq} = \exp \left[ \nu \left( k_p^2 + k_q^2 \right) t \right]$ is an integrating factor $$ \frac{\partial}{\partial t} \Big[ \Xi_{pq} \Omega_{pq} \Big] = -\underline{i} \Xi_{pq} \Big[ k_p \cdot \mathrm{fft2} \left(u_{ij} \omega_{ij} \right) + k_q \cdot \mathrm{fft2} \left( v_{ij} \omega_{ij} \right) \Big] $$ This secondary form is advantageous for numerical integration because we no longer need to consider both first (time derivative) and second-order (viscous term) stability, but brings with it a new problem.

The wavenumbers $k_p$ and $k_q$ are defined as follows (assuming a domain from $[0, 2 \pi]$, with $M$ points not including at $2 \pi$): $$ k_p = \begin{cases} p - 1, \quad \quad \quad \, p \leq \mathrm{floor} \left( \frac{M}{2} \right) + 1 \\ -M + p - 1, \quad p > \mathrm{floor} \left( \frac{M}{2} \right) + 1 \end{cases}$$ Because the wavenumbers are squared, our integrating factor $\Xi_{pq}$ is constantly increasing with time, and even though our viscosity constant $\nu$'s small size prevents an immediate blow-up, eventually we do start having overflow problems as time increases.

Is there any way to prevent this while still only maintaining first-order stability?

One idea I had was to use some type of logarithm-exponent trick, but that ran out of steam when I didn't know what to do with the LHS.

  • $\begingroup$ If you apply the product rule on the left and note that $\partial/\partial t \Xi = const \Xi$, you can see that you can just divide everything by $\Xi$ and the problem disappears. $\endgroup$ Jan 4 at 17:19
  • $\begingroup$ @WolfgangBangerth, that's just reverting back to the original form, where we have to deal with the second-order viscous stability. $\endgroup$ Jan 4 at 21:00
  • 1
    $\begingroup$ Integrate w.r.t. $t$ and divide through by $\Xi$ and you have the basic form of an exponential integrator. These are well-studied for these types of problems and only requires evaluating your integrating factor and related functions at times on the order of the integration step size na.math.kit.edu/download/papers/acta-final.pdf $\endgroup$
    – whpowell96
    Jan 4 at 22:44
  • $\begingroup$ @whpowell96, that seems promising, but I want to confirm. You are saying to transform the equation to and integrate the following? $$ \Omega_{pq} = \frac{1}{\Xi_{pq}(t)} \int_0^t - \underline{i} \Xi_{pq}(t) \left[ k_p \cdot \mathrm{fft2} (u_{ij} \omega_{ij} ) + k_q \cdot \mathrm{fft2} ( v_{ij} \omega_{ij}) \right] \, \, dt $$ I appreciate the link to the paper, but would it be possible to point me in the direction of a simple implementation? $\endgroup$ Jan 4 at 23:14
  • $\begingroup$ Equation 1.6 in the linked paper shows a very simple exponential Euler scheme and there is lots of literature (even on Wikipedia) on exponential Rosenbrock and Runge-Kutta schemes $\endgroup$
    – whpowell96
    Jan 5 at 20:16


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.