Nstrophy/src/main.c

/*
Copyright 2017-2024 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.
*/

#define VERSION "0.1"

#include <math.h>
#include <signal.h>
#include <string.h>
#include <stdlib.h>
#include <errno.h>
#include <wordexp.h>

#include "constants.cpp"
#include "driving.h"
#include "dstring.h"
#include "init.h"
#include "int_tools.h"
#include "lyapunov.h"
#include "navier-stokes.h"


// structure to store parameters, to make it easier to pass parameters to CLI functions
typedef struct nstrophy_parameters {
  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;
  int N1;
  int N2;
  double final_time;
  double nu;
  double delta;
  double L;
  double adaptive_tolerance;
  double adaptive_factor;
  double max_delta;
  unsigned int adaptive_cost;
  double print_freq;
  int seed;
  double starting_time;
  unsigned int driving;
  unsigned int init;
  unsigned int algorithm;
  unsigned int algorithm_lyapunov;
  bool keep_en_cst;
  FILE* initfile;
  FILE* drivingfile;
  double lyapunov_reset;
  unsigned int lyapunov_trigger;
  double D_epsilon;
  bool print_alpha;
} nstrophy_parameters;

// usage message
int print_usage();
// print parameters
int print_params(nstrophy_parameters parameters, char* initfile_str, char* drivingfile_str, FILE* file);

// read command line arguments
int read_args(int argc, const char* argv[], dstring* params, unsigned int* command, unsigned int* nthreads, dstring* savefile_str, dstring* utfile_str, dstring* resumefile_str);
int set_default_params(nstrophy_parameters* parameters);
int read_params(dstring param_str, nstrophy_parameters* parameters, dstring* initfile_str, dstring* drivingfile_str);
int set_parameter(char* lhs, char* rhs, nstrophy_parameters* parameters, bool* setN1, bool* setN2, dstring* initfile_str, dstring* drivingfile_str);
// read args from file
int args_from_file(dstring* params, unsigned int* command, unsigned int* nthreads, dstring* savefile_str, dstring* utfile_str, char* file_str);

// set driving force
_Complex double* set_driving(nstrophy_parameters parameters);
// set initial condition
_Complex double* set_init(nstrophy_parameters* parameters);

// signal handler
void sig_handler (int signo);

// global variable to handle interrupts
volatile bool g_abort = false;
// signal handler
void sig_handler (int signo){
  if (signo == SIGINT){
    fprintf(stderr,"received signal SIGINT, interrupting computation\n");
    g_abort = true;
  }
  else if (signo == SIGTERM){
    fprintf(stderr,"received signal SIGTERM, interrupting computation\n");
    g_abort = true;
  }
  else{
    fprintf(stderr,"received signal %d\n",signo);
  }
}


int main (
  int argc,
  const char* argv[]
){
  dstring param_str;
  nstrophy_parameters parameters;
  int ret;
  unsigned int command;
  unsigned int nthreads=1;
  _Complex double* u0;
  _Complex double *g;
  dstring savefile_str;
  dstring utfile_str;
  dstring initfile_str;
  dstring drivingfile_str;
  dstring resumefile_str;
  FILE* savefile=NULL;
  FILE* utfile=NULL;

  command=0;

  // init strings
  dstring_init(&param_str, 64);
  dstring_init(&savefile_str, 64);
  dstring_init(&utfile_str, 64);
  dstring_init(&initfile_str, 64);
  dstring_init(&drivingfile_str, 64);
  dstring_init(&resumefile_str, 64);

  // read command line arguments
  ret=read_args(argc, argv, &param_str, &command, &nthreads, &savefile_str, &utfile_str, &resumefile_str);
  if(ret<0){
    dstring_free(param_str);
    dstring_free(savefile_str);
    dstring_free(utfile_str);
    dstring_free(initfile_str);
    dstring_free(drivingfile_str);
    dstring_free(resumefile_str);
    return(ret);
  }

  // set default params
  set_default_params(&parameters);
  // read params
  ret=read_params(param_str, &parameters, &initfile_str, &drivingfile_str);
  if(ret<0){
    dstring_free(param_str);
    dstring_free(savefile_str);
    dstring_free(utfile_str);
    dstring_free(initfile_str);
    dstring_free(drivingfile_str);
    dstring_free(resumefile_str);
    return(ret);
  }

  // if command is 'resume', then read args from file
  if(command==COMMAND_RESUME){
    ret=args_from_file(&param_str, &command, &nthreads, &savefile_str, &utfile_str, dstring_to_str_noinit(&resumefile_str));
    if(ret<0){
      dstring_free(param_str);
      dstring_free(savefile_str);
      dstring_free(utfile_str);
      dstring_free(initfile_str);
      dstring_free(drivingfile_str);
      dstring_free(resumefile_str);
      return(ret);
    }
    // read params
    ret=read_params(param_str, &parameters, &initfile_str, &drivingfile_str);
    if(ret<0){
      dstring_free(param_str);
      dstring_free(savefile_str);
      dstring_free(utfile_str);
      dstring_free(initfile_str);
      dstring_free(drivingfile_str);
      dstring_free(resumefile_str);
      return(ret);
    }

    // reread arguments (to allow overrides from the command line, but do not override the command)
    unsigned int dummy_command;
    read_args(argc, argv, &param_str, &dummy_command, &nthreads, &savefile_str, &utfile_str, &resumefile_str);
    // reread params
    read_params(param_str, &parameters, &initfile_str, &drivingfile_str);
  }

  // free strings
  dstring_free(resumefile_str);

  // open initfile
  if(initfile_str.length!=0){
    parameters.initfile=fopen(dstring_to_str_noinit(&initfile_str),"r");
    if(parameters.initfile==NULL){
      fprintf(stderr,"Error opening file '%s' for reading: %s\n", dstring_to_str_noinit(&initfile_str), strerror(errno));
      dstring_free(param_str);
      dstring_free(savefile_str);
      dstring_free(utfile_str);
      dstring_free(initfile_str);
      dstring_free(drivingfile_str);
      return(-1);
    }
  }
  // open drivingfile
  if(drivingfile_str.length!=0){
    parameters.drivingfile=fopen(dstring_to_str_noinit(&drivingfile_str),"r");
    if(parameters.drivingfile==NULL){
      fprintf(stderr,"Error opening file '%s' for reading: %s\n", dstring_to_str_noinit(&drivingfile_str), strerror(errno));
      dstring_free(param_str);
      dstring_free(savefile_str);
      dstring_free(utfile_str);
      dstring_free(initfile_str);
      dstring_free(drivingfile_str);
      return(-1);
    }
  }

  // set driving force
  g=set_driving(parameters);
  // set initial condition
  u0=set_init(&parameters);

  // 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);
  }
  // close drivingfile
  if (parameters.drivingfile!=NULL){
    fclose(parameters.drivingfile);
  }

  // open savefile (do this after closing init file)
  if(savefile_str.length!=0){
    savefile=fopen(dstring_to_str_noinit(&savefile_str),"w");
    if(savefile==NULL){
      fprintf(stderr,"Error opening file '%s' for writing: %s\n", dstring_to_str_noinit(&savefile_str), strerror(errno));
      dstring_free(param_str);
      dstring_free(savefile_str);
      dstring_free(utfile_str);
      dstring_free(initfile_str);
      dstring_free(drivingfile_str);
      free(g);
      free(u0);
      return(-1);
    }
  }

  // open utfile (do this after closing init file)
  if(utfile_str.length!=0){
    utfile=fopen(dstring_to_str_noinit(&utfile_str),"w");
    if(utfile==NULL){
      fprintf(stderr,"Error opening file '%s' for writing: %s\n", dstring_to_str_noinit(&utfile_str), strerror(errno));
      dstring_free(param_str);
      dstring_free(savefile_str);
      dstring_free(utfile_str);
      dstring_free(initfile_str);
      dstring_free(drivingfile_str);
      free(g);
      free(u0);
      return(-1);
    }
  }

  // print parameters
  print_params(parameters, dstring_to_str_noinit(&initfile_str), dstring_to_str_noinit(&drivingfile_str), stderr);
  print_params(parameters, dstring_to_str_noinit(&initfile_str), dstring_to_str_noinit(&drivingfile_str), stdout);

  // free strings
  dstring_free(initfile_str);
  dstring_free(drivingfile_str);

  // 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_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_enstrophy, parameters.algorithm, parameters.print_freq, parameters.starting_time, parameters.print_alpha, nthreads, savefile, utfile, (char*)argv[0], dstring_to_str_noinit(&param_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_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.lyapunov_trigger, 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.algorithm_lyapunov, parameters.starting_time, nthreads);
  }
  else if(command==0){
    fprintf(stderr, "error: no command specified\n");
    print_usage();
  }

  free(g);
  free(u0);

  // free strings
  dstring_free(param_str);
  dstring_free(savefile_str);
  dstring_free(utfile_str);

  // close savefile
  if (savefile!=NULL){
    fclose(savefile);
  }

  // close utfile
  if (utfile!=NULL){
    fclose(utfile);
  }

  return(0);
}

// usage message
int print_usage(){
  fprintf(stderr, "usage:\n       nstrophy [-t nthreads] [-p parameters] [-s savefile] [-u u_outfile] <command>\n  where <command> is one of\n    enstrophy\n    uk\n    quiet\n    resume <resume_file>\n\n");
  return(0);
}

// print parameters
int print_params(
  nstrophy_parameters parameters,
  char* initfile_str,
  char* drivingfile_str,
  FILE* file
){
  fprintf(file,"# ");

  if (parameters.irreversible){
    fprintf(file,"equation=irreversible");
  } else {
    fprintf(file,"equation=reversible");
  }

  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:
    fprintf(file,", driving=test");
    break;
  case DRIVING_ZERO:
    fprintf(file,", driving=zero");
    break;
  case DRIVING_FILE:
    fprintf(file,", driving=file:%s", drivingfile_str);
    break;
  case DRIVING_FILE_TXT:
    fprintf(file,", driving=file_txt:%s", drivingfile_str);
    break;
  default:
    fprintf(file,", driving=test");
    break;
  }

  switch(parameters.init){
  case INIT_RANDOM:
    fprintf(file,", init=random");
    break;
  case INIT_GAUSSIAN:
    fprintf(file,", init=gaussian");
    break;
  case INIT_FILE:
    fprintf(file,", init=file:%s", initfile_str);
    break;
  case INIT_FILE_TXT:
    fprintf(file,", init=file_txt:%s", initfile_str);
    break;
  default:
    fprintf(file,", init=gaussian");
    break;
  }

  switch(parameters.algorithm){
  case ALGORITHM_RK4:
    fprintf(file,", algorithm=RK4");
    break;
  case ALGORITHM_RK2:
    fprintf(file,", algorithm=RK2");
    break;
  case ALGORITHM_RKF45:
    fprintf(file,", algorithm=RKF45, tolerance=%.15e",parameters.adaptive_tolerance);
    break;
  case ALGORITHM_RKDP54:
    fprintf(file,", algorithm=RKDP54, tolerance=%.15e",parameters.adaptive_tolerance);
    break;
  case ALGORITHM_RKBS32:
    fprintf(file,", algorithm=RKBS32, tolerance=%.15e",parameters.adaptive_tolerance);
    break;
  default:
    fprintf(file,", algorithm=RK4");
    break;
  }

  if(parameters.algorithm>ALGORITHM_ADAPTIVE_THRESHOLD){
    switch(parameters.adaptive_cost){
    case COST_L1:
      fprintf(file,", cost=L1");
      break;
    case COST_k3:
      fprintf(file,", cost=k3");
      break;
    case COST_k32:
      fprintf(file,", cost=k32");
      break;
    case COST_ENSTROPHY:
      fprintf(file,", cost=enstrophy");
      break;
    case COST_ALPHA:
      fprintf(file,", cost=alpha");
      break;
    }
  }

  fprintf(file,"\n");

  return 0;
}


// read command line arguments
#define CP_FLAG_PARAMS 1
#define CP_FLAG_NTHREADS 2
#define CP_FLAG_SAVEFILE 3
#define CP_FLAG_UTFILE 4
#define CP_FLAG_RESUME 5
int read_args(
  int argc,
  const char* argv[],
  dstring* params,
  unsigned int* command,
  unsigned int* nthreads,
  dstring* savefile_str,
  dstring* utfile_str,
  dstring* resumefile_str
){
  int i;
  int ret;
  // pointers
  char* ptr;
  // flag that indicates what argument is being read
  int flag=0;

  // loop over arguments
  for(i=1;i<argc;i++){
    // flag
    if(argv[i][0]=='-'){
      for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
	switch(*ptr){
	case 'p':
	  flag=CP_FLAG_PARAMS;
	  break;
	case 't':
	  flag=CP_FLAG_NTHREADS;
	  break;
	case 's':
	  flag=CP_FLAG_SAVEFILE;
	  break;
	case 'u':
	  flag=CP_FLAG_UTFILE;
	  break;
	default:
	  fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
	  print_usage();
	  return(-1);
	  break;
	}
      }
    }
    // params
    else if(flag==CP_FLAG_PARAMS){
      str_to_dstring_noinit((char*)argv[i], params);
      flag=0;
    }
    // nthreads
    else if(flag==CP_FLAG_NTHREADS){
      ret=sscanf(argv[i],"%u",nthreads);
      if(ret!=1){
	fprintf(stderr, "error: '-t' should be followed by an unsigned integer\n       got '%s'\n",argv[i]);
	return(-1);
      }
      flag=0;
    }
    // savefile
    else if(flag==CP_FLAG_SAVEFILE){
      str_to_dstring_noinit((char*)argv[i], savefile_str);
      flag=0;
    }
    // savefile
    else if(flag==CP_FLAG_UTFILE){
      str_to_dstring_noinit((char*)argv[i], utfile_str);
      flag=0;
    }
    // resume file
    else if(flag==CP_FLAG_RESUME){
      str_to_dstring_noinit((char*)argv[i], resumefile_str);
      flag=0;
    }
    // computation to run
    else{
      if(strcmp(argv[i], "uk")==0){
	*command=COMMAND_UK;
      }
      else if(strcmp(argv[i], "enstrophy")==0){
	*command=COMMAND_ENSTROPHY;
      }
      else if(strcmp(argv[i], "quiet")==0){
	*command=COMMAND_QUIET;
      }
      else if(strcmp(argv[i], "resume")==0){
	*command=COMMAND_RESUME;
	flag=CP_FLAG_RESUME;
      }
      else if(strcmp(argv[i], "lyapunov")==0){
	*command=COMMAND_LYAPUNOV;
      }
      else{
	fprintf(stderr, "error: unrecognized command: '%s'\n",argv[i]);
	return(-1);
      }
    }
  }

  // check that if the command is 'resume', then resumefile has been set
  if(*command==COMMAND_RESUME && resumefile_str->length==0){
    fprintf(stderr, "error: 'resume' command used, but no resume file\n");
    print_usage();
    return(-1);
  }

  return(0);
}

// set default parameters
int set_default_params(
  nstrophy_parameters* parameters
){
  // 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;
  //delta=2^-13
  parameters->delta=0.0001220703125;
  //nu=2^-11
  parameters->nu=0.00048828125;
  parameters->L=2*M_PI;
  parameters->adaptive_tolerance=1e-11;
  parameters->adaptive_factor=0.9;
  parameters->max_delta=1e-2;
  parameters->adaptive_cost=COST_L1;
  parameters->final_time=100000;
  parameters->print_freq=1;
  parameters->starting_time=0;
  parameters->seed=17;
  parameters->driving=DRIVING_TEST;
  parameters->drivingfile=NULL;
  parameters->init=INIT_GAUSSIAN;
  parameters->initfile=NULL;
  parameters->algorithm=ALGORITHM_RK4;
  parameters->keep_en_cst=false;
  parameters->algorithm_lyapunov=ALGORITHM_RK4;
  // default lyapunov_reset will be print_time (set later) for now set target to 0 to indicate it hasn't been set explicitly
  parameters->lyapunov_trigger=0;

  parameters->print_alpha=false;
  
  return(0);
}

// read parameters string
int read_params(
  dstring param_str,
  nstrophy_parameters* parameters,
  dstring* initfile_str,
  dstring* drivingfile_str
){
  int ret;
  unsigned int i;
  // buffer and associated pointer
  dstring buffer_lhs;
  dstring buffer_rhs;
  // whether N was set explicitly
  bool setN1=false;
  bool setN2=false;
  // which buffer to add to
  dstring* buffer=&buffer_lhs;

  // init
  dstring_init(&buffer_lhs, param_str.length);
  dstring_init(&buffer_rhs, param_str.length);

  for(i=0;i<param_str.length;i++){
    switch(param_str.string[i]){
    case '=':
      // reset buffer
      buffer_rhs.length=0;
      buffer=&buffer_rhs;
      break;
    case ';':
      //set parameter
      ret=set_parameter(dstring_to_str_noinit(&buffer_lhs), dstring_to_str_noinit(&buffer_rhs), parameters, &setN1, &setN2, initfile_str, drivingfile_str);
      if(ret<0){
	return ret;
      }
      // reset buffer
      buffer_lhs.length=0;
      buffer=&buffer_lhs;
      break;
    default:
      // add to buffer
      dstring_append(param_str.string[i],buffer);
      break;
    }
  }

  // set last param
  if (param_str.length!=0){
    ret=set_parameter(dstring_to_str_noinit(&buffer_lhs), dstring_to_str_noinit(&buffer_rhs), parameters, &setN1, &setN2, initfile_str, drivingfile_str);
    if(ret<0){
      return ret;
    }
  }

  // free vects
  dstring_free(buffer_lhs);
  dstring_free(buffer_rhs);

  // if N not set explicitly, set it to the smallest power of 2 that is >3*K+1 (the fft is faster on vectors whose length is a power of 2)
  if (!setN1){
    parameters->N1=smallest_pow2(3*(parameters->K1));
  }
  if (!setN2){
    parameters->N2=smallest_pow2(3*(parameters->K2));
  }

  // if lyapunov_reset or lyapunov_maxu are not set explicitly
  if(parameters->lyapunov_trigger==0){
    parameters->lyapunov_trigger=LYAPUNOV_TRIGGER_TIME;
    parameters->lyapunov_reset=parameters->print_freq;
  }

  return(0);
}


// set a parameter from the parameter string
int set_parameter(
  char* lhs,
  char* rhs,
  nstrophy_parameters* parameters,
  bool* setN1,
  bool* setN2,
  dstring* initfile_str,
  dstring* drivingfile_str
){
  int ret;

  if (strcmp(lhs,"equation")==0){
    if (strcmp(rhs,"irreversible")==0){
      parameters->irreversible=true;
    }
    else if (strcmp(rhs,"reversible")==0){
      parameters->irreversible=false;
    }
    else {
      fprintf(stderr, "error: 'equation' should be 'irreversible' or 'reversible'\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  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_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);
    }
  }
  else if (strcmp(lhs,"K1")==0){
    ret=sscanf(rhs,"%d",&(parameters->K1));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'K1' should be an integer\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"K2")==0){
    ret=sscanf(rhs,"%d",&(parameters->K2));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'K2' should be an integer\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"K")==0){
    ret=sscanf(rhs,"%d",&(parameters->K1));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'K' should be an integer\n       got '%s'\n",rhs);
      return(-1);
    }
    parameters->K2=parameters->K1;
  }
  else if (strcmp(lhs,"N1")==0){
    ret=sscanf(rhs,"%d",&(parameters->N1));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'N1' should be an integer\n       got '%s'\n",rhs);
      return(-1);
    }
    *setN1=true;
  }
  else if (strcmp(lhs,"N2")==0){
    ret=sscanf(rhs,"%d",&(parameters->N2));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'N2' should be an integer\n       got '%s'\n",rhs);
      return(-1);
    }
    *setN2=true;
  }
  else if (strcmp(lhs,"N")==0){
    ret=sscanf(rhs,"%d",&(parameters->N1));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'N' should be an integer\n       got '%s'\n",rhs);
      return(-1);
    }
    parameters->N2=parameters->N1;
    *setN1=true;
    *setN2=true;
  }
  else if (strcmp(lhs,"final_time")==0){
    ret=sscanf(rhs,"%lf",&(parameters->final_time));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'final_time' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"nu")==0){
    ret=sscanf(rhs,"%lf",&(parameters->nu));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'nu' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"delta")==0){
    ret=sscanf(rhs,"%lf",&(parameters->delta));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'delta' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"L")==0){
    ret=sscanf(rhs,"%lf",&(parameters->L));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'L' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"adaptive_tolerance")==0){
    ret=sscanf(rhs,"%lf",&(parameters->adaptive_tolerance));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'adaptive_tolerance' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"adaptive_factor")==0){
    ret=sscanf(rhs,"%lf",&(parameters->adaptive_factor));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'adaptive_factor' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"max_delta")==0){
    ret=sscanf(rhs,"%lf",&(parameters->max_delta));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'max_delta' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"adaptive_cost")==0){
    if (strcmp(rhs,"L1")==0){
      parameters->adaptive_cost=COST_L1;
    }
    else if (strcmp(rhs,"k3")==0){
      parameters->adaptive_cost=COST_k3;
    }
    else if (strcmp(rhs,"k32")==0){
      parameters->adaptive_cost=COST_k32;
    }
    else if (strcmp(rhs,"enstrophy")==0){
      parameters->adaptive_cost=COST_ENSTROPHY;
    }
    else if (strcmp(rhs,"alpha")==0){
      parameters->adaptive_cost=COST_ALPHA;
    }
    else{
      fprintf(stderr, "error: unrecognized adaptive_cost '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"print_freq")==0){
    ret=sscanf(rhs,"%lf",&(parameters->print_freq));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'print_freq' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"random_seed")==0){
    ret=sscanf(rhs,"%d",&(parameters->seed));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'random_seed' should be an 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);
    }
  }
  else if (strcmp(lhs,"D_epsilon")==0){
    ret=sscanf(rhs,"%lf",&(parameters->D_epsilon));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'D_epsilon' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"driving")==0){
    if (strcmp(rhs,"zero")==0){
      parameters->driving=DRIVING_ZERO;
    }
    else if (strcmp(rhs,"test")==0){
      parameters->driving=DRIVING_TEST;
    }
    // matches any argument that starts with 'file:'
    else if (strncmp(rhs,"file:",5)==0){
      parameters->driving=DRIVING_FILE;
      str_to_dstring_noinit(rhs+5,drivingfile_str);
    }
    // matches any argument that starts with 'file_txt:'
    else if (strncmp(rhs,"file_txt:",9)==0){
      parameters->driving=DRIVING_FILE_TXT;
      str_to_dstring_noinit(rhs+9,drivingfile_str);
    }
    else{
      fprintf(stderr, "error: unrecognized driving force '%s'\n",rhs);
      return(-1);
    }
  }
  // initial condition
  else if (strcmp(lhs,"init")==0){
    if (strcmp(rhs,"random")==0){
      parameters->init=INIT_RANDOM;
    }
    else if (strcmp(rhs,"gaussian")==0){
      parameters->init=INIT_GAUSSIAN;
    }
    // matches any argument that starts with 'file:'
    else if (strncmp(rhs,"file:",5)==0){
      parameters->init=INIT_FILE;
      str_to_dstring_noinit(rhs+5,initfile_str);
    }
    // matches any argument that starts with 'file_txt:'
    else if (strncmp(rhs,"file_txt:",9)==0){
      parameters->init=INIT_FILE_TXT;
      str_to_dstring_noinit(rhs+9,initfile_str);
    }
    else{
      fprintf(stderr, "error: unrecognized initial condition '%s'\n",rhs);
      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 if (strcmp(rhs,"RKF45")==0){
      parameters->algorithm=ALGORITHM_RKF45;
    }
    else if (strcmp(rhs,"RKDP54")==0){
      parameters->algorithm=ALGORITHM_RKDP54;
    }
    else if (strcmp(rhs,"RKBS32")==0){
      parameters->algorithm=ALGORITHM_RKBS32;
    }
    else{
      fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"keep_en_cst")==0){
    int tmp;
    ret=sscanf(rhs,"%d",&tmp);
    if(ret!=1 || (tmp!=0 && tmp!=1)){
      fprintf(stderr, "error: parameter 'keep_en_cst' should be 0 or 1\n       got '%s'\n",rhs);
      return(-1);
    }
    parameters->keep_en_cst=(tmp==1);
  }
  else if (strcmp(lhs,"print_alpha")==0){
    int tmp;
    ret=sscanf(rhs,"%d",&tmp);
    if(ret!=1 || (tmp!=0 && tmp!=1)){
      fprintf(stderr, "error: parameter 'print_alpha' should be 0 or 1\n       got '%s'\n",rhs);
      return(-1);
    }
    parameters->print_alpha=(tmp==1);
  }
  else if (strcmp(lhs,"lyapunov_reset")==0){
    if(parameters->lyapunov_trigger==LYAPUNOV_TRIGGER_SIZE){
      fprintf(stderr, "error: cannot use 'lyapunov_reset' and 'lyapunov_maxu' in the same run:\n  'lyapunov_reset' is to be used to renormalize the tangent flow at fixed times, 'lyapunov_maxu' is to be used to renormalize the tangent flow when the matrix exceeds a certain size.");
      return(-1);
    }
    parameters->lyapunov_trigger=LYAPUNOV_TRIGGER_TIME;
    ret=sscanf(rhs,"%lf",&(parameters->lyapunov_reset));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'lyapunov_reset' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  else if (strcmp(lhs,"lyapunov_maxu")==0){
    if(parameters->lyapunov_trigger==LYAPUNOV_TRIGGER_TIME){
      fprintf(stderr, "error: cannot use 'lyapunov_maxu' and 'lyapunov_reset' in the same run:\n  'lyapunov_reset' is to be used to renormalize the tangent flow at fixed times, 'lyapunov_maxu' is to be used to renormalize the tangent flow when the matrix exceeds a certain size.");
      return(-1);
    }
    parameters->lyapunov_trigger=LYAPUNOV_TRIGGER_SIZE;
    ret=sscanf(rhs,"%lf",&(parameters->lyapunov_reset));
    if(ret!=1){
      fprintf(stderr, "error: parameter 'lyapunov_maxu' should be a double\n       got '%s'\n",rhs);
      return(-1);
    }
  }
  // algorithm for tangent flow (for lyapunov)
  else if (strcmp(lhs,"algorithm_lyapunov")==0){
    if (strcmp(rhs,"RK4")==0){
      parameters->algorithm_lyapunov=ALGORITHM_RK4;
    }
    else if (strcmp(rhs,"RK2")==0){
      parameters->algorithm_lyapunov=ALGORITHM_RK2;
    }
    else if (strcmp(rhs,"RKF45")==0 || strcmp(rhs,"RKDP45")==0 || strcmp(rhs,"RKBS32")==0){
      fprintf(stderr, "error: cannot use an adaptove step algorithm for the tangent flow (Lyapunov exponents); must be one of 'RK4' or 'RK2', got:'%s'\n",rhs);
      return(-1);
    }
    else{
      fprintf(stderr, "error: unrecognized algorithm '%s'\n",rhs);
      return(-1);
    }
  }
  else{
    fprintf(stderr, "error: unrecognized parameter '%s'\n",lhs);
    return(-1);
  }

  return(0);
}

// read args from file
int args_from_file(
  dstring* params,
  unsigned int* command,
  unsigned int* nthreads,
  dstring* savefile_str,
  dstring* utfile_str,
  char* file_str
){
  dstring line;
  char c;

  // open file
  FILE* file=fopen(file_str,"r");
  if(file==NULL){
    fprintf(stderr,"Error opening file '%s' for reading: %s\n", file_str, strerror(errno));
    return(-1);
  }

  // allocate line buffer
  dstring_init(&line,64);

  while(1){
    c=fgetc(file);

    // end of file
    if (feof(file)){
      break;
    }

    // newline: read line and reset buffer
    if(c=='\n' || c=='\r'){
      // read entry
      // ignore lines with fewer than three characters
      if(line.length>=3){
	// find line starting with "#!"
	if(line.string[0]=='#' && line.string[1]=='!'){
	  wordexp_t pwordexp;
	  wordexp(line.string+2,&pwordexp,0);
	  // read arguments
	  read_args(pwordexp.we_wordc, (const char **)pwordexp.we_wordv, params, command, nthreads, savefile_str, utfile_str, NULL);
	  wordfree(&pwordexp);
	  // exit
	  break;
	}
      }
      // reset buffer
      line.length=0;
    }
    // add to buffer
    else {
      dstring_append(c,&line);
    }
  }
  dstring_free(line);
  // close file
  fclose(file);

  return(0);
}

// set driving force
_Complex double* set_driving(
  nstrophy_parameters parameters
){
  _Complex double* g=calloc(parameters.K1*(2*parameters.K2+1)+parameters.K2,sizeof(_Complex double));

  switch(parameters.driving){
  case DRIVING_ZERO:
    g_zero(g, parameters.K1, parameters.K2);
    break;
  case DRIVING_TEST:
    g_test(g, parameters.K1, parameters.K2);
    break;

  case DRIVING_FILE:
    init_file_bin(g, parameters.K1, parameters.K2, parameters.drivingfile);
    break;

  case DRIVING_FILE_TXT:
    init_file_txt(g, parameters.K1, parameters.K2, parameters.drivingfile);
    break;

  default:
    g_test(g, parameters.K1, parameters.K2);
    break;
  }
   return g;
}

// set initial condition
_Complex double* set_init(
  nstrophy_parameters* parameters
){
  _Complex double* u0=calloc(parameters->K1*(2*parameters->K2+1)+parameters->K2,sizeof(_Complex double));

  switch(parameters->init){
  case INIT_RANDOM:
    init_random(u0, parameters->K1, parameters->K2, parameters->seed);
    break;

  case INIT_GAUSSIAN:
    init_gaussian(u0, parameters->K1, parameters->K2);
    break;

  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);
    }
    break;

  case INIT_FILE_TXT:
    init_file_txt(u0, parameters->K1, parameters->K2, parameters->initfile);
    break;

  default:
    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;
}