eta_source.c

波浪数值模拟

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

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

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

  fclose(fp);

}





/* This function is adds on the directionally spread random wave forcing to d(ETA)/dt
     - it is essentially the same as the .._random1_add() function as all the directional
     stuff is wrapped up in the s_ytA s_ytB stuff                                */

void  rhs_calc_eta_source_f_random2_add(field2D *ETAT)
{


  double ggx;
  double rfsrc;

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

  double *etatd = REF2(ETAT, ixstart);

  double *romega = RW2->omega;
  double *ramp = RW2->amplitude;
  double *rsinwt = RW2->rsinwt;
  double *rcoswt = RW2->rcoswt;
  double *g_x = RW2->g_x;

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

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

   /* note no addition of phases here - = A(k) * sin[ omega(k)*t ] - random phases are in the ky component*/

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

  for (j=0;j<M;j++) {
    rfsrc = 0.0;
    for (k=0; k< RW2->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;

}




void  rhs_calc_eta_source_f_add(field2D *ETAT)
{

  if (eta_source_flag == OFF)
    return;  

  switch(wave_type) {
  case MONOCHROMATIC: 
    rhs_calc_eta_source_f_monochromatic_add(ETAT);
    break;
  case RANDOM1: 
    rhs_calc_eta_source_f_random1_add(ETAT);
    break;
  case RANDOM2: 
    rhs_calc_eta_source_f_random2_add(ETAT);
    break;
  }

}


void eta_source_report_monochromatic(monochromatic_wave_t *M)
{



  fprintf(stderr,"----------- Eta Source Function Report: MONOCHROMATIC Waves  -----------------\n");
  fprintf(stderr,"H = %.2lf (m),  freq = %.2lf (Hz),  Lwave = %.1lf (m) theta = %.3lf (deg)\n",M->Hwave, M->freq, M->L, M->theta*180/pi);
  fprintf(stderr,"h0 = %.3lf (m), kh_0 = %.2lf,   a/h = %.2lf,  LY/ly=%d , k_y = %.5lf (rad/m)\n",M->h0, M->kh0, M->aih0, M->n_wavenumber, M->ky);
  fprintf(stderr,"i_src = %d, i_width =%d,  x_src = %.1lf (m), x_width = %.1lf (m), delta=%.1f\n",ixcenter, num_gx, ixcenter*dx + 0.5*dx, num_gx*dx,delta);
  fprintf(stderr,"-------------------------------------------------------------------------\n");
}


char *random1_str[2] = {"narrow band", "file"};

void eta_source_report_random1(random1_wave_t *R1)
{

  char *str;

  switch (RW1->subtype) {
  case RANDOM1_NARROW_BAND:  
    str = random1_str[0];
    break;
  case RANDOM1_FILE:  
    str = random1_str[1];
    break;
  }

  /* add numfreq and fwidth here!  - add to R1 type */

  fprintf(stderr,"----------- Eta Source Function Report: RANDOM1 Waves; Subtype = %s  -----------------\n",str);
  fprintf(stderr,"Hsig = %.2lf (m),  freq0 = %.2lf (Hz),  fwidth =%.2lf (Hz), Lwave = %.1lf (m) theta0 = %.3lf (deg)\n",
                  R1->Hsig, R1->freq0, R1->fwidth, R1->L, R1->theta0*180/pi);
  fprintf(stderr,"h0 = %.3lf (m), kh_0 = %.2lf,   a/h = %.2lf,  LY/ly=%d , k_y = %.5lf (rad/m)\n",R1->h0, R1->kh0, R1->aih0, R1->n_wavenumber, R1->ky);
  fprintf(stderr,"i_src = %d, i_width =%d,  x_src = %.1lf (m), x_width = %.1lf (m), delta=%.1f\n",ixcenter, num_gx, ixcenter*dx + 0.5*dx, num_gx*dx,delta);
  fprintf(stderr,"-------------------------------------------------------------------------\n");
}



char *random2_str[2] = {"narrow band", "file_a"};

void eta_source_report_random2(random2_wave_t *R2)
{

  char *str;

  switch (R2->subtype) {
  case RANDOM2_NARROW_BAND:  
    str = random2_str[0];
    break;
  case RANDOM2_FILE_A:  
    str = random2_str[1];
    break;
  }

  /* add numfreq and fwidth here!  - add to R2 type */

  fprintf(stderr,"----------- Eta Source Function Report: RANDOM2 Waves; Subtype = %s  -----------------\n",str);
  fprintf(stderr,"Hsig = %.2lf (m),  freq0 = %.2lf (Hz),  fwidth =%.2lf (Hz), Lwave = %.1lf (m) theta0 = %.3lf (deg)\n",
                  R2->Hsig, R2->freq0, R2->fwidth, R2->L, R2->theta0*180/pi);
  fprintf(stderr,"h0 = %.3lf (m), kh_0 = %.2lf,   a/h = %.2lf,  LY/ly=%d , k_y = %.5lf (rad/m)\n",R2->h0, R2->kh0, R2->aih0, R2->n_wavenumber, R2->ky);
  fprintf(stderr,"i_src = %d, i_width =%d,  x_src = %.1lf (m), x_width = %.1lf (m), delta=%.1f\n",ixcenter, num_gx, ixcenter*dx + 0.5*dx, num_gx*dx,delta);
  fprintf(stderr,"-------------------------------------------------------------------------\n");
}


void eta_source_report()
{

  if (eta_source_flag == OFF)
    return;

  switch(wave_type) {
  case MONOCHROMATIC: 
    eta_source_report_monochromatic(MW);
    break;
  case RANDOM1:
    eta_source_report_random1(RW1);
    break;
  case RANDOM2:
    eta_source_report_random2(RW2);
    break;
  }

}


double get_average_wavemaker_depth(const field2D *DD, int ixcenter)
{
  int j;
  double havg = 0.0;
  int M = DD->M;

  for (j=0;j<M;j++) {
    havg += DR2(DD,ixcenter,j);
  }
  havg = havg/M;
  return havg;

}


/* Wave math functions */


double k_tanh(double x)
{
  double t;

  t = sinh(x)/cosh(x);
  return t;

}

double k_sech(double x)
{
  double t;

  t = 1.0/cosh(x);
  return t;

}


double get_wavenum(double omega, double depth)
{

  double f, dfdk;
  double k, kh;


  k = omega/ sqrt(gravity*depth);
  kh = k*depth;
  f = gravity * k * k_tanh(kh) - SQR(omega);

  while (fabs(f)>1e-10) {
    dfdk = gravity*kh * SQR( k_sech(kh) ) + gravity*k_tanh(kh);
    k = k - f/dfdk;
    kh = k * depth;
    f = gravity*k* k_tanh(kh) - SQR( omega );
  }

  return k;

}

double get_cg(double k, double depth)
{

  double denom;
  double numerator;
  double cg;
  double kh = depth*k;

  denom = sqrt(gravity*k*k_tanh(kh));
  numerator = 0.5*(gravity*k_tanh(kh)+gravity*kh*SQR(k_sech(kh)));

  cg = numerator/denom;

  return cg;

}


int is_odd(int n)
{

  if ( n & 1 )
    return n;
  else
    return n+1;

}

int k_round_odd(double a)
{
  double frac,ip;
  int ii;

  frac = modf(a,&ip);
  
  ii = (int ) ip;

  if (frac>=0.5) 
    ii++;

  return is_odd(ii);  /* make sure it is odd */

}


double init_random_phase_generator()
{
  int seed;
  seed = (int ) time(NULL);
  srand(seed);

}

double random_phase()       /* returns a random phase between 0...2*pi */
{

  double r;
  
  r = (double ) rand();
  return  2.0*pi * ( r/RAND_MAX);

}


/* This function returns a new thetea to ensure that the specified wave angle 
   actually fits the model domain!                           */

double calculate_fit_wave_angle(double k_tot, double theta, int ny, double *ky_new, int *n_wavenumber)
{
  double LY;            /* the alongshore domain width  */
  double kY;
  double ky;
  double theta_new;

  LY = ny*dy;            
  kY = 2*pi/LY;             /* the alongshore domain wavenumber */
  ky = sin(theta) * k_tot;  /* the wave alongshore wavenumber: so ky/kY must be an integer */

  /* What is closest integer of ky/KY ? */

  *n_wavenumber = (int ) round( ky/kY );

  *ky_new = *n_wavenumber * kY;

  theta_new = asin( *ky_new/k_tot );

  return theta_new;

}


void deallocate_monochromatic_wave_type()
{

  g_free(MW->g_x);
  g_free(MW->s_ytA);
  g_free(MW->s_ytB);

}


void deallocate_random1_wave_type()
{

  /* finish this */

  g_free(RW1->omega);
  g_free(RW1->phase);
  g_free(RW1->amplitude);
  g_free(RW1->theta);
  g_free(RW1->g_x);
  g_free(RW1->rsinwt);
  g_free(RW1->rcoswt);
  g_free(RW1->ktot);
  g_free(RW1->ky);

  deallocate_field2D(RW1->s_ytA);
  deallocate_field2D(RW1->s_ytB);

}


void deallocate_random2_wave_type()
{

  /* finish this */

  fprintf(stderr,"** Warning: not finished writing random2 wave type deallocation function\n");

  g_free(RW2->omega);
  g_free(RW2->amplitude);
  g_free(RW2->mtheta);
  g_free(RW2->g_x);
  g_free(RW2->rsinwt);
  g_free(RW2->rcoswt);
  g_free(RW2->ktot);
  g_free(RW2->ky);

  deallocate_field2D(RW1->s_ytA);
  deallocate_field2D(RW1->s_ytB);

}

void  deallocate_eta_source()
{

  if (eta_source_flag == ON) {
    switch (wave_type) {
    case MONOCHROMATIC:  
      deallocate_monochromatic_wave_type();
      break;
    case RANDOM1: 
      deallocate_random1_wave_type();
      break;
    case RANDOM2: 
      deallocate_random2_wave_type();
      break;
    }
  }

}




/* void  rhs_calc_eta_source_f_random2_add(field2D *ETAT) */
/* { */


/*   double ggx; */
/*   double rfsrc; */
/*   double ky, y, coskyy, sinkyy; */
/*   double phase; */
/*   double cforc, sforc; */
/*   double samp, camp; */

/*   double *kyd; */
/*   //  double dwt, dw, wt; */

/*   int M = ETAT->M; */
/*   int i,j,k,l; */

/*   double *etatd = REF2(ETAT, ixstart); */

/*   double *romega = RW2->omega; */
/*   double *rcampd =  REF2(RW2->camplitude, 0); */
/*   double *rsampd =  REF2(RW2->samplitude, 0); */
/*   double *rphased = REF2(RW2->phase, 0); */
/*   double *rsinwt = RW2->rsinwt; */
/*   double *rcoswt = RW2->rcoswt; */
/*   double *g_x = RW2->g_x; */



/*   for (j=0;j<M;j++) { */
/*     rfsrc = 0.0; */
/*     y = j*dy; */
/*     kyd = REF2(RW2->ky,0); */
/*     for (k=0; k< RW2->numfreq; k++) {   /\* loop over all frequencies *\/ */
/*       cforc = 0.0; */
/*       sforc = 0.0; */
/*       for (l=0; l< RW2->numky; l++) {   /\* loop over all alongshore wavenumbers *\/ */
/* 	ky = *kyd++; */
/*         phase = *rphased++; */
/*         camp = *rcampd++; */
/*         samp = *rsampd++; */
/*         sinkyy = sin(ky*y + phase); */
/*         coskyy = cos(ky*y + phase); */
/*         cforc += camp*coskyy;  */
/*         sforc += samp*sinkyy; */
/*       } */
/*       rfsrc += cforc * cos( romega[k]*local_time ) + sforc * sin( romega[k]*local_time ); */

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

/* } */