From 1616b6bbae0b30f8c2e94283f613df1adf3cbf95 Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian.jauslin@rutgers.edu>
Date: Wed, 5 Apr 2023 22:41:19 -0400
Subject: Fix initial enstrophy

---
 src/init.c | 38 +++++++++++++++++++++++++++++---------
 1 file changed, 29 insertions(+), 9 deletions(-)

(limited to 'src/init.c')

diff --git a/src/init.c b/src/init.c
index 5f3987b..4865ed9 100644
--- a/src/init.c
+++ b/src/init.c
@@ -7,9 +7,12 @@
 // random initial condition
 int init_random (
   _Complex double* u0,
+  double init_en,
   int K1,
   int K2,
-  int seed
+  double L,
+  int seed,
+  bool irreversible
 ){
   int kx,ky;
   double rescale;
@@ -29,16 +32,16 @@ int init_random (
     }
   }
 
-  // rescale to match with Gallavotti's initialization
-  rescale=0;
-  for(kx=-K1;kx<=K1;kx++){
-    for(ky=-K2;ky<=K2;ky++){
-      rescale=rescale+((__real__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])*(__real__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])+(__imag__ u0[klookup(kx,ky,2*K1+1,2*K2+1)])*(__imag__ u0[klookup(kx,ky,2*K1+1,2*K2+1)]))*(kx*kx+ky*ky);
-    }
+  if (irreversible) {
+    // fix the energy
+    rescale=compute_energy(u0, K1, K2);
+  } else {
+    // fix the enstrophy
+    rescale=compute_enstrophy(u0, K1, K2, L);
   }
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
-      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(1.54511597324389e+02/rescale);
+      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(init_en/rescale);
     }
   }
 
@@ -48,10 +51,14 @@ int init_random (
 // Gaussian initial condition
 int init_gaussian (
   _Complex double* u0,
+  double init_en,
   int K1,
-  int K2
+  int K2,
+  double L,
+  bool irreversible
 ){
   int kx,ky;
+  double rescale;
 
   for(kx=-K1;kx<=K1;kx++){
     for(ky=-K2;ky<=K2;ky++){
@@ -59,6 +66,19 @@ int init_gaussian (
     }
   }
 
+  if (irreversible) {
+    // fix the energy
+    rescale=compute_energy(u0, K1, K2);
+  } else {
+    // fix the enstrophy
+    rescale=compute_enstrophy(u0, K1, K2, L);
+  }
+  for(kx=-K1;kx<=K1;kx++){
+    for(ky=-K2;ky<=K2;ky++){
+      u0[klookup(kx,ky,2*K1+1,2*K2+1)]=u0[klookup(kx,ky,2*K1+1,2*K2+1)]*sqrt(init_en/rescale);
+    }
+  }
+
   return 0;
 }
 
-- 
cgit v1.2.3-54-g00ecf