diff --git a/README.md b/README.md index 1b87cd8..c0602f8 100644 --- a/README.md +++ b/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 - `equation=irreversible` or of the enstrophy if `equation=reversible`. +* `init_energy` (double, default is to not rescale the initial condition, is + 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. diff --git a/src/constants.cpp b/src/constants.cpp index ca27947..649d916 100644 --- a/src/constants.cpp +++ b/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 diff --git a/src/init.c b/src/init.c index 05283bb..97cceeb 100644 --- a/src/init.c +++ b/src/init.c @@ -17,21 +17,17 @@ limitations under the License. #include "init.h" #include "navier-stokes.h" #include "io.h" -#include #include +#include // random initial condition int init_random ( _Complex double* u0, - double init_en, int K1, int K2, - double L, - int seed, - bool irreversible + int seed ){ 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, - double L, - bool irreversible + int K2 ){ 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; } diff --git a/src/init.h b/src/init.h index dea15a9..e29395d 100644 --- a/src/init.h +++ b/src/init.h @@ -21,10 +21,10 @@ limitations under the License. #include // 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); diff --git a/src/main.c b/src/main.c index d1b045f..6ea2216 100644 --- a/src/main.c +++ b/src/main.c @@ -16,6 +16,7 @@ limitations under the License. #define VERSION "0.1" +#include #include #include #include @@ -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(¶meters); + // 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(¶m_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(¶m_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){ - ret=sscanf(rhs,"%lf",&(parameters->init_en)); + else if (strcmp(lhs,"init_enstrophy")==0){ + // 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; }