savefile and initfile

7 files changed, 231 insertions, 21 deletions

diff --git a/src/init.c b/src/init.c
index 9711db3..5f3987b 100644
--- a/src/init.c
+++ b/src/init.c
@@ -1,5 +1,6 @@
 #include "init.h"
 #include "navier-stokes.h"
+#include "io.h"
 #include <stdlib.h>
 #include <math.h>
 
@@ -60,3 +61,14 @@ int init_gaussian (
 
   return 0;
 }
+
+// Initialize from file
+int init_file (
+  _Complex double* u0,
+  int K1,
+  int K2,
+  FILE* initfile
+){
+  read_u(u0, K1, K2, initfile);
+  return 0;
+}
diff --git a/src/init.h b/src/init.h
index 8631b38..caa7e2e 100644
--- a/src/init.h
+++ b/src/init.h
@@ -1,10 +1,15 @@
 #ifndef INIT_H
 #define INIT_H
 
+#include <stdio.h>
+
 // random initial condition
 int init_random(_Complex double* u0, int K1, int K2, int seed);
 
 // Gaussian initial condition
 int init_gaussian(_Complex double* u0, int K1, int K2);
 
+// Initialize from file
+int init_file (_Complex double* u0, int K1, int K2, FILE* initfile);
+
 #endif
diff --git a/src/io.c b/src/io.c
new file mode 100644
index 0000000..5cdfc54
--- /dev/null
+++ b/src/io.c
@@ -0,0 +1,112 @@
+#include <stdbool.h>
+#include <stdlib.h>
+#include "io.h"
+#include "navier-stokes.h"
+
+// write final entry to file
+int write_u(_Complex double* u, int K1, int K2, FILE* file){
+  int kx,ky;
+
+  // do nothing if there is no file
+  if(file==NULL){
+    return 0;
+  }
+
+  for(kx=-K1;kx<=K1;kx++){
+    for (ky=-K2;ky<=K2;ky++){
+      fprintf(file,"%d %d % .15e % .15e\n",kx,ky,__real__ u[klookup(kx,ky,2*K1+1,2*K2+1)],__imag__ u[klookup(kx,ky,2*K1+1,2*K2+1)]);
+    }
+  }
+
+  return 0;
+}
+
+// read u from file
+int read_u(_Complex double* u, int K1, int K2, FILE* file){
+  int kx,ky;
+  double r,i;
+  char* line;
+  unsigned int len=256;
+  unsigned int pos=0;
+  char* line_realloc;
+  char c;
+  int ret;
+  unsigned int counter=0;
+  bool commented=false;
+
+  // error if there is no file (this should not happen)
+  if (file==NULL){
+    fprintf(stderr,"error reading u from file (this is a bug!)\n");
+    return -1;
+  }
+
+  // allocate line buffer
+  line=calloc(sizeof(char), len);
+
+  while(1){
+    c=fgetc(file);
+
+    // end of file
+    if (feof(file)){
+      break;
+    }
+
+    // newline: read line and reset buffer
+    if(c=='\n' || c=='\r'){
+      // increment line counter
+      counter++;
+
+      // read entry
+      // ignore empty lines
+      if(pos>0){
+	ret=sscanf(line, "%d %d %le %le", &kx, &ky, &r, &i);
+
+	// errors
+	if(ret!=4){
+	  fprintf(stderr, "warning: line %d does not match the input format: '%s'\n", counter, line);
+	}
+	else{
+	  if(kx>K1 || kx<-K1 || ky>K2 || ky<-K2){
+	    fprintf(stderr, "warning: reading line %d: kx or ky out of bounds: %d,%d\n", counter, kx, ky);
+	  }
+	  else{
+	    // set u
+	    u[klookup(kx, ky, 2*K1+1, 2*K2+1)]=r+i*I;
+	  }
+	}
+      }
+
+      // reset buffer
+      pos=0;
+      line[pos]='\0';
+      commented=false;
+    }
+    // comment: stop reading
+    else if (c=='#'){
+      commented=true;
+    }
+    // add to buffer (unless we are in a comment)
+    else if (!(commented)){
+      // check that there is room in buffer
+      if(pos==len){
+	// too short: reallocate
+	line_realloc=calloc(sizeof(char), 2*len);
+	for(pos=0;pos<len;pos++){
+	  line_realloc[pos]=line[pos];
+	}
+	free(line);
+	line=line_realloc;
+	len=2*len;
+      }
+
+      // add c to line
+      line[pos]=c;
+      pos++;
+      line[pos]='\0';
+    }
+  }
+
+  free(line);
+
+  return 0;
+}
diff --git a/src/io.h b/src/io.h
new file mode 100644
index 0000000..6a35707
--- /dev/null
+++ b/src/io.h
@@ -0,0 +1,12 @@
+#ifndef IO_H
+#define IO_H
+
+#include <stdio.h>
+
+// write u to file
+int write_u(_Complex double* u, int K1, int K2, FILE* file);
+
+// read u from file
+int read_u(_Complex double* u, int K1, int K2, FILE* file);
+
+#endif
diff --git a/src/main.c b/src/main.c
index e4152f0..92c53b9 100644
--- a/src/main.c
+++ b/src/main.c
@@ -6,6 +6,7 @@
 #include <string.h>
 #include <stdlib.h>
 #include <stdbool.h>
+#include <errno.h>
 #include "navier-stokes.h"
 #include "driving.h"
 #include "init.h"
@@ -28,17 +29,17 @@ typedef struct nstrophy_parameters {
 // usage message
 int print_usage();
 // print parameters
-int print_params(nstrophy_parameters parameters, unsigned int driving, unsigned int init, FILE* file);
+int print_params(nstrophy_parameters parameters, unsigned int driving, unsigned int init, char* initfile_str, FILE* file);
 
 // read command line arguments
-int read_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* init, unsigned int* nthreads);
+int read_args(int argc, const char* argv[], char** params, unsigned int* driving_force, unsigned int* command, unsigned int* init, unsigned int* nthreads, char** savefile_str, char** initfile_str);
 int read_params(char* param_str, nstrophy_parameters* parameters);
 int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2);
 
 // set driving force
 _Complex double* set_driving(unsigned int driving, nstrophy_parameters parameters);
 // set initial condition
-_Complex double* set_init(unsigned int init, nstrophy_parameters parameters);
+_Complex double* set_init(unsigned int init, nstrophy_parameters parameters, FILE* initfile);
 
 #define COMMAND_UK 1
 #define COMMAND_EEA 2
@@ -49,6 +50,7 @@ _Complex double* set_init(unsigned int init, nstrophy_parameters parameters);
 
 #define INIT_RANDOM 1
 #define INIT_GAUSSIAN 2
+#define INIT_FILE 3
 
 
 int main (
@@ -62,13 +64,17 @@ int main (
   unsigned int nthreads=1;
   _Complex double* u0;
   _Complex double *g;
+  char* savefile_str=NULL;
+  char* initfile_str=NULL;
+  FILE* savefile=NULL;
+  FILE* initfile=NULL;
 
   command=0;
   driving=0;
   init=0;
 
   // read command line arguments
-  ret=read_args(argc, argv, &param_str, &driving, &init, &command, &nthreads);
+  ret=read_args(argc, argv, &param_str, &driving, &init, &command, &nthreads, &savefile_str, &initfile_str);
   if(ret<0){
     return(-1);
   }
@@ -79,24 +85,47 @@ int main (
     return(-1);
   }
 
+  // open initfile
+  if(initfile_str!=NULL){
+    initfile=fopen(initfile_str,"r");
+    if(initfile==NULL){
+      fprintf(stderr,"Error opening file '%s' for reading: %s\n", initfile_str, strerror(errno));
+      return(-1);
+    }
+  }
+
   // set driving force
   g=set_driving(driving, parameters);
   // set initial condition
-  u0=set_init(init, parameters);
+  u0=set_init(init, parameters, initfile);
+
+  // close initfile (do this early, so that it is possible to use the same file for init and save)
+  if (initfile!=NULL){
+    fclose(initfile);
+  }
+
+  // open savefile (do this after closing init file)
+  if(savefile_str!=NULL){
+    savefile=fopen(savefile_str,"w");
+    if(savefile==NULL){
+      fprintf(stderr,"Error opening file '%s' for writing: %s\n", savefile_str, strerror(errno));
+      return(-1);
+    }
+  }
 
   // print parameters
-  print_params(parameters, driving, init, stderr);
-  print_params(parameters, driving, init, stdout);
+  print_params(parameters, driving, init, initfile_str, stderr);
+  print_params(parameters, driving, init, initfile_str, stdout);
 
   // 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.print_freq, nthreads);
+    uk(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads, savefile);
   }
   else if(command==COMMAND_EEA){
-    eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads);
+    eea(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, parameters.print_freq, nthreads, savefile);
   }
   else if(command==COMMAND_QUIET){
-    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, nthreads);
+    quiet(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, u0, g, nthreads, savefile);
   }
   else if(command==0){
     fprintf(stderr, "error: no command specified\n");
@@ -106,6 +135,11 @@ int main (
   free(g);
   free(u0);
 
+  // close savefile
+  if (savefile!=NULL){
+    fclose(savefile);
+  }
+
   return(0);
 }
 
@@ -120,6 +154,7 @@ int print_params(
   nstrophy_parameters parameters,
   unsigned int driving,
   unsigned int init,
+  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);
@@ -143,6 +178,9 @@ int print_params(
   case INIT_GAUSSIAN:
     fprintf(file,", init=gaussian");
     break;
+  case INIT_FILE:
+    fprintf(file,", init=file:%s", initfile_str);
+    break;
   default:
     fprintf(file,", init=gaussian");
     break;
@@ -158,6 +196,7 @@ int print_params(
 #define CP_FLAG_DRIVING 2
 #define CP_FLAG_NTHREADS 3
 #define CP_FLAG_INIT 4
+#define CP_FLAG_SAVEFILE 5
 int read_args(
   int argc,
   const char* argv[],
@@ -165,7 +204,9 @@ int read_args(
   unsigned int* driving_force,
   unsigned int* init,
   unsigned int* command,
-  unsigned int* nthreads
+  unsigned int* nthreads,
+  char** savefile_str,
+  char** initfile_str
 ){
   int i;
   int ret;
@@ -192,6 +233,9 @@ int read_args(
 	case 'i':
 	  flag=CP_FLAG_INIT;
 	  break;
+	case 's':
+	  flag=CP_FLAG_SAVEFILE;
+	  break;
 	default:
 	  fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
 	  print_usage();
@@ -236,12 +280,22 @@ int read_args(
       else if (strcmp(argv[i],"gaussian")==0){
 	*init=INIT_GAUSSIAN;
       }
+      // matches any argument that starts with 'file:'
+      else if (strncmp(argv[i],"file:",5)==0){
+	*init=INIT_FILE;
+	*initfile_str=(char*)argv[i]+5;
+      }
       else{
 	fprintf(stderr, "error: unrecognized initial condition '%s'\n",argv[i]);
 	return(-1);
       }
       flag=0;
     }
+    // savefile
+    else if(flag==CP_FLAG_SAVEFILE){
+      *savefile_str=(char*)argv[i];
+      flag=0;
+    }
     // computation to run
     else{
       if(strcmp(argv[i], "uk")==0){
@@ -495,7 +549,8 @@ _Complex double* set_driving(
 // set initial condition
 _Complex double* set_init(
   unsigned int init,
-  nstrophy_parameters parameters
+  nstrophy_parameters parameters,
+  FILE* initfile
 ){
   _Complex double* u0=calloc(sizeof(_Complex double),(2*parameters.K1+1)*(2*parameters.K2+1));
 
@@ -508,6 +563,10 @@ _Complex double* set_init(
     init_gaussian(u0, parameters.K1, parameters.K2);
     break;
 
+  case INIT_FILE:
+    init_file(u0, parameters.K1, parameters.K2, initfile);
+    break;
+
   default:
     init_gaussian(u0, parameters.K1, parameters.K2);
     break;
diff --git a/src/navier-stokes.c b/src/navier-stokes.c
index 7c4ecb3..12cba85 100644
--- a/src/navier-stokes.c
+++ b/src/navier-stokes.c
@@ -1,4 +1,5 @@
 #include "navier-stokes.h"
+#include "io.h"
 #include <math.h>
 #include <stdlib.h>
 
@@ -15,7 +16,8 @@ int uk(
   _Complex double* u0,
   _Complex double* g,
   unsigned int print_freq,
-  unsigned int nthreads
+  unsigned int nthreads,
+  FILE* savefile
 ){
   _Complex double* u;
   _Complex double* tmp1;
@@ -65,6 +67,9 @@ int uk(
     }
   }
 
+  // save final entry to savefile
+  write_u(u, K1, K2, savefile);
+
   ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
   return(0);
 }
@@ -82,7 +87,8 @@ int eea(
   _Complex double* u0,
   _Complex double* g,
   unsigned int print_freq,
-  unsigned int nthreads
+  unsigned int nthreads,
+  FILE* savefile
 ){
   _Complex double* u;
   _Complex double* tmp1;
@@ -126,6 +132,9 @@ int eea(
     }
   }
 
+  // save final entry to savefile
+  write_u(u, K1, K2, savefile);
+
   ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
   return(0);
 }
@@ -142,7 +151,8 @@ int quiet(
   double L,
   _Complex double* u0,
   _Complex double* g,
-  unsigned int nthreads
+  unsigned int nthreads,
+  FILE* savefile
 ){
   _Complex double* u;
   _Complex double* tmp1;
@@ -162,6 +172,9 @@ int quiet(
     ins_step(u, K1, K2, N1, N2, nu, delta, L, g, fft1, fft2, ifft, tmp1, tmp2, tmp3);
   }
 
+  // save final entry to savefile
+  write_u(u, K1, K2, savefile);
+
   ns_free_tmps(u, tmp1, tmp2, tmp3, fft1, fft2, ifft);
   return(0);
 }
diff --git a/src/navier-stokes.h b/src/navier-stokes.h
index 7841f8f..92ff754 100644
--- a/src/navier-stokes.h
+++ b/src/navier-stokes.h
@@ -13,16 +13,13 @@ typedef struct fft_vects {
 } fft_vect;
 
 // compute u_k
-int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreads);
-
-// compute the energy as a function of time
-int energy( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreads);
+int uk( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreadsl, FILE* savefile);
 
 // compute energy, enstrophy and alpha
-int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreads);
+int eea( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int print_freq, unsigned int nthreads, FILE* savefile);
 
 // 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, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int nthreads);
+int quiet( int K1, int K2, int N1, int N2, unsigned int nsteps, double nu, double delta, double L, _Complex double* u0, _Complex double* g, unsigned int nthreads, FILE* savefile);
 
 
 // initialize vectors for computation