compute convolution term in its own function

2 files changed, 38 insertions, 16 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 3fbce07..13d9531 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -329,6 +329,37 @@ int ins_rhs(
   int kx,ky;
   int i;
 
+  // compute convolution term
+  ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
+
+  for(i=0; i<(2*K1+1)*(2*K2+1); i++){
+    out[i]=0;
+  }
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      if(kx!=0 || ky!=0){
+        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
+      }
+    }
+  }
+
+  return(0);
+}
+
+// convolution term in right side of convolution equation
+int ns_T(
+  _Complex double* u,
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  fft_vect fft1,
+  fft_vect fft2,
+  fft_vect ifft
+){
+  int kx,ky;
+  int i;
+
   // F(px/|p|*u)*F(qy*|q|*u)
   // init to 0
   for(i=0; i<N1*N2; i++){
@@ -340,8 +371,8 @@ int ins_rhs(
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-        fft1.fft[klookup(kx,ky,N1,N2)]=kx/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)];
-        fft2.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)];
+        fft1.fft[klookup(kx,ky,N1,N2)]=kx/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N1;
+        fft2.fft[klookup(kx,ky,N1,N2)]=ky*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
       }
     }
   }
@@ -364,8 +395,8 @@ int ins_rhs(
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       if(kx!=0 || ky!=0){
-        fft1.fft[klookup(kx,ky,N1,N2)]=ky/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)];
-        fft2.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)];
+        fft1.fft[klookup(kx,ky,N1,N2)]=ky/sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N1;
+        fft2.fft[klookup(kx,ky,N1,N2)]=kx*sqrt(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]/N2;
       }
     }
   }
@@ -381,18 +412,6 @@ int ins_rhs(
   // inverse fft
   fftw_execute(ifft.fft_plan);
 
-  // write out
-  for(i=0; i<(2*K1+1)*(2*K2+1); i++){
-    out[i]=0;
-  }
-  for(kx=-K1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      if(kx!=0 || ky!=0){
-        out[klookup(kx,ky,2*K1+1,2*K2+1)]=-4*M_PI*M_PI*nu*(kx*kx+ky*ky)*u[klookup(kx,ky,2*K1+1,2*K2+1)]+(*g)(kx,ky)+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)]/N1/N2;
-      }
-    }
-  }
-
   return(0);
 }
 
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 5ac8d7d..2ff4e3f 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -31,6 +31,9 @@ int ins_step( _Complex double* u, int K1, int K2, int N1, int N2, double nu, dou
 // right side of Irreversible Navier-Stokes equation
 int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2, double nu, _Complex double (*g)(int,int), fft_vect fft1, fft_vect fft2, fft_vect ifft);
 
+// convolution term in right side of equation
+int ns_T( _Complex double* u, int K1, int K2, int N1, int N2, fft_vect fft1, fft_vect fft2, fft_vect ifft);
+
 // compute alpha
 _Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int));