Fix initial enstrophy

3 files changed, 53 insertions, 15 deletions

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;
 }
 
diff --git a/src/init.h b/src/init.h
index caa7e2e..ef59677 100644
--- a/src/init.h
+++ b/src/init.h
@@ -2,12 +2,13 @@
 #define INIT_H
 
 #include <stdio.h>
+#include <stdbool.h>
 
 // random initial condition
-int init_random(_Complex double* u0, int K1, int K2, int seed);
+int init_random(_Complex double* u0, double init_en, int K1, int K2, double L, int seed, bool irreversible);
 
 // Gaussian initial condition
-int init_gaussian(_Complex double* u0, int K1, int K2);
+int init_gaussian(_Complex double* u0, double init_en, int K1, int K2, double L, bool irreversible);
 
 // Initialize from file
 int init_file (_Complex double* u0, int K1, int K2, FILE* initfile);
diff --git a/src/main.c b/src/main.c
index f60a38e..9962c77 100644
--- a/src/main.c
+++ b/src/main.c
@@ -14,6 +14,7 @@
 
 // structure to store parameters, to make it easier to pass parameters to CLI functions
 typedef struct nstrophy_parameters {
+  double init_en; // initial value for the energy for ins and enstrophy for rns
   bool irreversible;
   int K1;
   int K2;
@@ -159,7 +160,15 @@ int print_params(
   char* initfile_str,
   FILE* file
 ){
-  fprintf(file,"# K1=%d, K2=%d, N1=%d, N2=%d, nu=%.15e, delta=%.15e, L=%.15e", parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nu, parameters.delta, parameters.L);
+  fprintf(file,"# ");
+
+  if (parameters.irreversible){
+    fprintf(file,"equation=irreversible");
+  } else {
+    fprintf(file,"equation=reversible");
+  }
+
+  fprintf(file,", K1=%d, K2=%d, N1=%d, N2=%d, nu=%.15e, delta=%.15e, L=%.15e, init_en=%.15e", parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nu, parameters.delta, parameters.L, parameters.init_en);
 
   switch(driving){
   case DRIVING_TEST:
@@ -338,6 +347,7 @@ int read_params(
   bool lhs=true;
 
   // defaults
+  parameters->init_en=1.54511597324389e+02;
   parameters->irreversible=true;
   parameters->K1=16;
   parameters->K2=parameters->K1;
@@ -442,6 +452,13 @@ int set_parameter(
       return(-1);
     }
   }
+  else if (strcmp(lhs,"init_en")==0){
+    ret=sscanf(rhs,"%lf",&(parameters->init_en));
+    if(ret!=1){
+      fprintf(stderr, "error: parameter 'init_en' should be a double\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
   else if (strcmp(lhs,"K1")==0){
     ret=sscanf(rhs,"%d",&(parameters->K1));
     if(ret!=1){
@@ -579,11 +596,11 @@ _Complex double* set_init(
 
   switch(init){
   case INIT_RANDOM:
-    init_random(u0, parameters.K1, parameters.K2, parameters.seed);
+    init_random(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.seed, parameters.irreversible);
     break;
 
   case INIT_GAUSSIAN:
-    init_gaussian(u0, parameters.K1, parameters.K2);
+    init_gaussian(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.irreversible);
     break;
 
   case INIT_FILE:
@@ -591,7 +608,7 @@ _Complex double* set_init(
     break;
 
   default:
-    init_gaussian(u0, parameters.K1, parameters.K2);
+    init_gaussian(u0, parameters.init_en, parameters.K1, parameters.K2, parameters.L, parameters.irreversible);
     break;
   }