Initial commit

1 files changed, 299 insertions, 0 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
new file mode 100644
index 0000000..3def963
--- /dev/null
+++ b/src/navier-stokes.c
@@ -0,0 +1,299 @@
+#include "navier-stokes.h"
+#include <math.h>
+
+#define M_PI 3.14159265358979323846
+
+#define CHK 1
+
+// next time step for Irreversible Navier-Stokes equation
+int ins_step(_Complex double* u, ns_params params, fft_vects vects, _Complex double* tmp1, _Complex double* tmp2, _Complex double* tmp3){
+  int kx,ky;
+
+  // k1
+  ins_rhs(tmp1, u, params, vects);
+  // add to output
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      tmp3[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/6*tmp1[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+
+  // u+h*k1/2
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/2*tmp1[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+  // k2
+  ins_rhs(tmp1, tmp2, params, vects);
+  // add to output
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      tmp3[KLOOKUP(kx,ky,params.S)]+=params.h/3*tmp1[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+
+  // u+h*k2/2
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h/2*tmp1[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+  // k3
+  ins_rhs(tmp1, tmp2, params, vects);
+  // add to output
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      tmp3[KLOOKUP(kx,ky,params.S)]+=params.h/3*tmp1[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+
+  // u+h*k3
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      tmp2[KLOOKUP(kx,ky,params.S)]=u[KLOOKUP(kx,ky,params.S)]+params.h*tmp1[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+  // k4
+  ins_rhs(tmp1, tmp2, params, vects);
+  // add to output
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      u[KLOOKUP(kx,ky,params.S)]=tmp3[KLOOKUP(kx,ky,params.S)]+params.h/6*tmp1[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+
+  return(0);
+}
+
+// right side of Irreversible Navier-Stokes equation
+int ins_rhs(_Complex double* out, _Complex double* u, ns_params params, fft_vects vects){
+  int kx,ky;
+
+#if CHK==0
+  // F(u/|p|)*F(q1*q2*u/|q|)
+  // init to 0
+  for(kx=0; kx<params.N*params.N; kx++){
+    vects.fft1[kx]=0;
+    vects.fft2[kx]=0;
+    vects.invfft[kx]=0;
+  }
+  // fill modes
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      if(kx!=0 || ky!=0){
+        vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+        vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*ky*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+      }
+    }
+  }
+  // fft
+  fftw_execute(vects.fft1_plan);
+  fftw_execute(vects.fft2_plan);
+  // write to invfft
+  for(kx=-2*params.K;kx<=2*params.K;kx++){
+    for(ky=-2*params.K;ky<=2*params.K;ky++){
+      vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+    }
+  }
+  // inverse fft
+  fftw_execute(vects.invfft_plan);
+  
+  // write out
+  for(kx=0; kx<params.S*params.S; kx++){
+    out[kx]=0;
+  }
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      if(kx!=0 || ky!=0){
+        out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*kx-ky*ky)/params.N/params.N;
+      }
+    }
+  }
+
+  // F(u/|p|)*F((q1*q1-q2*q2)*u/|q|)
+  // init to 0
+  for(kx=0; kx<params.N*params.N; kx++){
+    vects.fft2[kx]=0;
+    vects.invfft[kx]=0;
+  }
+  // fill modes
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      if(kx!=0 || ky!=0){
+        vects.fft2[KLOOKUP(kx,ky,params.N)]=(kx*kx-ky*ky)*u[KLOOKUP(kx,ky,params.S)]/sqrt(kx*kx+ky*ky);
+      }
+    }
+  }
+  // fft
+  fftw_execute(vects.fft2_plan);
+  // write to invfft
+  for(kx=-2*params.K;kx<=2*params.K;kx++){
+    for(ky=-2*params.K;ky<=2*params.K;ky++){
+      vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+    }
+  }
+  // inverse fft
+  fftw_execute(vects.invfft_plan);
+  
+  // write out
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      if(kx!=0 || ky!=0){
+        out[KLOOKUP(kx,ky,params.S)]=out[KLOOKUP(kx,ky,params.S)]-4*M_PI*M_PI/sqrt(kx*kx+ky*ky)*vects.invfft[KLOOKUP(kx,ky,params.N)]*(kx*ky)/params.N/params.N;
+      }
+    }
+  }
+#elif CHK == 1
+  // F(-p2/|p|*u)*F(q1*|q|*u)
+  // init to 0
+  for(kx=0; kx<params.N*params.N; kx++){
+    vects.fft1[kx]=0;
+    vects.fft2[kx]=0;
+    vects.invfft[kx]=0;
+  }
+  // fill modes
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      if(kx!=0 || ky!=0){
+        vects.fft1[KLOOKUP(kx,ky,params.N)]=-kx/sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
+        vects.fft2[KLOOKUP(kx,ky,params.N)]=kx*sqrt(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)];
+      }
+    }
+  }
+  // fft
+  fftw_execute(vects.fft1_plan);
+  fftw_execute(vects.fft2_plan);
+  // write to invfft
+  for(kx=-2*params.K;kx<=2*params.K;kx++){
+    for(ky=-2*params.K;ky<=2*params.K;ky++){
+      vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)] - vects.fft1[KLOOKUP(ky,kx,params.N)]*vects.fft2[KLOOKUP(ky,kx,params.N)];
+    }
+  }
+  
+  // inverse fft
+  fftw_execute(vects.invfft_plan);
+  
+  
+  // write out
+  for(kx=0; kx<params.S*params.S; kx++){
+    out[kx]=0;
+  }
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      if(kx!=0 || ky!=0){
+        out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu/params.S*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+2*M_PI/sqrt(kx*kx+ky*ky)/params.S*vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
+      }
+    }
+  }
+
+#elif CHK==2
+  // F(u)*F(q1*u)
+  // init to 0
+  for(kx=0; kx<params.N*params.N; kx++){
+    vects.fft1[kx]=0;
+    vects.fft2[kx]=0;
+    vects.invfft[kx]=0;
+  }
+  // fill modes
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      // u_k
+      vects.fft1[KLOOKUP(kx,ky,params.N)]=u[KLOOKUP(kx,ky,params.S)];
+      // 2i\pi k_x u_k
+      vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*kx*u[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+  // fft
+  fftw_execute(vects.fft1_plan);
+  fftw_execute(vects.fft2_plan);
+  // write to invfft
+  for(kx=-2*params.K;kx<=2*params.K;kx++){
+    for(ky=-2*params.K;ky<=2*params.K;ky++){
+      vects.invfft[KLOOKUP(kx,ky,params.N)]=vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+    }
+  }
+
+  // F(p1/p2*u)*F(q2*u)
+  // init to 0
+  for(kx=0; kx<params.N*params.N; kx++){
+    vects.fft1[kx]=0;
+    vects.fft2[kx]=0;
+  }
+  // fill modes
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      // k_x/k_y u_k
+      if(ky!=0){
+	vects.fft1[KLOOKUP(kx,ky,params.N)]=kx/ky*u[KLOOKUP(kx,ky,params.S)];
+      }
+      // 2i\pi k_y u_k
+      vects.fft2[KLOOKUP(kx,ky,params.N)]=2*M_PI*I*ky*u[KLOOKUP(kx,ky,params.S)];
+    }
+  }
+  // fft
+  fftw_execute(vects.fft1_plan);
+  fftw_execute(vects.fft2_plan);
+  // write to invfft
+  for(kx=-2*params.K;kx<=2*params.K;kx++){
+    for(ky=-2*params.K;ky<=2*params.K;ky++){
+      vects.invfft[KLOOKUP(kx,ky,params.N)]+=-vects.fft1[KLOOKUP(kx,ky,params.N)]*vects.fft2[KLOOKUP(kx,ky,params.N)];
+    }
+  }
+
+  // inverse fft
+  fftw_execute(vects.invfft_plan);
+
+  /*
+  // check: convolution instead of fft
+  for(kx=0; kx<params.S*params.S; kx++){
+    out[kx]=0;
+  }
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      for(px=-params.K;px<=params.K;px++){
+	for(py=-params.K;py<=params.K;py++){
+	  if(kx-px<=params.K && kx-px>=-params.K && ky-py<=params.K && ky-py>=-params.K){
+	    out[KLOOKUP(kx,ky,params.S)]+=2*M_PI*I*(u[KLOOKUP(px,py,params.S)]*(kx-px)*u[KLOOKUP(kx-px,ky-py,params.S)]-(py==0?0:px/py*u[KLOOKUP(px,py,params.S)]*(ky-py)*u[KLOOKUP(kx-px,ky-py,params.S)]));
+	  }
+	}
+      }
+      dd=(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])*(__real__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__real__ out[KLOOKUP(kx,ky,params.S)])+(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)])*(__imag__ vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N-__imag__ out[KLOOKUP(kx,ky,params.S)]);
+      if(dd>1e-25){
+	printf("%d %d % .8e\n",kx,ky, dd);
+      }
+    }
+  }
+  */
+
+
+  // write out
+  for(kx=0; kx<params.S*params.S; kx++){
+    out[kx]=0;
+  }
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      out[KLOOKUP(kx,ky,params.S)]=-4*M_PI*M_PI*params.nu*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]+params.g[KLOOKUP(kx,ky,params.S)]+vects.invfft[KLOOKUP(kx,ky,params.N)]/params.N/params.N;
+    }
+  }
+
+#endif
+  return(0);
+}
+
+
+// compute alpha
+_Complex double compute_alpha(_Complex double* u, ns_params params){
+  _Complex double num=0;
+  _Complex double denom=0;
+  int kx,ky;
+
+  for(kx=-params.K;kx<=params.K;kx++){
+    for(ky=-params.K;ky<=params.K;ky++){
+      denom+=(kx*kx+ky*ky)*(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]*conj(u[KLOOKUP(kx,ky,params.S)])*(1+(ky!=0?kx*kx/ky/ky:0));
+      num+=(kx*kx+ky*ky)*u[KLOOKUP(kx,ky,params.S)]*conj(params.g[KLOOKUP(kx,ky,params.S)])*(1+(ky!=0?kx*kx/ky/ky:0));
+    }
+  }
+
+  return(num/denom);
+}