lateral_mixing.c

波浪数值模拟

  for (i=1;i<N-1;i++) {
    for (j=0;j<M;j++) {
      utmp0 = (*xd) * 2.0;
      utmp1 =  *xdup++ - utmp0 + *xddn++;
      *lxd++ += idx2*utmp1;
      xd++; 
    }
  }

}


/* ***** BIHARMONIC 0 FUNCTIONS ******* */

void biharmonic0_u_all(field2D *BU, const field2D *LU)
{

  int i,j,im,M,MN;
  double tmp0, tmp1, tmp2;
  double *bud;
  double *lud, *ludup, *luddn;

  M = BU->M;
  MN = M-1;

  bud = BU->data;

  lud = &(LU->data[M]);
  ludup = &(LU->data[2*M]);
  luddn = LU->data;

  for (i=1;i<nx-1;i++) {
    im = i-1;
    tmp0 = (*lud) * 2.0;  // here we do j=0
    tmp1 =  *ludup++ - tmp0 + *luddn++;
    tmp2 = *(lud+1) - tmp0 + *(lud+MN);
    *bud++ = nu_bi*(idx2*tmp1 + idy2*tmp2);
    lud++;  
    for (j=1;j<ny-1;j++) {
      tmp0 = (*lud) * 2.0;
      tmp1 =  *ludup++ - tmp0 + *luddn++;
      tmp2 = *(lud+1) - tmp0 + *(lud-1);
      *bud++ = nu_bi*(idx2*tmp1 + idy2*tmp2);
      lud++;  
    }
    tmp0 = (*lud) * 2.0;  // here we do j=ny-1
    tmp1 =  *ludup++ - tmp0 + *luddn++;
    tmp2 = *(lud-MN) - tmp0 + *(lud-1);
    *bud++ = nu_bi*(idx2*tmp1 + idy2*tmp2);
    lud++; 
  }
}



void biharmonic0_v_all(field2D *BV, const field2D *LV)
{

  int i,j,im,M,MN;
  double tmp0, tmp1, tmp2;
  double *lvd;
  double *vd, *vdup, *vddn;


  M = BV->M;
  MN = M-1;

  lvd = BV->data;

  vd = &(LV->data[M]);
  vdup = &(LV->data[2*M]);
  vddn = LV->data;

  for (i=1;i<nx;i++) {
    im = i-1;
    tmp0 = (*vd) * 2.0;  // here we do j=0
    tmp1 =  *vdup++ - tmp0 + *vddn++;
    tmp2 = *(vd+1) - tmp0 + *(vd+MN);
    *lvd++ = nu_bi*(idx2*tmp1 + idy2*tmp2);
    vd++;  
    for (j=1;j<ny-1;j++) {
      tmp0 = (*vd) * 2.0;
      tmp1 =  *vdup++ - tmp0 + *vddn++;
      tmp2 = *(vd+1) - tmp0 + *(vd-1);
      *lvd++ = nu_bi*(idx2*tmp1 + idy2*tmp2);
      vd++;  
    }
    tmp0 = (*vd) * 2.0;  // here we do j=ny-1
    tmp1 =  *vdup++ - tmp0 + *vddn++;
    tmp2 = *(vd-MN) - tmp0 + *(vd-1);
    *lvd++ = nu_bi*(idx2*tmp1 + idy2*tmp2);
    vd++;  
  }
}



/* ***** BIHARMONIC 1 FUNCTIONS ******* */

void biharmonic1_u_all(field2D *BU, field2D *LU)
{

  bi_laplacian0_all(BU,LU);

}


void biharmonic1_v_all(field2D *BV, field2D *LV)
{


  bi_laplacian0_all(BV,LV);

}



/* void mult_laplacian0_depth(field2D *LU, field2D *LV) */
/* { */

/*   int i,j,im,M,N; */
/*   double *lud, *lvd, *lvd0, *lvd1, *lvd2; */
/*   double hbar, him; */

/*   M = LU->M; */
/*   N = LU->N; */

/*   lud = &(LU->data[M]); */
/*   lvd = &(LV->data[M]); */

/*   for (i=1;i<N-1;i++) { */
/*     him = HH[i-1]; */
/*     hbar = HBAR[i-1]; */
/*     for (j=0;j<M;j++) { */
/*       *lud++ *= hbar; */
/*       *lvd++ *= him; */
/*     } */
/*   } */

/*   him = HH[N-2]; */
/*   for (j=0;j<M;j++) */
/*       *lvd++ *= him; */

/*   /\* now deal with the boundary conditions *\/ */

/*   lvd0 = LV->data; */
/*   lvd1 = lvd0 +M; */
/*   lvd2 = lvd - M; */

/*   for (j=0;j<M;j++) { */
/* #ifdef SLIP */
/*     *lvd0++ = *lvd1++; */
/*     *lvd++ = *lvd2++; */
/* #endif  /\* SLIP *\/ */
/* #ifdef NOSLIP     /\* V=0,  h*lap(V) = 0 *\/ */
/*     *lvd0++ = - *lvd1++; */
/*     *lvd++ = - *lvd2++; */
/* #endif /\* NOSLIP *\/ */
/*   } */

/* } */


void biu_dxdy_all_add(field2D *BU, field2D *LV)
{

  int i,im,j, M, MN, N;
  double *lvd, *lvdm, *lvdup, *lvdupm;
  double *bud;
  double tmp1, tmp0;

  M = BU->M;
  N = LV->N;
  MN = M-1;

  bud = BU->data;
  lvd = &(LV->data[M]);
  lvdup = &(LV->data[2*M]);

  lvdm = lvd-1;
  lvdupm = lvdup-1;

  for (i=1;i<N-2;i++) {
    tmp0 = (*lvdup++ - *(lvdupm+M));
    tmp1 = (*lvd++ - *(lvdm+M));
    *bud++ += (tmp0-tmp1)*idx*idy;
    lvdm++;
    lvdupm++;
    for (j=1;j<M;j++) {
      tmp0 = (*lvdup++ - *lvdupm++);
      tmp1 = (*lvd++ - *lvdm++);
      *bud++ += (tmp0-tmp1)*idx*idy;
    }
  }
}


void biv_dxdy_all_add(field2D *BV, field2D *LU)
{

  int i,im,j, M, MN, N;
  double *lud, *ludp, *luddn, *luddnp;
  double *bvd;
  double tmp1, tmp0;

  M = BV->M;
  N = LU->N;
  MN = M-1;

  bvd = BV->data;
  lud = &(LU->data[M]);
  luddn = LU->data;

  ludp = lud+1;
  luddnp = luddn+1;

  for (i=1;i<N;i++) {
    for (j=0;j<M-1;j++) {
      tmp0 = (*luddnp++ - *luddn++);
      tmp1 = (*ludp++ - *lud++);
      *bvd++ += (tmp1-tmp0)*idx*idy;
    }
    tmp0 = (*(luddnp-M) - *luddn++);
    tmp1 = (*(ludp-M) - *lud++);
    *bvd++ += (tmp1-tmp0)*idx*idy;
    ludp++;
    luddnp++;

  }
}



void biharmonic2_u_all(field2D *BU, field2D *LU, field2D *LV)
{


  bi_dxx_all(BU,LU);
  field2D_mult_const(BU,2.0);
  bi_dyy_all_add(BU,LU);

  biu_dxdy_all_add(BU,LV);

  field2D_mult_const(BU,nu_bi);
}

void biharmonic2_v_all(field2D *BV, field2D *LU, field2D *LV)
{

  bi_dyy_all(BV,LV);
  field2D_mult_const(BV,2.0);
  bi_dxx_all_add(BV,LV);

  biv_dxdy_all_add(BV,LU);

  field2D_mult_const(BV,nu_bi);


}


/* Returns   nu_bi * \del^4 ( U, V)  
   where  DU = (nx-2,m) and DV = (nx-1,m)    */

void biharmonic0_mixing(field2D *DU, field2D *DV, const field2D *U, const field2D *V)
{

    laplacian0_all(LUU, LVV, U, V);
    biharmonic0_u_all(DU, LUU);
    biharmonic0_v_all(DV, LVV);

}

/* Returns   nu_newt * \laplacian( U, V)  
   where  DU = (nx-2,m) and DV = (nx-1,m)    */

void newtonian0_mixing(field2D *DU, field2D *DV, const field2D *U, const field2D *V)
{

  field2D_laplacian(DU, U,nu_newt*idx,nu_newt*idy);   /*  DU = nu_newt * \laplacian (U) */

  field2D_laplacian(DV, V,nu_newt*idx,nu_newt*idy);   /*  DV = nu_newt * \laplacian (V) */


}


/* Returns   nu_newt * \laplacian( U, V)  
   where  DU = (nx-2,m) and DV = (nx-1,m)    */

void newtonian1_mixing(field2D *DU, field2D *DV, const field2D *U, const field2D *V, const field2D_depth_var *DD)
{

  int NU = U->N -1;
  int M = U->M;

  static field2D *Ux=NULL, *Uy=NULL;
  static field2D *Vx=NULL, *Vy=NULL;

  if ( NOT_ALLOCATED( Ux ) ) {
    Ux = allocate_field2D_grid( NU, M, UGRID);
    Uy = allocate_field2D_grid( NU, M, UGRID);
    Vx = allocate_field2D_grid( NU, M, UGRID);
    Vy = allocate_field2D_grid( NU, M, UGRID);
  }

  funwaveC_error_message("** lateral_mixing.c: newtonian1_mixing() is broken!");

  field2D_diffx(Ux, U,idx);           /*  U_x = \delta_x (U)  */

  field2D_mult_const(DU, nu_newt);    /*  DU = nu_newt * \laplacian (U) */

  field2D_laplacian(DV, V,idx,idy);   /*  DV = \laplacian (V) */
  field2D_mult_const(DV, nu_newt);

}

void rhs_uvmom_lateral_mixing(field2D *DU, field2D *DV, const field2D *U, const field2D *V, const field2D_depth_var *DD)
{


  switch(lateral_mixing_type) {
  case NEWT0:  
    newtonian0_mixing(DU, DV, U, V);
    break;
  case NEWT1:  
    funwaveC_error_message("** lateral_mixing.c: newtonian1_mixing() is broken!");
 /*    newtonian1_u_all(DU); */
/*     newtonian1_v_all(DV); */
    break;
  case NEWT2:  
    funwaveC_error_message("** lateral_mixing.c: newtonian2_mixing() is broken!");
 /*    newtonian2_u_all(DU); */
/*     newtonian2_v_all(DV); */
    break;
  case BI0:  
    biharmonic0_mixing(DU, DV, U, V);
    break;
  case BI1:  
    funwaveC_error_message("** lateral_mixing.c: biharmonic1_mixing() is broken!");
 /*    laplacian0_all(LUU, LVV, U, V); */
/*     mult_laplacian0_depth(LUU,LVV); */
/*     biharmonic1_u_all(DU,LUU); */
/*     biharmonic1_v_all(DV,LVV); */
    break;
  case BI2:  
    funwaveC_error_message("** lateral_mixing.c: biharmonic2_mixing() is broken!");
 /*    laplacian0_all(LUU,LVV, U, V); */
/*     mult_laplacian0_depth(LUU,LVV); */
/*     biharmonic2_u_all(DU,LUU,LVV); */
/*     biharmonic2_v_all(DV,LUU,LVV); */
    break;
  case NONE:
    return;
  }
}



void init_lateral_mixing(int nx, int ny, mix_t DT, double nu1, double nu2)
{


  nu_newt = nu1;  nu_bi = nu2;
  lateral_mixing_type = DT;

  LUU = allocate_field2D(nx,ny);
  LVV = allocate_field2D(nx+1,ny);

  field2D_set_to_value(LUU,0.0);
  field2D_set_to_value(LVV,0.0);
}


void deallocate_lateral_mixing()
{

  deallocate_field2D(LUU);
  deallocate_field2D(LVV);

}


/* U is the max current, L is the domain length, and dt is the time step */

void lateral_mixing_stability_report( double Umax, double L, double dt )
{

  double Sbix, Rbix; 
  double Sbiy, Rbiy, RbiLy;
  double RbiL;

  double Sntx, Snty;
  double Rntx, Rnty, RntL;

    // Biharmonic Friction Sbi & Grid Reynolds #

    fprintf(stderr,"---------Lateral Mixing Stability Report ----------------------\n");

    switch(lateral_mixing_type) {
    case BI0:
    case BI1:
    case BI2:
      Sbix = fabs(nu_bi)*dt*idx*idx*idx*idx;     // < 1.2e-2 (AB2)  <0.8e-2 (AB3)
      Sbiy = fabs(nu_bi)*dt*idy*idy*idy*idy;     // < 1.2e-2 (AB2)  <0.8e-2 (AB3)
      fprintf(stderr,"Using Biharmonic Friction:  Sbi_x # = %g , Sbi_y # = %g  : ",Sbix, Sbiy);

      if ( (Sbix > 0.008) || (Sbiy > 0.008) )
          fprintf(stderr,"Warning: model may be Biharmonic Friction unstable. Sbi_[x,y] <  8e-3 (AB3)\n");
      else
          fprintf(stderr,"Biharmonic Friction OK\n");

      /* now grid reynolds numbers and then domain reynolds #'s */

      Rbix = fabs( Umax/(nu_bi*idx*idx*idx) );  // < 13 (AB2) - what about AB3?
      Rbiy = fabs( Umax/(nu_bi*idy*idy*idy) );  // < 13 (AB2)

      fprintf(stderr,"Grid Reynolds Number: R_bi-x # = %g , R_bi-y # = %g  : ",Rbix, Rbiy);
      if ( (Rbix>13) || (Rbiy>13) )
          fprintf(stderr,"Warning: Biharmonic Grid Re may be too small - not enough dissipation\n");
      else
          fprintf(stderr,"Biharmonic Grid Re OK\n");

      RbiL = fabs( Umax*L*L*L/nu_bi );
      fprintf(stderr,"Domain Scale Biharmonic Re # = %g\n",RbiL);
      break;

    case NEWT0:
    case NEWT1:
    case NEWT2:
      Sntx = nu_newt*dt*idx*idx;    // < ??????
      Snty = nu_newt*dt*idy*idy;    // < ??????
      fprintf(stderr,"Newtonian Mixing Stability Number: Snt_x # = %g,  Snt_y # =%g :   ",Sntx, Snty);

      if ( (Sntx>0.13) || (Snty>0.13) )
          fprintf(stderr,"Warning: model may be Newtonain Friction unstable. Snt < 0.13 (AB2) or ??? (AB3)\n");
      else
          fprintf(stderr,"Newtonian Friction OK\n");

      Rntx = Umax/(nu_newt*idx);    // < ??????
      Rnty = Umax/(nu_newt*idy);    // < ??????
      RntL = Umax*L/nu_newt;

      fprintf(stderr,"Grid Reynolds Number:  Rnt_x # = %g, Rnt_y # = %g \n",Rntx, Rnty);
      fprintf(stderr,"Domain Scale Newtonain Re # = %g  \n",RntL);
      break;