eta_source.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
 *
 */



/* eta_source.c   Falk Feddersen */       

#include <time.h>
#include <stdlib.h>
#include <math.h>

#include "eta_source.h"


/* What needs to be done here?  
    1) figure out the wave length L of the wave with the given frequency
    2) figure out where in x the source function is centered
    3) figure out the width of the source function:  say L/2
    4) figure out \beta from Kirby's thing
    5) Figure out the amplitude of the source function
*/

// Make these global variables to eta_source orginally

switch_t eta_source_flag;

wave_t wave_type;


typedef struct {
  double h0;         /* the average wavemaker water depth */
  double freq;       /* the peak frequency in Hz*/
  double omega;      /* the actual 2*pi*f array of radian frequencies */
  double Hwave, a0;  /* wave height and amplitude */
  double theta;      /* wave angle in radians */
  double ktot;       /* the total wavenumber */
  double kh0;        /* non-d dispersion parameter */
  double aih0;       /* a/h paramter */
  int    n_wavenumber; /* ratio Ly/ly */
  double ky;         /* the alongshore wavenumber at each frequency */
  double L;          /* The total wavelength */
  double *g_x;       /*  the g(x) function g(x) = exp(-beta (x-x_s)^2 ) */
  double *s_ytA, *s_ytB;     /* the alongshore part s(y) as the sin (A) and cos(B) parts  */
} monochromatic_wave_t;

monochromatic_wave_t *MW;


//double L;  // wavelength
//double omega;  // radian frequency
//double freq;   // frequency


//double h0;  // the water depth at wave maker
//double Hwave, a0;  // wave height and amplitude 
//double theta;      // wave angle in radians
//double k_tot;   // total wavenum
//double kh0;     // the value of kh0 at the wavemaker
//double aih0;    // the value of a/h at the wavemaker
 

//double k_y;    //  y-component of wavenumber
//int n_wavenumber;

double delta;   // this relates W = \delta L/2 (p. 16 funwave manual eq 37)
int ixcenter, ixstart;   // the indices of the center and stop of the source function
int num_gx;                // the number of points in the g(x) function
int icgx;                  // this is the center of num_gx in local coord  ic = (num_gx-1)/2
      
//double *g_x;               // the g(x) function g(x) = exp(-beta (x-x_s)^2 )


double local_time = 0.0;
double local_dt;

/* Local function declarations */

double get_cg(double k, double depth);
double get_wavenum(double omega, double depth);
int    k_round_odd(double a);    // returns the nearest odd integer to the double a
double get_average_wavemaker_depth(const field2D *DD, int ixcenter);
double calculate_fit_wave_angle(double k_tot, double theta, int ny, double *ky_new, int *n_wavenumber);
double init_random_phase_generator();
double random_phase();       /* returns a random phase between 0...2*pi */

double round( double x);

void init_monochromatic_wave_info(eta_source_info_t *S, double h0);
void init_random1_wave_narrow_band(eta_source_info_t *S, double h0);
void init_random1_wave_file(eta_source_info_t *S, double h0);

void  init_eta_source(eta_source_info_t *S, int nx, int ny, double dt)
{

  double h0;

  local_dt = dt;

  eta_source_flag = S->eta_source_flag;

  if (eta_source_flag == OFF)  /* if eta source is off then return */
    return;


  wave_type = S->wave_type;
  ixcenter = S->icenter;

  h0 = get_average_wavemaker_depth(H, ixcenter);  /* get the average depth out at the wavemaker */

  if (S->delta == -1.0)   /* here the width parameter of the source region is set */
    delta = 1.0;          /* delta=1.0 means source region is as wide as wavelength */
  else
    delta = S->delta;

  switch (wave_type) {
  case MONOCHROMATIC:  init_monochromatic_wave_info(S, h0 );
    break;
  case RANDOM1: 
    switch (S->subtype1) {
    case RANDOM1_NARROW_BAND:
      init_random1_wave_narrow_band(S, h0);
      break;
    case RANDOM1_FILE:
      init_random1_wave_file(S, h0);
      break;
    }
    break;
  case RANDOM2:
    funwaveC_error_message("RANDOM2 not yet implemented");
    switch (S->subtype2) {
    case RANDOM2_NARROW_BAND:
    case RANDOM2_FILE_A:
      break;
    }
  }


}

void init_monochromatic_wave_info(eta_source_info_t *S, double h0)
{

  int i,j;
  double x, y;
  double tmp;

  double D;     // the magnitude of the source function  (eq 34, p. 15, funwave manual)
  double I;     // the integral I (eq 35, p. 15, funwave manual)

  double W;  // width of source function
  double xc;  // x location of source function
  double beta2;  // FUNWAVE eta source shape paramter
  double k_x;   // x - component of wavenumber 
  const double alpha = -0.390;
  const double alpha1 = alpha+1.0/3.0;


  MW = (monochromatic_wave_t * ) g_malloc(sizeof(monochromatic_wave_t));

  MW->h0 = h0;
  MW->freq = S->frequency;
  MW->omega = MW->freq * 2 *pi;
  MW->Hwave = S->H;                    // the wave height
  MW->a0 = MW->Hwave/2.0;                  // amplitude of the wave
  MW->theta = pi * (S->theta)/180.0;   // convert this to radians 

  MW->aih0 = MW->a0/h0;                                    /* nominal value of a/h at the wavemaker */

  MW->ktot = get_wavenum(MW->omega,h0);  // h0 should be the offshore depth
  MW->kh0 = MW->ktot * h0;               // nominal value of kh0 at the wavemaker
  MW->L = 2*pi/(MW->ktot);                 // this is the overall wavelength

  /* Funtion to ensure that wave angle actually fits the domain! 
     returns wave angle in radians and also the y-wavenumber k_y  */

  MW->theta = calculate_fit_wave_angle(MW->ktot, MW->theta, ny, &(MW->ky), &(MW->n_wavenumber));  

  k_x = (MW->ktot) * cos((MW->theta));   // wavenumber in x direction


  W = delta * (MW->L) /2.0;     // this is the width of the source function,
                          //  but not the L not Lx = 2*pi/k_x !
  tmp = SQR( delta*(MW->L) ) ;
  beta2 = 80.0/ tmp;

  I = sqrt( pi/beta2 ) * exp( - SQR(k_x) / (4.0*beta2) );  // l = k_tot*cos(theta) (eq 35, p. 15, funwave manual)

  D = 2.0*(MW->a0) * cos((MW->theta)) * ( SQR(MW->omega) - alpha1*gravity* SQR( ((MW->ktot)*h0) ) * SQR( (MW->ktot) ) * h0 );
  D /= (MW->omega * (MW->ktot) * I ) * (1.0 - alpha* SQR( ((MW->ktot)*h0) ) );       // eq 34 in funwave manual


  /* Make the g(x0 function terms */
  /* num_gx is odd,  ic = (num_gx-1)/2, so if num_gx=3, ic = 1 etc. */

  num_gx = k_round_odd(W/dx);
  icgx = (num_gx - 1)/2;

  MW->g_x = (double * ) g_malloc( num_gx * sizeof(double) );
  for (i=0;i<num_gx;i++) {
    x = (i -icgx) * dx;
    MW->g_x[i] = exp( -beta2 * x * x );
  }
  ixstart = ixcenter - icgx;

  /* this is for the alongshore component - allocate but don't set up */

  MW->s_ytA = (double *) g_malloc( ny * sizeof(double) );  
  MW->s_ytB = (double *) g_malloc( ny * sizeof(double) );  

  for (j=0;j<ny;j++) {
    y = j*dy;
    MW->s_ytA[j] = D * sin( (MW->ky)*y );
    MW->s_ytB[j] = D * cos( (MW->ky)*y );
  }

}



void  rhs_calc_eta_source_f_monochromatic_add(field2D *ET)
{


  double sinwt, coswt;
  double ggx;
  double fsrc;

  int M = ET->M;
  int i,j;

  double *etatd = &(ET->data[ixstart*M]);
  double *g_x = MW->g_x;
  double *s_ytA = MW->s_ytA;
  double *s_ytB = MW->s_ytB;
  double omega = MW->omega;

  /* Check to see if this is turned on or not */

  if (eta_source_flag == OFF)
    return;  
  
  sinwt =  sin(omega*local_time);
  coswt =  cos(omega*local_time);


  for (i=0;i<num_gx;i++) {
    ggx = g_x[i];
    for (j=0;j<M;j++) {
      fsrc = ggx * (s_ytA[j]*coswt - s_ytB[j]*sinwt);
      *etatd++ += fsrc;
    }
  }

  local_time += local_dt;

}



typedef struct {
  random1_subtype  subtype;   /* this is the subtype of the random1 struct */
  int numfreq;     /* the total number of frequencies to include */
  double h0;       /* the average wavemaker water depth */
  double Hsig;     /* the significant wave height */
  double a0;  
  double freq0;    /* the peak frequency */
  double fwidth;
  double omega0;   /* the peak radian frequency */
  double theta0;   
  double ktot0;       /* the total wavenumber at peak frequency */
  double kh0;        /* non-d dispersion parameter */
  double aih0;       /* a/h paramter */
  double n_wavenumber; /* ratio Ly/ly */
  double ky0;         /* the alongshore wavenumber at each frequency */
  double L;          /* The total wavelength */
  double df;         /* ???? */

  double *omega;   /* the actual 2*pi*f array of radian frequencies */
  double *phase;   /* the random phase at each frequency */
  double *amplitude; /* the amplitude at this frequency */
  double *theta;   /* the angle at each frequency in radians */
  double *ktot;    /* the total wavenumber */
  double *ky;      /* the alongshore wavenumber at each frequency */
  double *g_x;

  double *rsinwt,  *rcoswt;   /* temporary variables used in ..._f_add() */
  field2D *s_ytA;
  field2D *s_ytB;

} random1_wave_t;

random1_wave_t *RW1;


void init_random1_wave_narrow_band(eta_source_info_t *S, double h0)
{


  int i, j, k;
  int rw1_n_wavenumber;
  int halfnumfreq;
  double freq, freqhalfwidth;
  double kynew;
  double k_x, W, tmp, beta2, I, D;
  double x, y;
  double sqramp = 0.0;
  double df;

  const double alpha = -0.390;
  const double alpha1 = alpha+1.0/3.0;

  RW1 = (random1_wave_t * ) g_malloc(sizeof(random1_wave_t));

  RW1->subtype = RANDOM1_NARROW_BAND;
  RW1->numfreq = S->numfreq;
  RW1->h0 = h0;
  RW1->Hsig = S->H;    /* is this right ?? */
  RW1->a0 = (RW1->Hsig)/2.0;  /* is this also right?? or should it be Hrms ? */
  RW1->aih0 = RW1->a0/h0;                                    /* nominal value of a/h at the wavemaker */
  RW1->freq0 = S->frequency;
  RW1->fwidth = S->fwidth;
  RW1->omega0 = 2.0*pi*RW1->freq0;
  RW1->fwidth = S->fwidth;
  RW1->theta0 = (pi/180.0)*S->theta;  /* convert from deg to radians */

  RW1->omega = (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW1->phase = (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW1->amplitude = (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW1->theta = (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW1->ktot =    (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW1->ky =    (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW1->rsinwt =    (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW1->rcoswt =    (double * ) g_malloc(sizeof(double)*S->numfreq);

  df = S->fwidth/(S->numfreq-1);
  halfnumfreq = RW1->numfreq/2;

  freqhalfwidth = S->fwidth/2.0;
  init_random_phase_generator();   /* need to initialize the random phases first */

  for (i=0;i<RW1->numfreq;i++) {
    freq = (RW1->freq0) + df*((i-halfnumfreq));
    RW1->omega[i] = 2*pi*freq;
    RW1->phase[i] = random_phase();  
    RW1->ktot[i] = get_wavenum(RW1->omega[i], h0);
    RW1->theta[i] = calculate_fit_wave_angle(RW1->ktot[i], RW1->theta0, ny, &kynew, &rw1_n_wavenumber);
    RW1->ky[i] = kynew;
    RW1->amplitude[i] = exp( -0.5 * SQR( (freq-(RW1->freq0))/freqhalfwidth ) );  /* needs to be normalized so sum squares =1 */
    sqramp += SQR( RW1->amplitude[i] );
  }

  for (i=0;i<RW1->numfreq; i++) {
    RW1->amplitude[i] /= sqrt( sqramp );
  }

  RW1->ktot0 = get_wavenum(RW1->omega0, h0);  // h0 should be the offshore depth
  RW1->kh0 = RW1->ktot0 * h0;               // nominal value of kh0 at the wavemaker
  RW1->L = 2*pi/(RW1->ktot0);                 // this is the overall wavelength

  /* Funtion to ensure that wave angle actually fits the domain! 
     returns wave angle in radians and also the y-wavenumber k_y  */

  k_x = (RW1->ktot0) * cos(RW1->theta0);   // peak wavenumber in x direction
  W = delta * (RW1->L) /2.0;     // this is the width of the source function,
                                //  but not the L not Lx = 2*pi/k_x !
  tmp = SQR( delta*(RW1->L) ) ;
  beta2 = 80.0/ tmp;

  I = sqrt( pi/beta2 ) * exp( - SQR(k_x) / (4.0*beta2) );  // l = k_tot*cos(theta) (eq 35, p. 15, funwave manual)

  D = 2.0*(RW1->a0) * cos((RW1->theta0)) * ( SQR(RW1->omega0) - alpha1*gravity* SQR( ((RW1->ktot0)*h0) ) * SQR( (RW1->ktot0) ) * h0 );
  D /= (RW1->omega0 * (RW1->ktot0) * I ) * (1.0 - alpha* SQR( ((RW1->ktot0)*h0) ) );       // eq 34 in funwave manual


  /* Make the g(x0 function terms */
  /* num_gx is odd,  ic = (num_gx-1)/2, so if num_gx=3, ic = 1 etc. */

  num_gx = k_round_odd(W/dx);
  icgx = (num_gx - 1)/2;

  RW1->g_x = (double * ) g_malloc( num_gx * sizeof(double) );
  for (i=0;i<num_gx;i++) {
    x = (i -icgx) * dx;
    RW1->g_x[i] = exp( -beta2 * x * x );
  }
  ixstart = ixcenter - icgx;

  /* this is for the alongshore component - allocate but don't set up */

  RW1->s_ytA = allocate_field2D(ny, RW1->numfreq);
  RW1->s_ytB = allocate_field2D(ny, RW1->numfreq);

  for (j=0;j<ny;j++) {
    y = j*dy;
    for (k=0;k<RW1->numfreq;k++) {
      DR2( RW1->s_ytA, j, k) = D * sin( (RW1->ky[k])*y );
      DR2( RW1->s_ytB, j, k) = D * cos( (RW1->ky[k])*y );
    }
  }


}

/* Counts the numbers of lines in a file - needed to return the number of frequencies in a file */

int count_line_numbers(char *filename)
{

  FILE *fp;
  int linenum = 0;
  char s[120];
  char *rt = NULL;

  if ((fp=fopen(s,"r"))==NULL) {
    fprintf(stderr,"Cannot open file for reading: %s\n",s);
    exit(-1);
  }

  while ((rt = fgets(s,120,fp))!=NULL) 
    linenum++;

  fclose(fp);
  return linenum;

}

/* With the correct numfreq, loads up the random1 wave spectrum file that looks like:
   columns of  freq (Hz) ampl (m???)  theta (deg)       
   OUTPUT: these quantities are now in RW1 and also Hsig, theta0, freq0 etc. are calculated */

void   random1_load_spectrum_file(random1_wave_t *RW1, char *spectrumfile)   // this should set Hsig, theta0, freq0 etc.
{
  FILE *fp;
  int i, num;
  double ff, aa, th;
  double sqramp = 0.0, thetasum=0.0, freqsum=0.0;
  double kynew;
  int rw1_n_wavenumber;
  double h0 = RW1->h0;

  if ((fp=fopen(spectrumfile,"r"))==NULL) {
    fprintf(stderr,"Cannot open file for reading: %s\n",spectrumfile);
    exit(-1);
  }

  /* Read in and process the spectrum file */

  for (i=0; i<RW1->numfreq; i++) {
    num = fscanf(fp,"%lf %lf %lf",&ff, &aa, &th);
    if (num!=3) {
      funwaveC_error_message("Random1 spectrum file did not read 3 elements");
    }
    RW1->omega[i] = 2.0*pi*ff;
    RW1->amplitude[i] = aa;
    RW1->theta[i] = (pi/180.0)*th;