From 199b8f0df5adeaaac4ca2afc2eff5237dfee36c3 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian@jauslin.org>
Date: Wed, 18 May 2022 22:22:42 +0200
Subject: Quiet mode

---
 src/main.c          |  7 +++++
 src/navier-stokes.c | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 src/navier-stokes.h |  7 +++++
 3 files changed, 97 insertions(+)

(limited to 'src')

diff --git a/src/main.c b/src/main.c
index 3383f10..904b5bf 100644
--- a/src/main.c
+++ b/src/main.c
@@ -20,6 +20,7 @@ int set_parameter(char* lhs, char* rhs, int* K1, int* K2, int* N1, int* N2, unsi
 
 #define COMMAND_UK 1
 #define COMMAND_ENSTROPHY 2
+#define COMMAND_QUIET 3
 
 #define DRIVING_TEST 1
 
@@ -69,6 +70,9 @@ int main (
   else if(command==COMMAND_ENSTROPHY){
     enstrophy(K1, K2, N1, N2, nsteps, nu, delta, g, print_freq);
   }
+  else if(command==COMMAND_QUIET){
+    quiet(K1, K2, N1, N2, nsteps, nu, delta, g);
+  }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
     print_usage();
@@ -145,6 +149,9 @@ int read_args(
       else if(strcmp(argv[i], "enstrophy")==0){
 	*command=COMMAND_ENSTROPHY;
       }
+      else if(strcmp(argv[i], "quiet")==0){
+	*command=COMMAND_QUIET;
+      }
       else{
 	fprintf(stderr, "error: unrecognized command: '%s'\n",argv[i]);
 	return(-1);
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 13d9531..87f5aba 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -105,6 +105,38 @@ int enstrophy(
   return(0);
 }
 
+// compute solution as a function of time, but do not print anything (useful for debugging)
+int quiet(
+  int K1,
+  int K2,
+  int N1,
+  int N2,
+  unsigned int nsteps,
+  double nu,
+  double delta,
+  _Complex double (*g)(int,int)
+){
+  _Complex double* u;
+  _Complex double* tmp1;
+  _Complex double* tmp2;
+  _Complex double* tmp3;
+  unsigned int t;
+  fft_vect fft1;
+  fft_vect fft2;
+  fft_vect ifft;
+
+  ns_init_tmps(&u, &tmp1, &tmp2, &tmp3, &fft1, &fft2, &ifft, K1, K2, N1, N2);
+  ns_init_u(u, K1, K2);
+
+  // iterate
+  for(t=0;t<nsteps;t++){
+    ins_step(u, K1, K2, N1, N2, nu, delta, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
+  }
+
+  ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
+  return(0);
+}
+
 
 // initialize vectors for computation
 int  ns_init_tmps(
@@ -332,6 +364,18 @@ int ins_rhs(
   // compute convolution term
   ns_T(u,K1,K2,N1,N2,fft1,fft2,ifft);
 
+
+  /*
+  // compare convolution term (store result in fft1.fft)
+  ns_T_nofft(fft1.fft, u, K1, K2, N1, N2);
+  double cmp=0.;
+  for(i=0;i<N1*N2;i++){
+    cmp+=(ifft.fft[i]-fft1.fft[i])*(ifft.fft[i]-fft1.fft[i]);
+  }
+  printf("% .15e\n",sqrt(cmp));
+  */
+
+
   for(i=0; i<(2*K1+1)*(2*K2+1); i++){
     out[i]=0;
   }
@@ -415,6 +459,45 @@ int ns_T(
   return(0);
 }
 
+// convolution term in right side of convolution equation, computed without fourier transform
+int ns_T_nofft(
+  _Complex double* out,
+  _Complex double* u,
+  int K1,
+  int K2,
+  int N1,
+  int N2
+){
+  int kx,ky;
+  int px,py;
+  int qx,qy;
+
+  // loop over K's (needs N1>=4*K1+1 and N2>=4*K2+1)
+  if (N1<4*K1+1 || N2<4*K2+1){
+    fprintf(stderr,"error: N1 and N2 need t be >= 4*K1+1 and 4*K2+1 respectively\n");
+    return(-1);
+  }
+  for(kx=-2*K1;kx<=2*K1;kx++){
+    for(ky=-2*K2;ky<=2*K2;ky++){
+      // init
+      out[klookup(kx,ky,N1,N2)]=0.;
+
+      for(px=-K1;px<=K1;px++){
+	for(py=-K2;py<=K2;py++){
+	  qx=kx-px;
+	  qy=ky-py;
+
+	  // cutoff in q
+	  if(qx>=-K1 && qx<=K1 && qy>=-K2 && qy<=K2 && qx*qx+qy*qy>0 && px*px+py*py>0){
+	    out[klookup(kx,ky,N1,N2)]+=(-qx*py+qy*px)*sqrt(qx*qx+qy*qy)/sqrt(px*px+py*py)*u[klookup(px,py,2*K1+1,2*K2+1)]*u[klookup(qx,qy,2*K1+1,2*K2+1)];
+	  }
+	}
+      }
+    }
+  }
+
+  return 0;
+}
 
 // compute alpha
 _Complex double compute_alpha(
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 2ff4e3f..d144fce 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -17,6 +17,10 @@ int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double d
 // compute enstrophy
 int enstrophy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int), unsigned int print_freq);
 
+// compute solution as a function of time, but do not print anything (useful for debugging)
+int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, _Complex double (*g)(int,int));
+
+
 // initialize vectors for computation
 int  ns_init_tmps( _Complex double **u, _Complex double ** tmp1, _Complex double **tmp2, _Complex double **tmp3, fft_vect* fft1, fft_vect *fft2, fft_vect *ifft, int K1, int K2, int N1, int N2);
 // release vectors
@@ -34,6 +38,9 @@ int ins_rhs( _Complex double* out, _Complex double* u, int K1, int K2, int N1, i
 // 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);
 
+// convolution term in right side of equation
+int ns_T_nofft( _Complex double* out, _Complex double* u, int K1, int K2, int N1, int N2);
+
 // compute alpha
 _Complex double compute_alpha( _Complex double* u, int K1, int K2, _Complex double (*g)(int,int));
 
-- 
cgit v1.2.3-54-g00ecf