From 5d92b56c01a79ecbb2df963b4732d0f12687cf75 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian@jauslin.org>
Date: Wed, 25 May 2022 21:42:21 -0400
Subject: Stabilize algorithm by enforcing real parts

---
 src/navier-stokes.c | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

(limited to 'src')

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 89bfbbd..e553022 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -495,7 +495,8 @@ int ns_T(
   fftw_execute(fft2.fft_plan);
   // write to ifft
   for(i=0;i<N1*N2;i++){
-    ifft.fft[i]=fft1.fft[i]*fft2.fft[i];
+    // control numerical truncation by taking imaginary part (fft1 and fft2 should be purely imaginary)
+    ifft.fft[i]=-(__imag__ fft1.fft[i])*(__imag__ fft2.fft[i]);
   }
 
   // F(py/|p|*u)*F(qx*|q|*u)
@@ -519,7 +520,8 @@ int ns_T(
   fftw_execute(fft2.fft_plan);
   // write to ifft
   for(i=0;i<N1*N2;i++){
-    ifft.fft[i]=ifft.fft[i]-fft1.fft[i]*fft2.fft[i];
+    // control numerical truncation by taking imaginary part (fft1 and fft2 should be purely imaginary)
+    ifft.fft[i]=ifft.fft[i]+(__imag__ fft1.fft[i])*(__imag__ fft2.fft[i]);
   }
 
   // inverse fft
-- 
cgit v1.2.3-54-g00ecf