bousinessq_dynamics.c

波浪数值模拟

void b_calc_uvmom_F1t_G1t(field2D  *F1t, field2D *G1t, const field2D *HBX, const field2D *HBY, const field2D *Ut, const field2D *Vt, const field2D *HUt, const field2D *HVt)
{

   b_calc_umom_F1(F1t, HBX, Vt, HVt);
   b_calc_vmom_G1(G1t, HBY, Ut, HUt);
}






void b_calc_uvmom_FG2(field2D *FG2, const field2D *ETA, const field2D *ZA, const field2D *U, const field2D *V, const field2D *HU, const field2D *HV)
{
  const int NU = U->N;
  const int NETA = ETA->N;
  const int M = U->M;

  static field2D *UXPVY_ETA=NULL, *HUXPVY_ETA=NULL;     /* workspace on an ETA grid point */
  static field2D *WORK_U=NULL, *WORK_V=NULL, *WORK_E=NULL, *WORK_E1=NULL, *WORK_E2=NULL, *WORK_E3=NULL, *WORK_E4=NULL;

  const double *zad, *etad;
  double *worke2d, *worke3d;
  double *uxpvyd, *huxpvyd, *fg2d;
  double za, eta;
  double tmp;
  int i,j;

  if ( NOT_ALLOCATED( UXPVY_ETA ) ) {
    UXPVY_ETA = allocate_field2D_grid( NU-1, M, ETAGRID);
    HUXPVY_ETA = allocate_field2D_grid( NU-1, M, ETAGRID);
    WORK_U = allocate_field2D_grid( NU, M, UGRID);
    WORK_V = allocate_field2D_grid( NETA, M, VGRID);
    WORK_E = allocate_field2D_grid( NETA, M, ETAGRID);
    WORK_E1 = allocate_field2D_grid( NETA, M, ETAGRID);
    WORK_E2 = allocate_field2D_grid( NETA, M, ETAGRID);
    WORK_E3 = allocate_field2D_grid( NETA, M, ETAGRID);
    WORK_E4 = allocate_field2D_grid( NETA, M, ETAGRID);
  }


  /* Regular code: WORK_E2 = 0.5*( z_a^2 - eta^2 ) 
                   WORK_E3 = ( z_a - eta )    */

  worke2d = WORK_E2->data;
  worke3d = WORK_E3->data;
  zad = &(ZA->data[M]);     /* ZA is on a full H grid - so start 1 in */
  etad = ETA->data;

  for (i=0;i<NETA;i++) {
    for (j=0;j<M;j++) {
      za = *zad++;
      eta = *etad++;
      *worke2d++ = 0.5 * ( za*za - eta*eta);
      *worke3d++ =  ( za - eta);
    }
  }


  /* This does the first line of FG2:  in brackets of equation 15 */

  b_calc_ux_plus_vy(UXPVY_ETA, U, V, idx, idy);  /* (u_x + v_y) on ETA point */
  field2D_diffx(WORK_U, UXPVY_ETA, idx);  /*  (u_x + v_y)_x on U point */
  field2D_self_multiply(WORK_U, U);       /*  WORK_U = u (u_x + v_y)_x */
  field2D_barx(WORK_E, WORK_U);            /*  WORK_E = \bar^x{u (u_x + v_y)_x}  on an eta point  */
  field2D_diffy_eta_to_v(WORK_V, UXPVY_ETA, idy);  /*  WORK_V = { (u_x + v_y)_y}  on an V point  */
  field2D_self_multiply(WORK_V, V);               /* WORK_V = { v(u_x + v_y)_y}  on an V point  */
  field2D_bary_v_to_eta(WORK_E1, WORK_V);         /* WORK_E1 = \bar^y { v(u_x + v_y)_y}  but on ETA  points  */
  field2D_self_add(WORK_E, WORK_E1);              /* WORK_E =  \bar^x{u (u_x + v_y)_x} + \bar^y { v(u_x + v_y)_y} on ETA POINTS */
  field2D_multiply(FG2, WORK_E, WORK_E2);         /* FG = (1/2)(z^2-eta^2) * \bar^x{u (u_x + v_y)_x} + \bar^y { v(u_x + v_y)_y} on ETA POINTS */


  /* This does the 2nd line  line of FG2:  in brackets of equation 15 */

  b_calc_ux_plus_vy(HUXPVY_ETA, HU, HV, idx, idy);   /* (hu_x + hv_y) on ETA point */
  field2D_diffx(WORK_U, HUXPVY_ETA, idx);  /*  (hu_x + hv_y)_x on U point */
  field2D_self_multiply(WORK_U, U);       /*  WORK_U = u (hu_x + hv_y)_x */
  field2D_barx(WORK_E, WORK_U);            /*  WORK_E = \bar^x{u (hu_x + hv_y)_x}  on an eta point  */
  field2D_diffy_eta_to_v(WORK_V, HUXPVY_ETA, idy);  /*  WORK_V = { (hu_x + hv_y)_y}  on an V point  */
  field2D_self_multiply(WORK_V, V);                /* WORK_V = { v(hu_x + hv_y)_y}  on an V point  */
  field2D_bary_v_to_eta(WORK_E1, WORK_V);          /* WORK_E1 = \bar^y { v(hu_x + hv_y)_y}  but on ETA  points  */
  field2D_self_add(WORK_E, WORK_E1);               /* WORK_E =  \bar^x{u (hu_x + hv_y)_x} + \bar^y { v(hu_x + hv_y)_y} on ETA POINTS */
  field2D_multiply_add(FG2, WORK_E, WORK_E3);      /* FG = FG + (z-eta) * ( \bar^x{u (hu_x + hv_y)_x} + \bar^y { v(hu_x + hv_y)_y} on ETA POINTS */

  /* This does the 3rd line of FG2: in brackets of equation 15 of funwave manual */

  huxpvyd = HUXPVY_ETA->data;
  uxpvyd  = UXPVY_ETA->data;
  etad =    ETA->data;
  fg2d =     FG2->data;

  for (i=0;i<NETA;i++) {
    for (j=0;j<M;j++) {
      tmp =  *huxpvyd++ + (*etad++) * (*uxpvyd++);
      *fg2d++ += 0.5 * tmp*tmp;
    }
  }

}


/* This routinue calculates F2 and G2 on the RHS of the U & V momentum equations.
   These terms are the gradient of a scalar defined on ETA points.   
    1) Calculate the scalar FG2 on the ETA points
         - this scalar field FG2 is passed to this function as a WORK array
    2) calculate the \delta_x and \delta_y  of the scalar to get F2 and G2  */

/* void b_calc_uvmom_F2_G2(field2D *F2, field2D *G2, const field2D *ETA,  const field2D *ZA, const field2D *U, const field2D *V,  */
/*                         const field2D *HU, const field2D *HV) */
/* { */

/*   int M = F2->M; */
/*   int N = ETA->N; */

/*   static field2D *FG2=NULL; */

/*   if ( NOT_ALLOCATED( FG2 ) ) { */
/*     FG2 = allocate_field2D_grid( N, M, ETAGRID);  /\* FG2  need to be on an ETA grid for gradient to be on a U, V grid *\/ */
/*   } */

/*   b_calc_uvmom_FG2(FG2, ETA, ZA, U, V, HU, HV ); */
/*   field2D_diffx(F2, FG2, -idx*gamma);                  /\* F2 = - \delta_x( FG2 ) - note the minus sign! Also gamma *\/ */
/*   field2D_diffy_eta_to_v(G2, FG2, -idy*gamma);         /\* G2 = - \delta_y( FG2 ) - note the minus sign! Also gamma *\/ */

/* } */




/* Calculates field2D *CCT:  the precursor of F^t and G^t in the momentum equations.  The gradient of this
   term goes into the momentum equation:  ie  F^t = \delta_x (CCT)                        
   CCT = (1/2)*\eta^2 * [ (dU/dt)_x + (dV/dt)_y] + \eta*[ (h dU/dt)_x + (h dV/dt)_y ]

   numerically:  CCT is on a eta points.

   CCT = (1/2)*\eta^2 * [ \delta_x(dU/dt) + \delta_y(dV/dt)] + \eta*[ \delta_x (\barx{h} dU/dt) + \delta_y (\bary{h} dV/dt) ]

   ** dU/dt must be on same grid as U:   dV/dt needs only be on N-2 grid relative to V.   
   ** HBX - needs to be on full U grid.   HBY on an inside grid (N-2)

                  */

void b_calc_uvmom_CCT(field2D  *CCT, const field2D *ETA, const field2D *UT, const field2D *VT, const field2D *HUT, const field2D *HVT )
{
  int M = CCT->M;
  int MN = M-1;
  int NETA = ETA->N;

  const double *etad = &(ETA->data[0]);
  const double *utd = &(UT->data[0]);
  const double *vtd = &(VT->data[M]);
  const double *hutd = &(HUT->data[0]);
  const double *hvtd = &(HVT->data[M]);

  double *cctd = &(CCT->data[0]);
  double tmp1, tmp2, eta, eta2;
  int i,j;

  for (i=0;i<NETA;i++) {
    eta = *etad++;
    eta2 = eta*eta;
    tmp1 = 0.5*eta2*( idx*( *(utd+M)-*utd) + idy*( *vtd - *(vtd+MN) ) );
    tmp2 = eta* ( idx*( *(hutd+M) - *hutd)  + idy*( (*hvtd) - (*(hvtd+MN))));
    *cctd++ = tmp1 + tmp2; 
    utd++;
    vtd++;
    hutd++;
    hvtd++;
    //    hbxd++;
    //    hbyd++;

    for (j=1;j<M;j++) {
      eta = *etad++;
      eta2 = eta*eta;
      tmp1 = 0.5*eta2*( idx*( *(utd+M)-*utd) + idy*( *vtd - *(vtd-1) ) );
      tmp2 = eta* ( idx*( *(hutd+M) - *hutd)  + idy*( (*hvtd) - (*(hvtd-1))));
      *cctd++ = tmp1 + tmp2; 
      utd++;
      vtd++;
    hutd++;
    hvtd++;
    //    hbxd++;
    //    hbyd++;
    }
  }

}


void b_calc_uvmom_CCT_add(field2D  *CCT, const field2D *ETA, const field2D *UT, const field2D *VT, const field2D *HUT, const field2D *HVT )
{
  int M = CCT->M;
  int MN = M-1;
  int NETA = ETA->N;

  const double *etad = &(ETA->data[0]);
  const double *utd = &(UT->data[0]);
  const double *vtd = &(VT->data[M]);
  const double *hutd = &(HUT->data[0]);
  const double *hvtd = &(HVT->data[M]);

  double *cctd = &(CCT->data[0]);
  double tmp1, tmp2, eta, eta2;
  int i,j;

  for (i=0;i<NETA;i++) {
    eta = *etad++;
    eta2 = eta*eta;
    tmp1 = 0.5*eta2*( idx*( *(utd+M)-*utd) + idy*( *vtd - *(vtd+MN) ) );
    tmp2 = eta* ( idx*( *(hutd+M) - *hutd)  + idy*( (*hvtd) - (*(hvtd+MN))));
    *cctd++ += tmp1 + tmp2; 
    utd++;
    vtd++;
    hutd++;
    hvtd++;
    //    hbxd++;
    //    hbyd++;

    for (j=1;j<M;j++) {
      eta = *etad++;
      eta2 = eta*eta;
      tmp1 = 0.5*eta2*( idx*( *(utd+M)-*utd) + idy*( *vtd - *(vtd-1) ) );
      tmp2 = eta* ( idx*( *(hutd+M) - *hutd)  + idy*( (*hvtd) - (*(hvtd-1))));
      *cctd++ += tmp1 + tmp2; 
      utd++;
      vtd++;
    hutd++;
    hvtd++;
    //    hbxd++;
    //    hbyd++;
    }
  }

}


/* Calculates the F^t and G^t terms in eq 17,18 of funwave manual.  These terms are a grad of a scalar,
   denoted CCT.    So first calculate the scalar and then do the gradient.   Relatively simple         
   Note that CCT is a dummy field2D*   and needs to be passed a memory allocation when it is called     */


/* void b_calc_uvmom_Ft_Gt(field2D *Ft, field2D *Gt, const field2D *ETA, const field2D *UT, const field2D *VT, const field2D *HUT, const field2D *HVT) */
/* { */

/*   int M = Ft->M; */
/*   int N = Ft->N + 1; */

/*   static field2D *CCT=NULL; */

/*   if ( NOT_ALLOCATED( CCT ) ) { */
/*     CCT = allocate_field2D_grid( N, M, ETAGRID);  /\* CCT need to be on an ETA grid for gradient to be on a U, V grid *\/ */
/*   } */

/*   b_calc_uvmom_CCT(CCT, ETA, UT, VT, HUT, HVT ); */
/*   field2D_diffx(Ft, CCT, idx);                   /\* Ft = \delta_x (CCT) - note no minus sign here - w/ gamma *\/ */
/*   field2D_diffy_eta_to_v(Gt, CCT, idy);          /\* Gt = \delta_y (CCT) - note no minus sign here - w/ gamma *\/ */

/* } */


void b_calc_uvmom_FG_WK(field2D *FWK, field2D *GWK, const field2D *ETA, const field2D *ZA, 
                        const field2D *U, const field2D *V, const field2D *HU, const field2D *HV,
			const field2D *UT, const field2D *VT, const field2D *HUT, const field2D *HVT)
{

  int M = FWK->M;
  int N = ETA->N;

  static field2D *FGWK=NULL;

  if ( NOT_ALLOCATED( FGWK ) ) {
    FGWK = allocate_field2D_grid( N, M, ETAGRID);  /* FG2  need to be on an ETA grid for gradient to be on a U, V grid */
  }

  b_calc_uvmom_FG2(FGWK, ETA, ZA, U, V, HU, HV );
  field2D_mult_const(FGWK, -1.0);
  b_calc_uvmom_CCT_add(FGWK, ETA, UT, VT, HUT, HVT );

  /* FWK = - \delta_x( FGWK )    GWK = - \delta_x( GWK )  */

  field2D_gradient_eta_to_uv(FWK, GWK, FGWK, idx, idy); 


}



/* Bousinessq model calculate the shallow water pressure gradient ( g \grad{ \eta } ) */
/* Works for both U and V pressure gradients */

void b_calc_pressure_gradient(field2D *PGU, field2D *PGV, const field2D *ETA)
{

   int i,j;

   static const double grav = 9.81;
   const int M = PGU->M;
   //   const int MN = M-1;
   const int NETA = ETA->N;

   const double *etad1 = &(ETA->data[0]);
   const double *etad2 = &(ETA->data[M]);

   double etap, etam;

   double *pgud = &(PGU->data[0]);
   double *pgvd = &(PGV->data[0]);


   /* First calculate grav * d(eta)/dx */

   for (i=1;i<NETA;i++) {
     for (j=0;j<M;j++) 
       *pgud++ = -idx*grav*(*etad2++ - *etad1++);
   }

   /* Now calculate grav * d(eta)/dy */

   etad1 = &(ETA->data[0]);
 
   for (i=0;i<NETA;i++) {
     etam = *etad1++;
     for (j=0;j<M-1;j++) {
       etap = *etad1++;
       *pgvd++  = -idy*grav*( etap - etam);
       etam = etap;
     }
     *pgvd++ = -idy*grav*( *(etad1 - M) -  etam );
   }

}







/* This subroutinue calculates the first part of d(eta)/dt = E + ....   
   where E = EP1 + ....       EP1 = (1/kappa) * [ { (h+\eta)*u }_x + { (h+\eta)*v }_y ]
   not that in the funwave manual  (h+\eta) is replaced by \Lambda for the slot technique stuff  

   NUMERICS:    (1/kappa) [ \delta_x (h^{bx}* u) + \delta_y (h^{by} * v) )
   Assumes that h^{bx} and h^{by}  are already calculated   

   Note that the depth input can be either h or (h+\eta).   It needs to be h for linear solutions

   */

void b_calc_eta_div_hu(field2D *ET, const field2D *U, const field2D *V, const field2D *DEPTHBX, const field2D *DEPTHBY)
{

  static double ikappa = 1.0;
  int M = U->M;
  int MN = M-1;
  int i,j;
  int NET = ET->N;

  const double *ud = &(U->data[0]);
  const double *vd = &(V->data[M]);
  const double *depthbxd = REF2(DEPTHBX, 0);
  const double *depthbyd = REF2(DEPTHBY, 0);
  
  double *ep1d = &(ET->data[0]);
  double hux, hvy;
  double depthx_up, depthx_dn;
  double depthy_pl, depthy_mn;

  for (i=0;i<NET;i++) {

    /* here is j = 0 */

    depthx_up = *(depthbxd+M);
    depthx_dn = *depthbxd;
    depthy_mn = *(depthbyd+MN);
    depthy_pl = *depthbyd++;


    hux =  idx * ( depthx_up * (*(ud+M) ) -  depthx_dn  * (*ud) );
    hvy =  idy * ( depthy_pl * (*vd ) - depthy_mn * (*(vd+MN)) );
    *ep1d++ = -ikappa * (hux + hvy);

    ud++;                            /* advance the pointers */
    vd++;
    depthbxd++;