From 0cdb914b5764f692189ed2bc395e3b09ead758e4 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Fri, 27 May 2022 16:09:17 -0400 Subject: savefile and initfile --- src/init.c | 12 ++++++ src/init.h | 5 +++ src/io.c | 112 ++++++++++++++++++++++++++++++++++++++++++++++++++++ src/io.h | 12 ++++++ src/main.c | 83 ++++++++++++++++++++++++++++++++------ src/navier-stokes.c | 19 +++++++-- src/navier-stokes.h | 9 ++--- 7 files changed, 231 insertions(+), 21 deletions(-) create mode 100644 src/io.c create mode 100644 src/io.h (limited to 'src') 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 #include @@ -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 + // 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 +#include +#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 + +// 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 #include #include +#include #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, ¶m_str, &driving, &init, &command, &nthreads); + ret=read_args(argc, argv, ¶m_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 #include @@ -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 -- cgit v1.2.3-54-g00ecf