eta_source.c

波浪数值模拟

    RW1->phase[i] = random_phase();  
    RW1->ktot[i] = get_wavenum(RW1->omega[i], RW1->h0);
    RW1->theta[i] = calculate_fit_wave_angle(RW1->ktot[i], RW1->theta[i], ny, &kynew, &rw1_n_wavenumber);
    RW1->ky[i] = kynew;

    sqramp += aa*aa;
    thetasum += aa*aa*(th*pi/180.0);
    freqsum +=  aa*aa*ff;
  }

  RW1->freq0 = freqsum/sqramp;
  RW1->theta0 = thetasum/sqramp;
  RW1->Hsig =  4.0*sqrt(sqramp);  /* check if this is 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->fwidth = -1.0;
  RW1->omega0 = 2.0*pi*RW1->freq0;
  
  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

  /* Now normalize the RW1->amplitudes so that their square sum = 1 */

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

  fclose(fp);

}




void init_random1_wave_file(eta_source_info_t *S, double h0)
{


  int i, j, k;
  int rw1_n_wavenumber;
  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_FILE;
  RW1->numfreq = count_line_numbers(S->spectrum_file->str);
  RW1->h0 = h0;


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

  init_random_phase_generator();   /* need to initialize the random phases first */

  /* Now load up  omega, amplitude, theta */

  random1_load_spectrum_file(RW1, S->spectrum_file->str);  // this should set Hsig, theta0, freq0 etc. */


  /* 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 );
    }
  }


}



void  rhs_calc_eta_source_f_random1_add(field2D *ETAT)
{


  double ggx;
  double rfsrc;

  //  double dwt, dw, wt;

  int M = ETAT->M;
  int i,j,k;

  double *etatd = REF2(ETAT, ixstart);

  double *romega = RW1->omega;
  double *ramp = RW1->amplitude;
  double *rphase = RW1->phase;
  double *rsinwt = RW1->rsinwt;
  double *rcoswt = RW1->rcoswt;
  double *g_x = RW1->g_x;

  double *rs_ytA = REF2(RW1->s_ytA, 0);
  double *rs_ytB = REF2(RW1->s_ytB, 0);

  for (k=0;k < RW1->numfreq; k++) {

   /* note addition of phases here - = A(k) * sin[ omega(k)*t + phase(t) ] */

    rsinwt[k] =  ramp[k]*sin(romega[k]*local_time+rphase[k]);   
    rcoswt[k] =  ramp[k]*cos(romega[k]*local_time+rphase[k]);
  }

  for (j=0;j<M;j++) {
    rfsrc = 0.0;
    for (k=0; k< RW1->numfreq; k++) {   /* allows for 1 wave angle at each freq */
      rfsrc += ( (*rs_ytA++)  * rcoswt[k] - (*rs_ytB++) * rsinwt[k]);
    }

    etatd = REF2(ETAT, ixstart) + j;  /* set up the d(eta)/dt pointer at the right j */

    for (i=0;i<num_gx;i++) {
      *(etatd+i*M) += rfsrc * g_x[i];
    }

  }

  local_time += local_dt;

}


#define NUMTHETA 7    /* default number of ky at each frequency bin */

typedef struct {
  random2_subtype  subtype;   /* this is the subtype of the random1 struct */

  int numfreq;     /* the total number of frequencies to include */
  int numky;       /* the total number of ky at each frequency */
  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;   /* the  average over freq and direction wave angle */
  double spread0;   /* the  average over freq and direction wave angle */
  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 *amplitude;   /* the actual 2*pi*f array of radian frequencies */
  double *mtheta;   /* the angle at each frequency in radians */
  double *spread;   /* the angle at each frequency in radians */
  double *ktot;     /* the total wavenumber at each frequency */

  field2D *kamp;       /* the (omega,ky) amplitudes:  size (numfreq,NUMTHETA)   */
  field2D *ky;      /* the alongshore wavenumber at each frequency size (numfreq,NUMTHETA) */

  double *g_x;

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

} random2_wave_t;

random2_wave_t *RW2;


/* 
   This function takes in the set of ky and kamp and calculates s_ytA and s_ytB:
   s_ytA(y,freq) = \Sum_j  d_{ij} cos[ k_{ij}y + phase ]  : j = wavenumber index, i = freq index
   s_ytB(y,freq) = \Sum_j  d_{ij} sin[ k_{ij}y + phase ]  : j = wavenumber index, i = freq index 
     - one of the .init functions needs to call this
     - assumes that the field2D of  ky and kamp are already set up and allocated

*/


void random2_calculate_s_ytAB(random2_wave_t *R2)
{

  int i, j, l;
  double stmp1, stmp2;
  double y;
  double dij, kky, ph;

  field2D *sytA;
  field2D *sytB;

  sytA = allocate_field2D(ny, R2->numfreq);
  sytB = allocate_field2D(ny, R2->numfreq);

  for (l=0;l<ny;l++) {
    y = l*dy;
    for (i=0;i<R2->numfreq;i++) {
      stmp1 = 0.0;
      stmp2 = 0.0;
      for (j=0;j<R2->numky;j++) {
	dij = DR2(R2->kamp, i, j);
        kky = DR2(R2->ky, i, j);
	ph = random_phase();
	stmp1 += dij * cos( kky*y +ph);
	stmp2 += dij * sin( kky*y +ph);
      }

      DR2( sytA, l, i) = stmp1;    //D * sin( (RW1->ky[k])*y );
      DR2( sytB, l, i) = stmp2;   //D * cos( (RW1->ky[k])*y );
    }
  }

  R2->s_ytA = sytA;
  R2->s_ytB = sytB;


}


/* 
   Function takes in mean angle and spread (in radians) and # of total ky and returns the ky and the amplitudes 
   normalized to sum = 1                
*/

void calculate_ky_from_dir_spread_gaussian(double mtheta, double spread, double ktot, int nky, double *kampd, double *kyd)
{

  int i;
  int n_wavenumber;
  double kynew;
  double theta, theta_new;
  double dtheta = 3.0*spread/(nky-1);
  int nkyhalf = nky/2;
  double tot = 0.0;

  for (i=0;i<nky;i++) {
    theta = mtheta + (i-nkyhalf)*dtheta;
    theta_new = calculate_fit_wave_angle(ktot, theta, ny, &kynew, &n_wavenumber);
    kyd[i] = kynew;
    kampd[i] = exp(  -0.5 * SQR( (theta_new-mtheta)/spread ) );
    tot += kampd[i];
  }

  /* Now normalize so that the sum of kampd = 1 */

  tot /= nky;
  for (i=0;i<nky;i++) {
    kampd[i] /= tot;
  }

}





void init_random2_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;

  RW2 = (random2_wave_t * ) g_malloc(sizeof(random2_wave_t));

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

  RW2->omega = (double * ) g_malloc(sizeof(double)*S->numfreq);
  //  RW2->phase = (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW2->amplitude = (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW2->mtheta = (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW2->ktot =    (double * ) g_malloc(sizeof(double)*S->numfreq);

  RW2->ky =    allocate_field2D(S->numfreq, RW2->numky) ;
  RW2->kamp =    allocate_field2D(S->numfreq, RW2->numky) ;

  RW2->rsinwt =    (double * ) g_malloc(sizeof(double)*S->numfreq);
  RW2->rcoswt =    (double * ) g_malloc(sizeof(double)*S->numfreq);

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

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

  for (i=0;i<RW2->numfreq;i++) {
    freq = (RW2->freq0) + df*((i-halfnumfreq));
    RW2->omega[i] = 2*pi*freq;
    //    RW2->phase[i] = random_phase();  
    RW2->ktot[i] = get_wavenum(RW2->omega[i], h0);
    RW2->mtheta[i] = RW2->theta0;
    RW2->spread[i] = RW2->spread0;
   
    //    calculate_ky_from_dir_spread_gaussian(double mtheta, double spread, double ktot, int nky, double *kampd, double *kyd)
    //    RW2->ky[i] = kynew; FIX!!!
    RW2->amplitude[i] = exp( -0.5 * SQR( (freq-(RW2->freq0))/freqhalfwidth ) );  /* needs to be normalized so sum squares =1 */
    sqramp += SQR( RW2->amplitude[i] );
  }

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

  RW2->ktot0 = get_wavenum(RW2->omega0, h0);  // h0 should be the offshore depth
  RW2->kh0 = RW2->ktot0 * h0;               // nominal value of kh0 at the wavemaker
  RW2->L = 2*pi/(RW2->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 = (RW2->ktot0) * cos(RW2->theta0);   // peak wavenumber in x direction
  W = delta * (RW2->L) /2.0;     // this is the width of the source function,
                                //  but not the L not Lx = 2*pi/k_x !
  tmp = SQR( delta*(RW2->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*(RW2->a0) * cos((RW2->theta0)) * ( SQR(RW2->omega0) - alpha1*gravity* SQR( ((RW2->ktot0)*h0) ) * SQR( (RW2->ktot0) ) * h0 );
  D /= (RW2->omega0 * (RW2->ktot0) * I ) * (1.0 - alpha* SQR( ((RW2->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;

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

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

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

  for (j=0;j<ny;j++) {
    y = j*dy;
    for (k=0;k<RW2->numfreq;k++) {
      //      DR2( RW2->s_ytA, j, k) = D * sin( (RW2->ky[k])*y );
      //      DR2( RW2->s_ytB, j, k) = D * cos( (RW2->ky[k])*y );  NEED TO FIX THIS BADLY
    }
  }


}



/* 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 RW2 and also Hsig, theta0, freq0 etc. are calculated */


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

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

  /* Read in and process the random2 spectrum file : the file format is
        freq (hz)  amp (?)  theta(deg)  spread(deg)  - assumes a cos^(2n) type theta distribution?? */

  for (i=0; i<R2->numfreq; i++) {
    num = fscanf(fp,"%lf %lf %lf %lf",&ff, &aa, &th, &sp);
    if (num!=4) {
      funwaveC_error_message("Random2 spectrum file did not read 4 elements");
    }
    R2->omega[i] = 2.0*pi*ff;
    R2->amplitude[i] = aa;
    R2->mtheta[i] = (pi/180.0)*th;   // convert angle to radians
    R2->spread[i] = (pi/180.0)*sp;   // convert spread to radians
    R2->ktot[i] = get_wavenum(R2->omega[i], R2->h0);

    /* do this part of calculating the angles later */


    //    R2->ky[i] = kynew;

    sqramp += aa*aa;
    thetasum += aa*aa*(th*pi/180.0);
    freqsum +=  aa*aa*ff;
  }


  R2->freq0 = freqsum/sqramp;
  R2->theta0 = thetasum/sqramp;
  R2->Hsig =  4.0*sqrt(sqramp);  /* check if this is right! */

  R2->a0 = (R2->Hsig)/2.0;  /* is this also right?? or should it be Hrms ? */
  R2->aih0 = R2->a0/h0;                                    /* nominal value of a/h at the wavemaker */