Replace init_en with init_enstrophy and init_energy

README.md

@@ -102,8 +102,13 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are
   force, excluding kx<0 and kx=0&&ky<=0. The plaintext file format is
   `kx ky real_part imag_part`.
 
-* `init_en` (double, default 1.54511597324389e+02): initial value of the energy if
+* `init_energy` (double, default is to not rescale the initial condition, is
-  `equation=irreversible` or of the enstrophy if `equation=reversible`.
+  incompatible with `init_enstrophy`): enforce a value for the energy of the
+  initial condition.
+
+* `init_enstrophy` (double, default is to not rescale the initial condition, is
+  incompatible with `init_energy`): enforce a value for the enstrophy of the
+  initial condition.
 
 * `random_seed` (int, default ): seed for random initialization.
 

src/constants.cpp

@@ -45,3 +45,6 @@ limitations under the License.
 #define COST_k32 3
 #define COST_ENSTROPHY 4
 #define COST_ALPHA 5
+
+#define FIX_ENSTROPHY 1
+#define FIX_ENERGY 2

src/init.c

@@ -17,21 +17,17 @@ limitations under the License.
 #include "init.h"
 #include "navier-stokes.h"
 #include "io.h"
-#include <stdlib.h>
 #include <math.h>
+#include <stdlib.h>
 
 // random initial condition
 int init_random (
   _Complex double* u0,
-  double init_en,
   int K1,
   int K2,
-  double L,
+  int seed
-  int seed,
-  bool irreversible
 ){
   int kx,ky;
-  double rescale;
   double x,y;
 
   srand(seed);
@@ -45,33 +41,16 @@ int init_random (
     }
   }
 
-  if (irreversible) {
-    // fix the energy
-    rescale=compute_energy(u0, K1, K2);
-  } else {
-    // fix the enstrophy
-    rescale=compute_enstrophy(u0, K1, K2, L);
-  }
-  for(kx=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
-    }
-  }
-
   return 0;
 }
 
 // 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=0;kx<=K1;kx++){
     for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
@@ -79,19 +58,6 @@ 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=0;kx<=K1;kx++){
-    for(ky=(kx>0 ? -K2 : 1);ky<=K2;ky++){
-      u0[klookup_sym(kx,ky,K2)]=u0[klookup_sym(kx,ky,K2)]*sqrt(init_en/rescale);
-    }
-  }
-
   return 0;
 }
 

src/init.h

@@ -21,10 +21,10 @@ limitations under the License.
 #include <stdbool.h>
 
 // random initial condition
-int init_random(_Complex double* u0, double init_en, int K1, int K2, double L, int seed, bool irreversible);
+int init_random(_Complex double* u0, int K1, int K2, int seed);
 
 // Gaussian initial condition
-int init_gaussian(_Complex double* u0, double init_en, int K1, int K2, double L, bool irreversible);
+int init_gaussian(_Complex double* u0, int K1, int K2);
 
 // Initialize from file
 int init_file_txt (_Complex double* u0, int K1, int K2, FILE* initfile);

src/main.c

@@ -16,6 +16,7 @@ limitations under the License.
 
 #define VERSION "0.1"
 
+#include <math.h>
 #include <signal.h>
 #include <string.h>
 #include <stdlib.h>
@@ -33,7 +34,9 @@ limitations under the License.
 
 // 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
+  double init_energy;
+  double init_enstrophy;
+  unsigned int init_enstrophy_or_energy; //whether to fix the initial enstrophy or energy
   bool irreversible;
   int K1;
   int K2;
@@ -221,6 +224,15 @@ int main (
   // set initial condition
   u0=set_init(&parameters);
 
+  // if init_enstrophy is not set in the parameters, then compute it from the initial condition
+  if (parameters.init_enstrophy_or_energy!=FIX_ENSTROPHY){
+    parameters.init_enstrophy=compute_enstrophy(u0, parameters.K1, parameters.K2, parameters.L);
+  }
+  // same with init_energy
+  if (parameters.init_enstrophy_or_energy!=FIX_ENERGY){
+    parameters.init_energy=compute_energy(u0, parameters.K1, parameters.K2);
+  }
+
   // close initfile (do this early, so that it is possible to use the same file for init and save)
   if (parameters.initfile!=NULL){
     fclose(parameters.initfile);
@@ -272,19 +284,19 @@ int main (
 
   // run command
   if (command==COMMAND_UK){
-    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
+    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
   }
   else if(command==COMMAND_ENSTROPHY){
     // register signal handler to handle aborts
     signal(SIGINT, sig_handler);
     signal(SIGTERM, sig_handler);
-    enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.print_freq, parameters.starting_time, parameters.print_alpha, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
+    enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, parameters.print_alpha, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_str), dstring_to_str_noinit(&savefile_str), dstring_to_str_noinit(&utfile_str));
   }
   else if(command==COMMAND_QUIET){
-    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, nthreads, savefile);
+    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, parameters.starting_time, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, nthreads, savefile);
   }
   else if(command==COMMAND_LYAPUNOV){
-    lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.nu, parameters.D_epsilon, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_en, parameters.algorithm, parameters.starting_time, nthreads);
+    lyapunov(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.final_time, parameters.lyapunov_reset, parameters.nu, parameters.D_epsilon, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, parameters.max_delta, parameters.adaptive_cost, u0, g, parameters.irreversible, parameters.keep_en_cst, parameters.init_enstrophy, parameters.algorithm, parameters.starting_time, nthreads);
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
@@ -333,7 +345,12 @@ int print_params(
     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);
+  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);
+  if (parameters.init_enstrophy_or_energy==FIX_ENSTROPHY){
+    fprintf(file,", init_enstrophy=%.15e", parameters.init_enstrophy);
+  } else if (parameters.init_enstrophy_or_energy==FIX_ENERGY){
+    fprintf(file,", init_energy=%.15e", parameters.init_energy);
+  }
 
   switch(parameters.driving){
   case DRIVING_TEST:
@@ -535,7 +552,8 @@ int read_args(
 int set_default_params(
   nstrophy_parameters* parameters
 ){
-  parameters->init_en=1.54511597324389e+02;
+  // no choice
+  parameters->init_enstrophy_or_energy=0; // do not set init_enstrophy or init_energy here: do that after reading init
   parameters->irreversible=true;
   parameters->K1=16;
   parameters->K2=parameters->K1;
@@ -658,10 +676,29 @@ int set_parameter(
       return(-1);
     }
   }
-  else if (strcmp(lhs,"init_en")==0){
+  else if (strcmp(lhs,"init_enstrophy")==0){
-    ret=sscanf(rhs,"%lf",&(parameters->init_en));
+    // check that the energy was not also set
+    if (parameters->init_enstrophy_or_energy==FIX_ENERGY){
+      fprintf(stderr, "error: cannot use init_enstrophy if init_energy has been used\n");
+      return(-1);
+    }
+    parameters->init_enstrophy_or_energy=FIX_ENSTROPHY;
+    ret=sscanf(rhs,"%lf",&(parameters->init_enstrophy));
     if(ret!=1){
-      fprintf(stderr, "error: parameter 'init_en' should be a double\n       got '%s'\n",rhs);
+      fprintf(stderr, "error: parameter 'init_enstrophy' should be a double\n       got '%s'\n",rhs);
+      return(-1);
+    }
+  }
+  else if (strcmp(lhs,"init_energy")==0){
+    // check that the enstrophy was not also set
+    if (parameters->init_enstrophy_or_energy==FIX_ENSTROPHY){
+      fprintf(stderr, "error: cannot use init_energy if init_enstrophy has been used\n");
+      return(-1);
+    }
+    parameters->init_enstrophy_or_energy=FIX_ENERGY;
+    ret=sscanf(rhs,"%lf",&(parameters->init_energy));
+    if(ret!=1){
+      fprintf(stderr, "error: parameter 'init_energy' should be a double\n       got '%s'\n",rhs);
       return(-1);
     }
   }
@@ -1009,11 +1046,11 @@ _Complex double* set_init(
 
   switch(parameters->init){
   case INIT_RANDOM:
-    init_random(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->seed, parameters->irreversible);
+    init_random(u0, parameters->K1, parameters->K2, parameters->seed);
     break;
 
   case INIT_GAUSSIAN:
-    init_gaussian(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->irreversible);
+    init_gaussian(u0, parameters->K1, parameters->K2);
     break;
 
   case INIT_FILE:
@@ -1031,9 +1068,30 @@ _Complex double* set_init(
     break;
 
   default:
-    init_gaussian(u0, parameters->init_en, parameters->K1, parameters->K2, parameters->L, parameters->irreversible);
+    init_gaussian(u0, parameters->K1, parameters->K2);
     break;
   }
 
+  // rescale energy or enstrophy
+  if (parameters->init_enstrophy_or_energy==FIX_ENERGY) {
+    // fix the energy
+    double rescale=compute_energy(u0, parameters->K1, parameters->K2);
+    int kx,ky;
+    for(kx=0;kx<=parameters->K1;kx++){
+      for(ky=(kx>0 ? -parameters->K2 : 1);ky<=parameters->K2;ky++){
+	u0[klookup_sym(kx,ky,parameters->K2)]=u0[klookup_sym(kx,ky,parameters->K2)]*sqrt(parameters->init_energy/rescale);
+      }
+    }
+  } else if (parameters->init_enstrophy_or_energy==FIX_ENSTROPHY) {
+    // fix the enstrophy
+    double rescale=compute_enstrophy(u0, parameters->K1, parameters->K2, parameters->L);
+    int kx,ky;
+    for(kx=0;kx<=parameters->K1;kx++){
+      for(ky=(kx>0 ? -parameters->K2 : 1);ky<=parameters->K2;ky++){
+	u0[klookup_sym(kx,ky,parameters->K2)]=u0[klookup_sym(kx,ky,parameters->K2)]*sqrt(parameters->init_enstrophy/rescale);
+      }
+    }
+  }
+  
   return u0;
 }