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/navier-stokes.c | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 83 insertions(+)

(limited to 'src/navier-stokes.c')

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(
-- 
cgit v1.2.3-54-g00ecf