add L parameter

1 files changed, 15 insertions, 11 deletions

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 99e41ab..89bfbbd 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -2,8 +2,6 @@
 #include <math.h>
 #include <stdlib.h>
 
-#define M_PI 3.14159265358979323846
-
 // compute solution as a function of time
 int uk(
   int K1,
@@ -13,6 +11,7 @@ int uk(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -44,7 +43,7 @@ int uk(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
     
     if(t%print_freq==0){
       fprintf(stderr,"%d % .8e ",t,t*delta);
@@ -77,6 +76,7 @@ int energy(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -97,7 +97,7 @@ int energy(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
     
     if(t%print_freq==0){
 
@@ -126,6 +126,7 @@ int enstrophy(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int print_freq,
   unsigned int nthreads
@@ -150,7 +151,7 @@ int enstrophy(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
     alpha=compute_alpha(u, K1, K2, g);
     
     // running average
@@ -177,6 +178,7 @@ int quiet(
   unsigned int nsteps,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   unsigned int nthreads
 ){
@@ -194,7 +196,7 @@ int quiet(
 
   // iterate
   for(t=0;t<nsteps;t++){
-    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+    ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
   }
 
   ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
@@ -342,6 +344,7 @@ int ins_step(
   int N2,
   double nu,
   double delta,
+  double L,
   _Complex double (*g)(int,int),
   fft_vect fft1,
   fft_vect fft2,
@@ -353,7 +356,7 @@ int ins_step(
   int kx,ky;
 
   // k1
-  ins_rhs(tmp1, u, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, u, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -368,7 +371,7 @@ int ins_step(
     }
   }
   // k2
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -383,7 +386,7 @@ int ins_step(
     }
   }
   // k3
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -398,7 +401,7 @@ int ins_step(
     }
   }
   // k4
-  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, g, fft1, fft2, ifft);
+  ins_rhs(tmp1, tmp2, K1, K2, N1, N2, nu, L, g, fft1, fft2, ifft);
   // add to output
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -418,6 +421,7 @@ int ins_rhs(
   int N1,
   int N2,
   double nu,
+  double L,
   _Complex double (*g)(int,int),
   fft_vect fft1,
   fft_vect fft2,
@@ -447,7 +451,7 @@ int ins_rhs(
   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)];
+        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)(kx,ky)+4*M_PI*M_PI/L/L/sqrt(kx*kx+ky*ky)*ifft.fft[klookup(kx,ky,N1,N2)];
       }
     }
   }