Enforce T(u,-k)=T(u,k)^* in ns_T

1 files changed, 15 insertions, 2 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 1ab4a44..1980207 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -383,8 +383,7 @@ int ins_rhs(
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-	// enforce the reality of u by adding ifft.fft(k) and the conjugate of ifft.fft(-k)
-        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+g[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
+        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI/L/L*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+g[klookup(kx,ky,2*K1+1,2*K2+1)]+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
       }
     }
   }
@@ -458,6 +457,13 @@ int ns_T(
   // inverse fft
   fftw_execute(ifft.fft_plan);
 
+  // enforce T(u,-k)=T(u,k)^*
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      ifft.fft[klookup(kx,ky,N1,N2)]=(ifft.fft[klookup(kx,ky,N1,N2)]+conj(ifft.fft[klookup(-kx,-ky,N1,N2)]))/2;
+    }
+  }
+
   return(0);
 }
 
@@ -498,6 +504,13 @@ int ns_T_nofft(
     }
   }
 
+  // enforce T(u,-k)=T(u,k)^*
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      out[klookup(kx,ky,N1,N2)]=(out[klookup(kx,ky,N1,N2)]+conj(out[klookup(-kx,-ky,N1,N2)]))/2;
+    }
+  }
+
   return 0;
 }