bousinessq_dynamics.c

波浪数值模拟

    for (j=1;j<M;j++) {

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


      hux =  idx * (  depthx_up * (*(ud+M) ) -  depthx_dn * (*ud) );
      hvy =  idy * (  depthy_pl * (*vd )     -  depthy_mn * (*(vd-1)) );
      *ep1d++ = -ikappa * (hux + hvy);
      ud++;                            /* advance the pointers */
      vd++;
      depthbxd++;
    }
  }

}


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

/* void b_calc_eta_EP1(field2D *EP1, 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 NEP1 = EP1->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 = &(EP1->data[0]); */
/*   double hux, hvy; */
/*   double depthx_up, depthx_dn; */
/*   double depthy_pl, depthy_mn; */

/*   for (i=0;i<NEP1;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++; */



/*     for (j=1;j<M;j++) { */

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


/*       hux =  idx * (  depthx_up * (*(ud+M) ) -  depthx_dn * (*ud) ); */
/*       hvy =  idy * (  depthy_pl * (*vd )     -  depthy_mn * (*(vd-1)) ); */
/*       *ep1d++ = -ikappa * (hux + hvy); */
/*       ud++;                            /\* advance the pointers *\/ */
/*       vd++; */
/*       depthbxd++; */
/*     } */
/*   } */

/* } */


/* This does the 2nd line of E (eq 9 in funwave manual) 
   -{ a1 h^3 (u_xx + v_xy) + a2 h^2 ........ }_x   
   discretized on a cgrid                            */

void b_calc_eta_EP2(field2D *EP2, const field2D *HBX, const field2D *UXX, const field2D *VXY, const field2D *HUXX, const field2D *HVXY)
{

  const int NEP2 = EP2->N;
  const int M = EP2->M;

  int i,j;
  
  const double *uxxd = &(UXX->data[0]);
  const double *vxyd = &(VXY->data[0]);
  const double *huxxd = &(HUXX->data[0]);
  const double *hvxyd = &(HVXY->data[0]);
  const double *hbxd = &(HBX->data[0]);
  
  double *ep2d = &(EP2->data[0]);
  double hup,hdn,hup2,hup3,hdn2,hdn3;
  double tmp1, tmp0;

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

    for (j=0;j<M;j++) {

      hup = *(hbxd+M);
      hdn = *hbxd;
      hup2 = hup*hup;
      hdn2 = hdn*hdn;
      hup3 = hup*hup2;  /* h^3 up */
      hdn3 = hdn*hdn2;  /* h^3 down */

      tmp1 =  a1*hup3 * (*(uxxd+M) + *(vxyd+M) ) + a2* hup2 * (*(huxxd+M) + *(hvxyd+M) );
      tmp0 =  a1*hdn3 * (*uxxd + *vxyd ) + a2 * hdn2 * (*huxxd + *hvxyd );

      *ep2d++ += -idx*(tmp1-tmp0);

      uxxd++;                            /* advance the pointers */
      vxyd++;
      huxxd++;
      hvxyd++;
      hbxd++;
    }
  }


}

/* This does the 3nd line of E (eq 9 in funwave manual) 
   -{ a1 h^3 (u_xy + v_yy) + a2 h^2 ........ } _y 
   discretized on a cgrid                            */

void b_calc_eta_EP3(field2D *EP3, const field2D *HBY, const field2D *VYY, const field2D *UXY, const field2D *HVYY, const field2D *HUXY)
{

  int M = EP3->M;
  int MN = M-1;
  int NETA = EP3->N;

  int i,j;

  /* VYY, UXY, HVYY, and HUXY are all on V  grids so start data pointers at &( ..[0]);
     HBY is on a V grid to start at  &(...[M]);  */

  const double *vyyd = &(VYY->data[0]);
  const double *uxyd = &(UXY->data[0]);
  const double *hvyyd = &(HVYY->data[0]);
  const double *huxyd = &(HUXY->data[0]);
  const double *hbyd = &(HBY->data[M]);
  
  double *ep3d = &(EP3->data[0]);

  double hp,hm,hp2,hp3,hm2,hm3;
  double tmp1, tmp0;

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

    /* j = 0 */

    hm = *(hbyd+MN);
    hp = *hbyd++;

    hp2 = hp*hp;
    hm2 = hm*hm;
    hp3 = hp*hp2;  /* plus h^3 */
    hm3 = hm*hm2;  /* minus h^3 */

    tmp1 =  a1*hp3 * (*(vyyd) + *(uxyd) ) + a2* hp2 * (*(hvyyd) + *(huxyd) );
    tmp0 =  a1*hm3 * (*(vyyd+MN) + *(uxyd+MN) ) + a2 * hm2 * (*(hvyyd+MN) + *(huxyd+MN) );

    *ep3d++ += -idy*(tmp1-tmp0);

    vyyd++;                            /* advance the pointers */
    uxyd++;
    hvyyd++;
    huxyd++;

    tmp0 = tmp1;
    hm = hp;

    for (j=1;j<M;j++) {

      hp = *hbyd++;
      hp2 = hp*hp;
      hm2 = hm*hm;
      hp3 = hp*hp2;
      hm3 = hm*hm2;

      tmp1 =  a1*hp3 * (*(vyyd) + *(uxyd) ) + a2* hp2 * (*(hvyyd) + *(huxyd) );

      *ep3d++ += -idy*(tmp1-tmp0);

      tmp0 = tmp1;
      hm = hp;

      vyyd++;                            /* advance the pointers */
      uxyd++;
      hvyyd++;
      huxyd++;

    }
  }

}


/* Calculates the E term on the rhs of the eta equation.  It is broken down into 3 components EP1, EP2, and EP3 */

void b_calc_eta_E(field2D *E, const field2D *U, const field2D *V, const field2D *DEPTHBX, const field2D *DEPTHBY,
                  const field2D *HBX, const field2D *HBY, const field2D *UXX, const field2D *VXY, const field2D *HUXX, const field2D *HVXY,
                  const field2D *UXY, const field2D *VYY, const field2D *HUXY, const field2D *HVYY)
{


  b_calc_eta_div_hu(E, U, V, DEPTHBX, DEPTHBY);  /* E = -(1/kappa) div (h [u,v])  */
  b_calc_eta_EP2(E, HBX, UXX, VXY, HUXX, HVXY);     /* E += - ( )_x  */
  b_calc_eta_EP3(E, HBY, VYY, UXY, HVYY, HUXY);     /* E += - ( )_y  */


}


void b_calc_eta_E2PX(field2D *E2, const field2D *ETABX, const field2D *HBX, const field2D *UXX, const field2D *VXY, 
                     const field2D *HUXX, const field2D *HVXY)
{

  int NU = UXX->N;
  int M = E2->M;
  int i,j;

  static field2D *WORK_U=NULL;

  const double *hbxd = HBX->data;        /* Assumes on a U grid */
  const double *etabxd = ETABX->data;    /* Assumes on a U grid */
  const double *vxyd = VXY->data;        /* Assumes on a U grid */
  const double *uxxd = UXX->data;        /* Assumes on a U grid */
  const double *hvxyd = HVXY->data;      /* Assumes on a U grid */
  const double *huxxd = HUXX->data;      /* Assumes on a U grid */
  
  double *work_ud;
  double h, eta, tmp1, tmp2;

  if ( NOT_ALLOCATED( WORK_U ) ) {
    WORK_U = allocate_field2D_grid( NU, M, UGRID);
  }

  work_ud = WORK_U->data;

  for (i=0;i<NU;i++) {
    for (j=0;j<M;j++) {
      h = *hbxd++;
      eta = *etabxd++;
      tmp1 = a1*h*h*eta + (1.0/6.0) *eta*(h*h - eta*eta);
      tmp2 = a2*h*eta - 0.5*eta*(h+eta);
      *work_ud++ = tmp1 * ( *uxxd++ + *vxyd++) + tmp2 * ( *huxxd++ + *hvxyd++);
    }
  }

  field2D_diffx_add(E2, WORK_U, idx);  /* take the x gradient,w/ -idx,  of WORK_U_IE2 on U points to put onto ETA grid */

}



void b_calc_eta_E2PY(field2D *E2, const field2D *ETABY, const field2D *HBY, const field2D *VYY, const field2D *UXY, 
                     const field2D *HVYY, const field2D *HUXY)
{

  int NV = VYY->N;
  int NE2 = E2->N;
  int M = E2->M;
  int i,j;

  static field2D *WORK_V=NULL;   /* This is a V (n-2,m) grid workspace */
  //  static field2D *TE2=NULL;     /*  Temp workspace on a eta grid */

  const double *hbyd = &(HBY->data[M]);        /* Assumes on a V grid */
  const double *etabyd = ETABY->data;    /* Assumes on a V grid */
  const double *vyyd = &(VYY->data[0]);  /* Assumes on a V grid */
  const double *uxyd = UXY->data;        /* Assumes on a V grid */
  const double *hvyyd = &(HVYY->data[0]);      /* Assumes on a V grid */
  const double *huxyd = HUXY->data;      /* Assumes on a V grid */
  
  double *work_vd; 
  double h, eta, tmp1, tmp2;

  if ( NOT_ALLOCATED( WORK_V ) ) {
    WORK_V = allocate_field2D_grid( NE2, M, VGRID);
    //    TE2 = allocate_field2D_grid( NE2, M, ETAGRID);
  }



  work_vd = WORK_V->data;

  for (i=0;i<NV;i++) {
    for (j=0;j<M;j++) {
      h = *hbyd++;
      eta = *etabyd++;
      tmp1 = a1*h*h*eta + (1.0/6.0) *eta*(h*h - eta*eta);
      tmp2 = a2*h*eta - 0.5*eta*(h+eta);
      *work_vd++ = tmp1 * ( *uxyd++ + *vyyd++) + tmp2 * ( *huxyd++ + *hvyyd++);
    }
  }

                                                    /* Note -idy*gamma: gamma = 0 for weakly nonlinear case */
  field2D_diffy_v_to_eta_add(E2, WORK_V, -idy);   /* take the y gradient of WORK_U_IE2 on U points to put onto ETA grid */

  //  field2D_diffy_v_to_eta(TE2, WORK_V, -idy*gamma);   /* take the y gradient of WORK_U_IE2 on U points to put onto ETA grid */
  //  field2D_self_add(E2, TE2);                        /* E2 += TE2 */


}

/* Calculates the WEI_KIRBY E2 term, assumes that it will be added onto d(ETA)/dt */

void b_calc_eta_E2(field2D *E2, const field2D *ETABX, const field2D *ETABY,
                  const field2D *HBX, const field2D *HBY, const field2D *UXX, const field2D *VXY, const field2D *HUXX, const field2D *HVXY,
                  const field2D *UXY, const field2D *VYY, const field2D *HUXY, const field2D *HVYY)
{

   b_calc_eta_E2PX(E2, ETABX, HBX, UXX, VXY, HUXX, HVXY);  
   b_calc_eta_E2PY(E2, ETABY, HBY, VYY, UXY, HVYY, HUXY);


}



/* Used to calculate the initial condition of UUtmp and VVtmp */

void b_calc_UVtmp_from_UV(field2D *UUtmp, field2D *VVtmp, const field2D *HBX, const field2D *HBY, 
                const field2D *U, field2D *UXX, field2D *HUXX, const field2D *V, field2D *VYY, field2D *HVYY)
{

  field2D *U1;
  field2D *V1;

  U1 = allocate_field2D_grid(nx, ny, UGRID);
  V1 = allocate_field2D_grid(nx+1, ny, VGRID);
  field2D_set_to_value(V1, 0.0);

  b_calc_uxx(UXX, U);
  b_calc_huxx(HUXX, U, HBX);

  b_calc_vyy(VYY, V);
  b_calc_hvyy(HVYY, V, HBY);

  field2D_copy(UUtmp, U);
  field2D_copy(VVtmp, V);

  field2D_multiply(U1, UXX, HBX);          /* U1 = h^{bx} u_xx */
  field2D_mult_const(U1, b1);              /* U1 = b1 * h^{bx} u_xx */
  field2D_self_multiply(U1, HBX);          /* U1 = b1 * (h^{bx})^2 u_xx */
  field2D_self_add(UUtmp, U1);             /* UUtmp = U + b1 * (h^{bx})^2 u_xx */

  field2D_multiply(U1, HUXX, HBX);         /* U1 = h^{bx} (hu)_xx */
  field2D_mult_const(U1, b2);              /* U1 = b1 * h^{bx} u_xx */
  field2D_self_add(UUtmp, U1);             /* U1 = b1 * (h^{bx})^2 u_xx */

  field2D_apply_zero_bc(UUtmp);

  field2D_multiply_sub1(V1, VYY, HBY);          /* V1 = h^{by} v_yy */
  field2D_mult_const(V1, b1);              /* V1 = b1 * h^{by} v_yy */
  field2D_self_multiply(V1, HBY);          /* V1 = b1 * (h^{by})^2 v_yy */
  field2D_self_add(VVtmp, V1);             /* VVtmp = V + b1 * (h^{by})^2 v_yy */

  field2D_multiply_sub1(V1, HVYY, HBY);         /* V1 = h^{by} (hv)_yy */
  field2D_mult_const(V1, b2);              /* V1 = b1 * h^{by} v_yy */
  field2D_self_add(VVtmp, V1);             /* V1 = b1 * (h^{by})^2 v_yy */
  field2D_apply_nogradient_bc(VVtmp);

  //  deallocate_field2D(U1);
  //  deallocate_field2D(V1);
}