diff --git a/README.md b/README.md index fa2f6f0..ba664a9 100644 --- a/README.md +++ b/README.md @@ -57,8 +57,8 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are * `N2` (int, default `N`): same as `N` but only for y component. -* `nsteps` (long int, default 10000000): number of steps in the computation. - Set to 0 to keep on going forever. +* `final_time` (double, default 100000): time at which to end the computation. + Set to <0 to keep on going forever. * `nu` (double, default 0.00048828125): viscosity. @@ -66,11 +66,10 @@ should be a `;` sperated list of `key=value` pairs. The possible keys are * `L` (double, default 2pi): size of box. -* `print_freq` (long int, default 1000): only print every `print_freq` steps. +* `print_freq` (double, default 1): only print when time crosses integer + multiples of `print_freq`. -* `starting_step` (long int, default 0): start the computation at this step. - -* `starting_time` (double, default starting_time*delta): start the computation at this time. +* `starting_time` (double, default 0): starting time. * `driving`: either `zero` for no driving, `test` (default) for a testing driving force or `file:` or `file_txt:` to read the diff --git a/src/main.c b/src/main.c index bc4cc0c..8cea261 100644 --- a/src/main.c +++ b/src/main.c @@ -23,7 +23,6 @@ limitations under the License. #include #include #include -#include #include #include "constants.cpp" @@ -42,15 +41,14 @@ typedef struct nstrophy_parameters { int K2; int N1; int N2; - unsigned int nsteps; + double final_time; double nu; double delta; double L; double adaptive_tolerance; double adaptive_factor; - uint64_t print_freq; + double print_freq; int seed; - uint64_t starting_step; double starting_time; unsigned int driving; unsigned int init; @@ -67,7 +65,7 @@ int print_params(nstrophy_parameters parameters, char* initfile_str, char* drivi // read command line arguments int read_args(int argc, const char* argv[], char** params, unsigned int* command, unsigned int* nthreads, char** savefile_str); int read_params(char* param_str, nstrophy_parameters* parameters, char** initfile_str, char** drivingfile_str); -int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, bool* set_starting_time, char** initfile_str, char** drivingfile_str); +int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, char** initfile_str, char** drivingfile_str); // set driving force _Complex double* set_driving(nstrophy_parameters parameters); @@ -172,15 +170,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, parameters.adaptive_tolerance, parameters.adaptive_factor, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_step, 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, u0, g, parameters.irreversible, 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); - enstrophy(parameters.K1, parameters.K2, parameters.N1, parameters.N2, parameters.nsteps, parameters.nu, parameters.delta, parameters.L, parameters.adaptive_tolerance, parameters.adaptive_factor, u0, g, parameters.irreversible, parameters.algorithm, parameters.print_freq, parameters.starting_step, parameters.starting_time, nthreads, savefile, (char*)argv[0], param_str, savefile_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, 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.adaptive_tolerance, parameters.adaptive_factor, parameters.starting_step, u0, g, parameters.irreversible, 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.starting_time, u0, g, parameters.irreversible, parameters.algorithm, nthreads, savefile); } else if(command==0){ fprintf(stderr, "error: no command specified\n"); @@ -377,7 +375,6 @@ int read_params( // whether N was set explicitly bool setN1=false; bool setN2=false; - bool set_starting_time=false; // whether lhs (false is rhs) bool lhs=true; @@ -393,10 +390,9 @@ int read_params( parameters->L=2*M_PI; parameters->adaptive_tolerance=1e-11; parameters->adaptive_factor=0.9; - parameters->nsteps=10000000; - parameters->print_freq=1000; + parameters->final_time=100000; + parameters->print_freq=1; parameters->starting_time=0; - parameters->starting_step=0; parameters->seed=17; parameters->driving=DRIVING_TEST; parameters->drivingfile=NULL; @@ -423,7 +419,7 @@ int read_params( break; case ';': //set parameter - ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &set_starting_time, initfile_str, drivingfile_str); + ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str, drivingfile_str); if(ret<0){ return ret; } @@ -450,7 +446,7 @@ int read_params( // set last param if (*param_str!='\0'){ - ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, &set_starting_time, initfile_str, drivingfile_str); + ret=set_parameter(buffer_lhs, buffer_rhs, parameters, &setN1, &setN2, initfile_str, drivingfile_str); if(ret<0){ return ret; } @@ -469,10 +465,6 @@ int read_params( parameters->N2=smallest_pow2(3*(parameters->K2)); } - // if starting_time not set explicitly and algorithm is not adaptive, set it to delta*starting_step - if (!set_starting_time && parameters->algorithmstarting_time=parameters->delta*parameters->starting_step; - } return(0); } @@ -484,7 +476,6 @@ int set_parameter( nstrophy_parameters* parameters, bool* setN1, bool* setN2, - bool* set_starting_time, char** initfile_str, char** drivingfile_str ){ @@ -557,10 +548,10 @@ int set_parameter( *setN1=true; *setN2=true; } - else if (strcmp(lhs,"nsteps")==0){ - ret=sscanf(rhs,"%u",&(parameters->nsteps)); + else if (strcmp(lhs,"final_time")==0){ + ret=sscanf(rhs,"%lf",&(parameters->final_time)); if(ret!=1){ - fprintf(stderr, "error: parameter 'nsteps' should be an unsigned integer\n got '%s'\n",rhs); + fprintf(stderr, "error: parameter 'final_time' should be a double\n got '%s'\n",rhs); return(-1); } } @@ -600,9 +591,9 @@ int set_parameter( } } else if (strcmp(lhs,"print_freq")==0){ - ret=sscanf(rhs,"%lu",&(parameters->print_freq)); + ret=sscanf(rhs,"%lf",&(parameters->print_freq)); if(ret!=1){ - fprintf(stderr, "error: parameter 'print_freq' should be an unsigned integer\n got '%s'\n",rhs); + fprintf(stderr, "error: parameter 'print_freq' should be a double\n got '%s'\n",rhs); return(-1); } } @@ -613,20 +604,12 @@ int set_parameter( return(-1); } } - else if (strcmp(lhs,"starting_step")==0){ - ret=sscanf(rhs,"%lu",&(parameters->starting_step)); - if(ret!=1){ - fprintf(stderr, "error: parameter 'starting_step' should be a long unsigned integer\n got '%s'\n",rhs); - return(-1); - } - } else if (strcmp(lhs,"starting_time")==0){ ret=sscanf(rhs,"%lf",&(parameters->starting_time)); if(ret!=1){ fprintf(stderr, "error: parameter 'starting_time' should be a double\n got '%s'\n",rhs); return(-1); } - *set_starting_time=true; } else if (strcmp(lhs,"driving")==0){ if (strcmp(rhs,"zero")==0){ @@ -747,11 +730,11 @@ _Complex double* set_init( case INIT_FILE: init_file_bin(u0, parameters->K1, parameters->K2, parameters->initfile); + // read start time + fread(&(parameters->starting_time), sizeof(double), 1, parameters->initfile); if(parameters->algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){ // read delta fread(&(parameters->delta), sizeof(double), 1, parameters->initfile); - // read start time - fread(&(parameters->starting_time), sizeof(double), 1, parameters->initfile); } break; diff --git a/src/navier-stokes.c b/src/navier-stokes.c index 8594e5c..12a2944 100644 --- a/src/navier-stokes.c +++ b/src/navier-stokes.c @@ -17,8 +17,8 @@ limitations under the License. #include "constants.cpp" #include "io.h" #include "navier-stokes.h" -#include "statistics.h" #include +#include #include #include @@ -28,7 +28,7 @@ int uk( int K2, int N1, int N2, - uint64_t nsteps, + double final_time, double nu, double delta, double L, @@ -38,13 +38,11 @@ int uk( _Complex double* g, bool irreversible, unsigned int algorithm, - uint64_t print_freq, - uint64_t starting_step, + double print_freq, double starting_time, unsigned int nthreads, FILE* savefile ){ - double time=starting_time; _Complex double* u; _Complex double* tmp1; _Complex double* tmp2; @@ -53,7 +51,7 @@ int uk( _Complex double* tmp5; _Complex double* tmp6; _Complex double* tmp7; - uint64_t t; + double time; fft_vect fft1; fft_vect fft2; fft_vect ifft; @@ -65,18 +63,23 @@ int uk( // print column headers printf("# 1:i 2:t "); - t=3; + unsigned int i=3; for(kx=-K1;kx<=K1;kx++){ for (ky=-K2;ky<=K2;ky++){ - printf(" %6lu:(%4d,%4d)r ",t,kx,ky); - t++; - printf(" %6lu:(%4d,%4d)i ",t,kx,ky); - t++; + printf(" %6u:(%4d,%4d)r ",i,kx,ky); + i++; + printf(" %6u:(%4d,%4d)i ",i,kx,ky); + i++; } } + // period + // add 0.1 to ensure proper rounding + uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1); + // iterate - for(t=starting_step;nsteps==0 || t(n+1)*print_freq){ + n++; + + fprintf(stderr,"% .8e ",time); + printf("% .15e ",time); for(kx=-K1;kx<=K1;kx++){ for (ky=-K2;ky<=K2;ky++){ @@ -120,7 +125,7 @@ int enstrophy( int K2, int N1, int N2, - uint64_t nsteps, + double final_time, double nu, double delta, double L, @@ -130,8 +135,7 @@ int enstrophy( _Complex double* g, bool irreversible, unsigned int algorithm, - uint64_t print_freq, - uint64_t starting_step, + double print_freq, double starting_time, unsigned int nthreads, FILE* savefile, @@ -148,11 +152,11 @@ int enstrophy( _Complex double* tmp5; _Complex double* tmp6; _Complex double* tmp7; - double time=starting_time; + double time; double alpha, enstrophy; + double prevtime; double avg_a,avg_en,avg_en_x_a; // index - uint64_t t; fft_vect fft1; fft_vect fft2; fft_vect ifft; @@ -166,12 +170,15 @@ int enstrophy( avg_a=0; avg_en=0; avg_en_x_a=0; + prevtime=starting_time; - // special first case when starting_time is not a multiple of print_freq - uint64_t first_box = print_freq - (starting_step % print_freq); + // period + // add 0.1 to ensure proper rounding + uint64_t n=(uint64_t)((starting_time-fmod(starting_time, print_freq))/print_freq+0.1); // iterate - for(t=starting_step;nsteps==0 || t(n+1)*print_freq){ + n++; + + // adjust duration of interval to its actual value + avg_a*=print_freq/(time-prevtime); + avg_en*=print_freq/(time-prevtime); + avg_en_x_a*=print_freq/(time-prevtime); - if(t>starting_step && t%print_freq==0){ // print to stderr so user can follow along if(algorithm==ALGORITHM_RKF45){ - fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta); + fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy, delta); } else { - fprintf(stderr,"%lu % .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); + fprintf(stderr,"% .8e % .8e % .8e % .8e % .8e % .8e % .8e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); } // print to stdout - printf("%8lu % .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",t,time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); + printf("% .15e % .15e % .15e % .15e % .15e % .15e % .15e\n",time, avg_a, avg_en_x_a, avg_en, alpha, alpha*enstrophy, enstrophy); } + // reset averages + avg_a=0; + avg_en=0; + avg_en_x_a=0; + prevtime=time; + // catch abort signal if (g_abort){ // print u to stderr if no savefile @@ -223,22 +244,20 @@ int enstrophy( if(params_string!=NULL) { char* params=calloc(sizeof(char), strlen(params_string)+1); strcpy(params, params_string); - remove_entry(params, "starting_step"); remove_entry(params, "starting_time"); remove_entry(params, "init"); - remove_entry(params, "nsteps"); if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){ remove_entry(params, "delta"); } - fprintf(savefile," -p \"%s;starting_step=%lu;nsteps=%lu;init=file:%s", params, t+1, (nsteps+starting_step < t+1 ? 0 : nsteps+starting_step-t-1), savefile_string); + fprintf(savefile," -p \"%s;init=file:%s", params, savefile_string); free(params); // write delta if adaptive, and not writing binary if(algorithm>ALGORITHM_ADAPTIVE_THRESHOLD && (savefile==stderr || savefile==stdout)){ - fprintf(savefile,";delta=%.15e;starting_time=%.15e", delta, starting_time); + fprintf(savefile,";delta=%.15e", delta); } - // write starting_time if not adaptive - if(algorithmALGORITHM_ADAPTIVE_THRESHOLD){ - // first binary entry: delta fwrite(&delta, sizeof(double), 1, savefile); - // ssecond binary entry: starting time - fwrite(&starting_time, sizeof(double), 1, savefile); } } } @@ -269,13 +287,13 @@ int quiet( int K2, int N1, int N2, - uint64_t nsteps, + double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, - uint64_t starting_step, + double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, @@ -291,7 +309,7 @@ int quiet( _Complex double* tmp5; _Complex double* tmp6; _Complex double* tmp7; - uint64_t t; + double time; fft_vect fft1; fft_vect fft2; fft_vect ifft; @@ -301,7 +319,8 @@ int quiet( copy_u(u, u0, K1, K2); // iterate - for(t=starting_step;nsteps==0 || t #include -#include #include #define M_PI 3.14159265358979323846 @@ -34,13 +33,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, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_step, double starting_time, unsigned int nthreadsl, FILE* savefile); +int uk( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double starting_time, unsigned int nthreadsl, FILE* savefile); // compute enstrophy and alpha -int enstrophy( int K1, int K2, int N1, int N2, uint64_t nsteps, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, uint64_t print_freq, uint64_t starting_step, double starting_time, unsigned int nthreads, FILE* savefile, char* cmd_string, char* params_string, char* savefile_string); +int enstrophy( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, double print_freq, double 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, double adaptive_tolerance, double adaptive_factor, uint64_t starting_step, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile); +int quiet( int K1, int K2, int N1, int N2, double final_time, double nu, double delta, double L, double adaptive_tolerance, double adaptive_factor, double starting_time, _Complex double* u0, _Complex double* g, bool irreversible, unsigned int algorithm, unsigned int nthreads, FILE* savefile); // initialize vectors for computation diff --git a/src/statistics.c b/src/statistics.c deleted file mode 100644 index 2acfca4..0000000 --- a/src/statistics.c +++ /dev/null @@ -1,44 +0,0 @@ -/* -Copyright 2017-2023 Ian Jauslin - -Licensed under the Apache License, Version 2.0 (the "License"); -you may not use this file except in compliance with the License. -You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - -Unless required by applicable law or agreed to in writing, software -distributed under the License is distributed on an "AS IS" BASIS, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the License for the specific language governing permissions and -limitations under the License. -*/ - -#include "statistics.h" - -// run this at each step to compute the running average -double average_step( - // the value at the step - double val, - // the average computed so far - double avg, - uint64_t t, - uint64_t starting_time, - uint64_t print_freq, - uint64_t first_box -){ - - // running average - // reset averages - if(t % print_freq == 1){ - return(0); - } - - // compute average - // different computationin first box if starting_time is not a multiple of print_freq - if(t < starting_time + first_box){ - return(avg+val/first_box); - } else { - return(avg+val/print_freq); - } -} diff --git a/src/statistics.h b/src/statistics.h deleted file mode 100644 index ace97c5..0000000 --- a/src/statistics.h +++ /dev/null @@ -1,25 +0,0 @@ -/* -Copyright 2017-2023 Ian Jauslin - -Licensed under the Apache License, Version 2.0 (the "License"); -you may not use this file except in compliance with the License. -You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - -Unless required by applicable law or agreed to in writing, software -distributed under the License is distributed on an "AS IS" BASIS, -WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -See the License for the specific language governing permissions and -limitations under the License. -*/ - -#ifndef STATISTICS_H -#define STATISTICS_H - -#include - -// run this at each step to compute the running average -double average_step( double val, double avg, uint64_t t, uint64_t starting_time, uint64_t print_freq, uint64_t first_box); - -#endif