parameters to choose rk2

4 files changed, 73 insertions, 26 deletions

diff --git a/src/constants.cpp b/src/constants.cpp
new file mode 100644
index 0000000..03daea7
--- /dev/null
+++ b/src/constants.cpp
@@ -0,0 +1,17 @@
+#define COMMAND_UK 1
+#define COMMAND_EEA 2
+#define COMMAND_QUIET 3
+
+#define DRIVING_ZERO 1
+#define DRIVING_TEST 2
+#define DRIVING_FILE 3
+#define DRIVING_FILE_TXT 4
+
+#define INIT_RANDOM 1
+#define INIT_GAUSSIAN 2
+#define INIT_FILE 3
+#define INIT_FILE_TXT 4
+
+#define ALGORITHM_RK4 4
+#define ALGORITHM_RK2 2
+
diff --git a/src/main.c b/src/main.c
index d7d30bb..4ad182e 100644
--- a/src/main.c
+++ b/src/main.c
@@ -9,12 +9,14 @@
 #include <stdbool.h>
 #include <stdint.h>
 #include <errno.h>
-#include "navier-stokes.h"
+
+#include "constants.cpp"
 #include "driving.h"
 #include "init.h"
 #include "int_tools.h"
-
 #include "io.h"
+#include "navier-stokes.h"
+
 
 // structure to store parameters, to make it easier to pass parameters to CLI functions
 typedef struct nstrophy_parameters {
@@ -33,6 +35,7 @@ typedef struct nstrophy_parameters {
   uint64_t starting_time;
   unsigned int driving;
   unsigned int init;
+  unsigned int algorithm;
   FILE* initfile;
   FILE* drivingfile;
 } nstrophy_parameters;
@@ -55,21 +58,6 @@ _Complex double* set_init(nstrophy_parameters parameters);
 // signal handler
 void sig_handler (int signo);
 
-#define COMMAND_UK 1
-#define COMMAND_EEA 2
-#define COMMAND_QUIET 3
-
-#define DRIVING_ZERO 1
-#define DRIVING_TEST 2
-#define DRIVING_FILE 3
-#define DRIVING_FILE_TXT 4
-
-#define INIT_RANDOM 1
-#define INIT_GAUSSIAN 2
-#define INIT_FILE 3
-#define INIT_FILE_TXT 4
-
-
 // global variable to handle interrupts
 volatile bool g_abort = false;
 // signal handler
@@ -165,15 +153,15 @@ int main (
 
   // run command
   if (command==COMMAND_UK){
-    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile);
+    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile);
   }
   else if(command==COMMAND_EEA){
     // register signal handler to handle aborts
     signal(SIGINT, sig_handler);
-    eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.print_freq, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_str);
+    eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_str);
   }
   else if(command==COMMAND_QUIET){
-    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, parameters.starting_time, u0, g, parameters.irreversible, nthreads, savefile);
+    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, parameters.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile);
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
@@ -250,6 +238,18 @@ int print_params(
     break;
   }
 
+  switch(parameters.algorithm){
+  case ALGORITHM_RK4:
+    fprintf(file,", algorithm=RK4");
+    break;
+  case ALGORITHM_RK2:
+    fprintf(file,", algorithm=RK2");
+    break;
+  default:
+    fprintf(file,", algorithm=RK4");
+    break;
+  }
+
   fprintf(file,"\n");
 
   return 0;
@@ -376,6 +376,7 @@ int read_params(
   parameters->drivingfile=NULL;
   parameters->init=INIT_GAUSSIAN;
   parameters->initfile=NULL;
+  parameters->algorithm=ALGORITHM_RK4;
 
   if (param_str!=NULL){
     // init
@@ -623,6 +624,19 @@ int set_parameter(
       return(-1);
     }
   }
+  // algorithm
+  else if (strcmp(lhs,"algorithm")==0){
+    if (strcmp(rhs,"RK4")==0){
+      parameters->algorithm=ALGORITHM_RK4;
+    }
+    else if (strcmp(rhs,"RK2")==0){
+      parameters->algorithm=ALGORITHM_RK2;
+    }
+    else{
+      fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs);
+      return(-1);
+    }
+  }
   else{
     fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs);
     return(-1);
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 755e829..6387fcf 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -1,3 +1,4 @@
+#include "constants.cpp"
 #include "io.h"
 #include "navier-stokes.h"
 #include "statistics.h"
@@ -18,6 +19,7 @@ int uk(
   _Complex double* u0,
   _Complex double* g,
   bool irreversible,
+  unsigned int algorithm,
   uint64_t print_freq,
   uint64_t starting_time,
   unsigned int nthreads,
@@ -51,7 +53,11 @@ int uk(
 
   // iterate
   for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
-    ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+    if(algorithm==ALGORITHM_RK2){
+      ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
+    } else {
+      ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+    }
     
     if(t%print_freq==0){
       fprintf(stderr,"%lu % .8e ",t,t*delta);
@@ -91,6 +97,7 @@ int eea(
   _Complex double* u0,
   _Complex double* g,
   bool irreversible,
+  unsigned int algorithm,
   uint64_t print_freq,
   uint64_t starting_time,
   unsigned int nthreads,
@@ -127,7 +134,11 @@ int eea(
 
   // iterate
   for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
-    ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+    if(algorithm==ALGORITHM_RK2){
+      ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
+    } else {
+      ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+    }
 
     energy=compute_energy(u, K1, K2);
     alpha=compute_alpha(u, K1, K2, g, L);
@@ -196,6 +207,7 @@ int quiet(
   _Complex double* u0,
   _Complex double* g,
   bool irreversible,
+  unsigned int algorithm,
   unsigned int nthreads,
   FILE* savefile
 ){
@@ -214,7 +226,11 @@ int quiet(
 
   // iterate
   for(t=starting_time;nsteps==0 || t<starting_time+nsteps;t++){
-    ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+    if(algorithm==ALGORITHM_RK2){
+      ns_step_rk2(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, irreversible);
+    } else {
+      ns_step_rk4(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3, irreversible);
+    }
   }
 
   // save final entry to savefile
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 77b4cdb..933fe1d 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -18,13 +18,13 @@ typedef struct fft_vects {
 } fft_vect;
 
 // compute u_k
-int uk( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, uint64_t print_freq, uint64_t starting_time, unsigned int nthreadsl, FILE* savefile);
+int uk( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_time, unsigned int nthreadsl, FILE* savefile);
 
 // compute energy, enstrophy and alpha
-int eea( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, uint64_t print_freq, uint64_t starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_string);
+int eea( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_string);
 
 // 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, uint64_t nsteps, double nu, double delta, double L, uint64_t starting_time, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int nthreads, FILE* savefile);
+int quiet( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, uint64_t starting_time, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile);
 
 
 // initialize vectors for computation