input.c

波浪数值模拟

/*
 * Copyright (c) 2001--2005 Falk Feddersen
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */


/*  input_new.c -- funcs for loading IC's and init.dat file */ 

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


#include "input.h"

const double pi = 3.14159265358979;  // is needed in input.c for peturbing the forcing - get rid of it eventually
int lineno=0;

void parse_error_message(char *s)
{

  fprintf(stderr,"Error parsing init file at line number: %d\n",lineno);
  fprintf(stderr,"%s\n",s);
  exit(-1);

}


/* --------- input file parsing function -------------
   These should really be in a different file           */


char* get_next_line(FILE *fp, char *s)
{
  char c1 = '%';
  char *rt = NULL;

  while (c1 == '%') {
    lineno++;
    rt=fgets(s,120,fp);
    if (rt==NULL)
      return rt;
    c1 = s[0];
  }

  return rt;
}

/* Process timing routines */


char *units[4]={"(sec)","(min)","(hr)","(day)"};


void timing_report(timing *T)
{

  fprintf(stderr,"------------- * Timing Report * --------------\n");
  fprintf(stderr,"Total Run Time = %9.3g (s)    dt = %7.3g (s)\n",T->total_run_time,T->dt);
  fprintf(stderr,"Level1 Run Time = %9.3g (s)   Level2 Run Time= %9.3g (s)\n",T->level1_time,T->level2_time);
  fprintf(stderr,"Level3 # = %d   Level2 # = %d   Level1 # = %d\n",T->level3_iter,T->level2_iter,T->level1_iter);
  fprintf(stderr,"Default time units are %s\n",T->units_ptr);
  fprintf(stderr,"-----------------------------------------------\n");
}


/* returns index for which units are represented by the input string */

int fix_units(char *unit)
{
  int i;

  for (i=0;i<4;i++)
    if (strcmp(unit,units[i])==0)
      return i;
  
  fprintf(stderr,"Syntax Error: Don't recognize the time units %s.\n",unit);
  exit(-1);
}

/* Convert the units into seconds - the default units */

gdouble change_units(gdouble t, int u)
{
  switch (u) {
  case 0: return t;
  case 1: return 60.0*t;
  case 2: return 3600.0*t;
  case 3: return 24.0*3600.0*t;
  default: parse_error_message("Bad flag to change_units!");
  }

  return -1.0;   // this should never happen!
}

gdouble get_time_in_current_units(timing *T)
{
  gdouble time;

  time = T->time_index * T->dt;

  switch (T->units_flag) {
  case 0: return time;
  case 1: return time/60.0;
  case 2: return time/3600.0;
  case 3: return time/(24.0*3600.0);
  default: parse_error_message("Bad flag to change_units_back!");
  }

  return -1.0;  // this should never happen!
}


gdouble fmodfalk(gdouble num, gdouble denom)
{

   gdouble t1;
   gdouble t2=0.0, t3=0.0;

   t1 = num/denom;
   t2 = modf(t1,&t3);

   return t2;

}

int fdivfalk(gdouble num, gdouble denom)
{

   gdouble t1;
   gdouble t2=0.0, t3=0.0;
   int tt3;

   t1 = num/denom;
   t2 = modf(t1,&t3);

   tt3 = (int ) t3;

   return tt3;

}


timing* process_timing(FILE *fp)
{

  char s[120], tstring[120];
  gdouble time1,time2,time3,dt;
  gdouble t1, t2, t3, dt2;
  int num;

  char units1[12], units2[12], units3[12], units4[12];
  int un1, un2, un3, un4;

  timing *T;

  /* load the time info string  */

  if (get_next_line(fp,s)==NULL) 
    parse_error_message("Syntax Error: Timing information missing from init file.");

  /* allocate the timing */

  T = (timing * ) g_malloc(sizeof(timing));

  T->time_index = 0;
  T->finished = 0;
  T->current_time = 0.0;

  /* load individual component pieces */

  if ((num=sscanf(s,"%s %lf  %s  %lf %s %lf %s %lf %s",tstring, &t1,units1,&t2,units2,&t3,units3,&dt2, units4))!=9) 
    parse_error_message("Syntax Error: Timing information incorrect - expecting 9 - please check.");

  if (strcmp(tstring,"timing")!=0) {
    parse_error_message("Syntax Error: Timing information incorrect - expecting 'timing' as first token.");
  }

  un1 = fix_units(units1);        un2 = fix_units(units2);
  un3 = fix_units(units3);        un4 = fix_units(units4);

  T->units_ptr = units[un1];
  T->units_flag = un1;
 
  time1 = change_units(t1,un1);     time3 = change_units(t3,un3);
  time2 = change_units(t2,un2);     dt = change_units(dt2,un4);

  /* check for AB3 time stepping */

#ifdef AB3
     dt2 = dt;
#else
     dt2 = 2.0*dt;
#endif

   if (fmodfalk(time1,time2)!=0.0) 
     parse_error_message("Syntax Error: Time1 not integral to Time2.  Check init file.");
   else 
     T->level3_iter = fdivfalk(time1,time2);

   if (fmodfalk(time2,time3)!=0.0) 
     parse_error_message("Syntax Error: Time2 not integral to Time3.  Check init file.");
   else 
     T->level2_iter = fdivfalk(time2,time3);
   if (fmodfalk(time3,dt2)!=0.0) 
     parse_error_message("Syntax Error: Time3 not integral to dt.  Check init file.");
   else 
     T->level1_iter = fdivfalk(time3,dt2);

   T->total_run_time = time1;
   T->level1_time = time3;
   T->level2_time = time2;
   T->dt = dt;

   return T;
}


mix_t mix_t_index(char *a)
{
  
  if (!strcasecmp(a,"newtonian0"))
    return NEWT0;

  if (!strcasecmp(a,"newt0"))
    return NEWT0;

  if (!strcasecmp(a,"newtonian1"))
    return NEWT1;

  if (!strcasecmp(a,"newt1"))
    return NEWT1;

  if (!strcasecmp(a,"newtonian2"))
    return NEWT2;

  if (!strcasecmp(a,"newt2"))
    return NEWT2;

  if (!strcasecmp(a,"biharmonic0"))
    return BI0;

  if (!strcasecmp(a,"bi0"))
    return BI0;

  if (!strcasecmp(a,"biharmonic1"))
    return BI1;

  if (!strcasecmp(a,"bi1"))
    return BI1;

  if (!strcasecmp(a,"biharmonic2"))
    return BI2;

  if (!strcasecmp(a,"bi2"))
    return BI2;

  if (!strcasecmp(a,"none"))
    return NONE;

  parse_error_message("** Syntax Error: Wrong lateral mixing type");

}


dynamics_t dynamics_t_index(char *d)
{
  
  if (!strcasecmp(d,"linear"))
    return LINEAR_DYNAMICS;

  if (!strcasecmp(d,"nswe"))
    return NSWE_DYNAMICS;

  if (!strcasecmp(d,"peregrine"))
    return PEREGRINE_DYNAMICS;

  if (!strcasecmp(d,"nwogu"))
    return NWOGU_DYNAMICS;

  if (!strcasecmp(d,"wei_kirby"))
    return WEI_KIRBY_DYNAMICS;

  parse_error_message("** Syntax Error: Wrong dynamics type");

}




void parse_dynamics_info(char *s, inputinfo *I)
{
  char dyn[32], dyntype[32];
  char bub[32];
  int num;

  num = sscanf(s,"%s %s %s",bub,dyn,dyntype);  
  
  if (num != 3) {
    parse_error_message("** Syntax error: not enough arguments  on 'funwaveC dynamics [..]' line **");
  }

  if (strcmp(bub,"funwaveC")!=0) {
    parse_error_message("** Syntax error on 'funwaveC dynamics [..]' line: expecting 'funwaveC' as 1st argument  **");
  }

  if (strcasecmp(dyn,"dynamics")!=0) {
    parse_error_message("** Syntax error on 'funwaveC dynamics [..]' line: expecting 'dynamics' as 2nd argument  **");
  }

  I->Dynamics = dynamics_t_index(dyntype);

  return;

}

/* Depth parsing stuff */

void parse_bathymetry(char *s, inputinfo *I)
{

  char type[32];
  char bub[32];
  char fname[80];

  gdouble d0,d1;
  depthinfo *D;

  D = (depthinfo *) g_malloc(sizeof(depthinfo));


  sscanf(s,"%s %s",bub,type);  

  if (strcmp(bub,"bathymetry")!=0) 
    parse_error_message("** Syntax error on bathymetry line **");

  if (strcasecmp(type,"flat")==0) {
    sscanf(s,"%*s %*s %lf",&d0);
    D->dtype = FLAT;
    D->h0 = d0;
    I->D = D;
    return;
  }

  if (strcasecmp(type,"planar")==0) {
    sscanf(s,"%*s %*s %lf %lf",&d0,&d1);
    D->dtype = PLANAR;
    D->h0 = d0;
    D->slope = d1;
    I->D = D;
    return;
  }

  if (strcasecmp(type,"file1D")==0.0) {
    sscanf(s,"%*s %*s %s",fname);
    D->dtype = HFILE1D;
    D->depth_file = g_string_new(fname);
    I->D = D;
    return;
  }

  if (strcasecmp(type,"file2D")==0.0) {
    sscanf(s,"%*s %*s %s",fname);
    D->dtype = HFILE2D;
    D->depth_file = g_string_new(fname);
    I->D = D;
    return;
  }

  parse_error_message("Syntax Error: Unknown Type for Bathymetry (flat, planar, file1D, file2D).");

}



void  parse_dimensions(char *s, inputinfo *I)
{

  char bub[120];
  int num, nx, ny;
  double dx, dy;

  num = sscanf(s,"%s %d %d %lf  %lf",bub,&nx, &ny, &dx, &dy);

  if (strcasecmp(bub,"dimension")) {
    parse_error_message("Bad first line for dimension input");
  }

  if (num!=5)
    parse_error_message("Not enough inputs for <dimensions nx ny dx dy> input.");


  I->nx = nx;   I->ny = ny;
  I->dx = dx;   I->dy = dy;

}


void  parse_bottom_stress(char *s, inputinfo *I)
{

  char bub[120];
  int num; 
  double cd; 

  num = sscanf(s,"%s %lf ",bub, &cd);

  if (strcasecmp(bub,"bottomstress")) {
    parse_error_message("Bad first line for bottom stress input");
  }

  if ((num <2)||(num >3))
    parse_error_message("Not enough inputs for bottom stress input");

  /* next check to make sure coefficients they are > 0 */

  if (cd<0.0) {
    printf("** shallow.c: cd<0, cd=%lf\n",cd);
    parse_error_message("Specified value for cd < 0.0");
  }

  I->cd = cd;


}



void  parse_lateral_mixing(char *s, inputinfo *I)
{

  char bub[120];
  char diss_str[120];
  int num; 
  double viscosity;
  mix_t DT;

  num = sscanf(s,"%s %s %lf",bub, diss_str,  &viscosity);

  if (strcasecmp(bub,"mixing")) {
    parse_error_message("Bad first line for lateral mixing input");
  }

  if ((num <3)||(num >3))
    parse_error_message("Not enough inputs for lateral mixing: 'lateral diss_str viscosity'");

  /* next check to make sure coefficients they are > 0 */

  DT = mix_t_index(diss_str);
  I->DT = DT;

  if (viscosity<0) {
    fprintf(stderr,"** shallow.c: viscosity=%lf\n",viscosity);  
    parse_error_message("Specified value for newtonian or hyper viscosity < 0");
  }

  
  switch (DT) {
  case NEWT0:
  case NEWT1:
  case NEWT2:  I->nu_newt = viscosity;  I->nu_bi = 0.0; break;
  case BI0:
  case BI1:
  case BI2:    I->nu_bi = -1.0*viscosity;  I->nu_newt = 0.0;  break;
  case NONE:   I->nu_bi = 0.0;   I->nu_newt = 0.0;  break;
  }
  

}


void  parse_eta_source(char *s, inputinfo *I)
{

  char bub[120];
  char e_str[120], on_off[120], filename[120];
  int num; 
  double  freq, H, theta, delta, fwidth;
  int nf;
  int icenter;
  eta_source_info_t *ES;


  /* this input line should look like:
     "eta_source [on,off] [monochromatic,random]  H (m)  freq (Hz)  theta (deg),  icenter "  */

  num = sscanf(s,"%s %s %s %lf %lf %lf  %d  %lf",bub, on_off, e_str, &H, &freq, &theta, &icenter, &delta);

  if (strcasecmp(bub,"eta_source")) {
    parse_error_message("Bad first line for eta_source input");
  }

  if ((num <2)||(num >8))
    parse_error_message("Not enough inputs for eta_source: 'eta_source [on,off] [monochromatic,random] H f theta icenter'");

  /* next check to make sure coefficients they are > 0 */

  ES = (eta_source_info_t * ) g_malloc( sizeof(eta_source_info_t) );

  /* check if eta source is on or off */

  if (!strcasecmp(on_off,"on")) {
    if (num<9) 
      ES->eta_source_flag = ON;
    else
      parse_error_message("** Syntax Error: Missing information arguments on eta_source input line");
  }
  else {
    if (!strcasecmp(on_off,"off")) {
      ES->eta_source_flag = OFF;
      I->ES = ES;
      return;
    }
    else
      parse_error_message("** Syntax Error: Bad 2nd option for eta_source input line, should be [on,off]");
  }

  /* check if eta source is monochromatic or random */

  if (!strcasecmp(e_str,"monochromatic")) {