From 4b7c89fb354790064d4d823e3b657e43c38d6ce8 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian@jauslin.org>
Date: Wed, 25 May 2022 10:23:53 -0400
Subject: Same initial u as Gallavotti

---
 src/navier-stokes.c | 34 +++++++++++++++-------------------
 1 file changed, 15 insertions(+), 19 deletions(-)

(limited to 'src')

diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 0a2842c..99e41ab 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -278,31 +278,26 @@ int ns_init_u(
   int K2
 ){
   int kx,ky;
-  /*
-  double rescale;
 
   srand(17);
 
   // random init (set half, then the other half are the conjugates)
-  for(ky=0;ky<=K2;ky++){
-    u[klookup(0,ky,2*K1+1,2*K2+1)]=(-RAND_MAX*0.5+rand())*1.0/RAND_MAX+(-RAND_MAX*0.5+rand())*1.0/RAND_MAX*I;
-  }
-  for(kx=1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=(-RAND_MAX*0.5+rand())*1.0/RAND_MAX+(-RAND_MAX*0.5+rand())*1.0/RAND_MAX*I;
-    }
-  }
-  // conjugates
-  for(ky=-K2;ky<=-1;ky++){
-    u[klookup(0,ky,2*K1+1,2*K2+1)]=conj(u[klookup(0,-ky,2*K1+1,2*K2+1)]);
-  }
-  for(kx=-K1;kx<=-1;kx++){
+  for(kx=0;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=conj(u[klookup(-kx,-ky,2*K1+1,2*K2+1)]);
+      if (kx==0 && ky<=0){
+	u[klookup(kx,ky,2*K1+1,2*K2+1)]=0.;
+      }
+      else{
+	double x=-0.5+((double) rand())/RAND_MAX;
+	double y=-0.5+((double) rand())/RAND_MAX;
+	u[klookup(kx,ky,2*K1+1,2*K2+1)]=x+y*I;
+	u[klookup(-kx,-ky,2*K1+1,2*K2+1)]=conj(u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+      }
     }
   }
 
-  // rescale: 1/k decay
+  // rescale to match with Gallavotti's initialization
+  double rescale;
   rescale=0;
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -311,10 +306,9 @@ int ns_init_u(
   }
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      u[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(155.1/rescale);
+      u[klookup(kx,ky,2*K1+1,2*K2+1)]=u[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(1.54511597324389e+02/rescale);
     }
   }
-  */
 
 
   /*
@@ -326,12 +320,14 @@ int ns_init_u(
   }
   */
 
+  /*
   // exponentially decaying init
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
       u[klookup(kx,ky,2*K1+1,2*K2+1)]=exp(-sqrt(kx*kx+ky*ky));
     }
   }
+  */
 
   return 0;
 }
-- 
cgit v1.2.3-54-g00ecf