bousinessq_dynamics.c

波浪数值模拟

/*
 * Copyright (c) 2001-2005 Falk Feddersen
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* ---  bousinessq_dynamics.c  */

#include "bousinessq_dynamics.h"
#include "bousinessq_operations.h"   /* for b_calc_ux_plus_vy */


/* #ifdef LINEAR */
/* const double beta = -0.531;         /\*  betea =  za/h *\/ */
/* const double a1 = 0; //0.5*beta*beta - (1.0/6.0);   /\* a1 and a2 are bousinessq parameters in the d(eta)/dt = ... equation *\/ */
/* const double a2 = 0; //beta + 0.5; */
/* const double b1 = 0.0; //0.5* beta*beta;        /\* b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation *\/ */
/* const double b2 = 0.0; //beta; */
/* const double gamma = 0.0; */
/* #endif */

/* #ifdef NSWE */
/* const double beta = -0.531;         /\*  betea =  za/h *\/ */
/* const double a1 = 0; //0.5*beta*beta - (1.0/6.0);   /\* a1 and a2 are bousinessq parameters in the d(eta)/dt = ... equation *\/ */
/* const double a2 = 0; //beta + 0.5; */
/* const double b1 = 0.0; //0.5* beta*beta;        /\* b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation *\/ */
/* const double b2 = 0.0; //beta; */
/* const double gamma = 0.0; */
/* #endif */

/* #ifdef PEREGRINE   /\* peregrine depth-averaged dynamics *\/ */
/* const double beta = -0.5; */
/* const double a1 = 0.0;  */
/* const double a2 = 0.0;  */
/* const double b1 = 1.0/6.0; /\* not sure why b1 is what it is *\/ */
/* const double b2 = -0.5;    /\* or b2 but these are the peregrine JFM (1967) values *\/ */
/* const double gamma = 0.0; */
/* #endif  /\* PEREGRINE *\/ */

/* #ifdef WEI_KIRBY */
/* const double beta = -0.531;         /\*  betea =  za/h *\/ */
/* const double a1 =  -2.568616666666665e-02;  /\* a1 = 0.5*beta_h *beta_h - (1.0/6.0);  a1 and a2 are bous parameters in d(eta)/dt =  *\/ */
/* const double a2 = -0.031;                     /\* a2 = beta + 0.5;  *\/ */
/* const double b1 =   1.409805000000000e-01;    /\* b1 = 0.5* beta_h*beta_h;   b1 and b2 are bousinessq parameters in the d(u,v)/dt = *\/ */
/* const double b2 =  -0.531; */
/* const double gamma = 1.0; */
/* #endif */

/* #ifdef NWOGU */
/* const double beta = -0.531;         /\*  betea =  za/h *\/ */
/* const double a1 =  -2.568616666666665e-02;    /\* a1 = 0.5*beta_h *beta_h - (1.0/6.0);  a1 and a2 are bous parameters in d(eta)/dt = ... equation *\/ */
/* const double a2 = -0.031;                     /\* a2 = beta + 0.5;  *\/ */
/* const double b1 =   1.409805000000000e-01;    /\* b1 = 0.5* beta_h*beta_h;   b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation *\/ */
/* const double b2 =  -0.531; */
/* const double gamma = 0.0; */
/* #endif */


double beta;         /*  bete =  za/h */
double a1;
double a2;
double b1;
double b2;
//double gamma;

/* ---- Constant Initialization Routine --------- */

void b_init_constants(dynamics_t Dynamics)
{

  switch (Dynamics) {
  case WEI_KIRBY_DYNAMICS:
    beta = -0.531;         /*  betea =  za/h */
    a1 =  -2.568616666666665e-02;  /* a1 = 0.5*beta_h *beta_h - (1.0/6.0);  a1 and a2 are bous parameters in d(eta)/dt =  */
    a2 = -0.031;                     /* a2 = beta + 0.5;  */
    b1 =   1.409805000000000e-01;    /* b1 = 0.5* beta_h*beta_h;   b1 and b2 are bousinessq parameters in the d(u,v)/dt = */
    b2 =  -0.531;
    //    gamma = 1.0;
    break;
  case NWOGU_DYNAMICS:
    beta = -0.531;         /*  betea =  za/h */
    a1 =  -2.568616666666665e-02;    /* a1 = 0.5*beta_h *beta_h - (1.0/6.0);  a1 and a2 are bous parameters in d(eta)/dt = ... equation */
    a2 = -0.031;                     /* a2 = beta + 0.5;  */
    b1 =   1.409805000000000e-01;    /* b1 = 0.5* beta_h*beta_h;   b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation */
    b2 =  -0.531;
    //    gamma = 0.0;
    break;
  case PEREGRINE_DYNAMICS:
    beta = -0.5;
    a1 = 0.0; 
    a2 = 0.0; 
    b1 = 1.0/6.0; /* not sure why b1 is what it is */
    b2 = -0.5;    /* or b2 but these are the peregrine JFM (1967) values */
    //    gamma = 0.0;
    break;
  case NSWE_DYNAMICS:
    beta = -0.531;         /*  betea =  za/h */
    a1 = 0; //0.5*beta*beta - (1.0/6.0);   /* a1 and a2 are bousinessq parameters in the d(eta)/dt = ... equation */
    a2 = 0; //beta + 0.5;
    b1 = 0.0; //0.5* beta*beta;        /* b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation */
    b2 = 0.0; //beta;
    //    gamma = 0.0;
    break;
  case LINEAR_DYNAMICS:
    beta = -0.531;         /*  betea =  za/h */
    a1 = 0; //0.5*beta*beta - (1.0/6.0);   /* a1 and a2 are bousinessq parameters in the d(eta)/dt = ... equation */
    a2 = 0; //beta + 0.5;
    b1 = 0.0; //0.5* beta*beta;        /* b1 and b2 are bousinessq parameters in the d(u,v)/dt = ... equation */
    b2 = 0.0; //beta;
    //    gamma = 0.0;
    break;
  }


}
/* ----------  NUMERICAL MODEL BOUSINESSQ DYNAMICS ROUTINUES ------------------- */



void b_calc_u_nonlinear(field2D *BUNL, const field2D *U, const field2D *V)
{

  int M = V->M;
  int MN = M-1;
  int N_BUNL = BUNL->N;

  const double *ud = &(U->data[M]);
  const double *vd = &(V->data[M]);
  double *unld = BUNL->data;

  const double hidx = 0.5*idx;
  const double hidy = 0.5*idy;

  double nluux, nlvuy, vbxy;

  int i,j;

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

    /* First do j=0 */

    /* Here go the nonlinear terms 1) u \delta_x (u^{bx}) +  v^{bxy} \delta_y  */

    nluux = hidx * ( *ud * ( *(ud+M) - *(ud-M) ) ) ;
    vbxy =  0.25 * (*vd  + *(vd+M) + *(vd+MN) + *(vd+MN+M) );   // DR2(V,i,j) + DR2(V,i,jm) + DR2(V,i+1,j) + DR2(V,i+1,jm) 
    nlvuy = hidy * ( vbxy * ( *(ud+1) - *(ud+MN) ) ) ;

    *unld++ = -(nluux + nlvuy);
    ud++;
    vd++;

       /* Done with j=0 */

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

      /* Here go the nonlinear terms 1) u \delta_x (u^{bx}) +  v^{bxy} \delta_y  */
 
      nluux = hidx * ( *ud * ( *(ud+M) - *(ud-M) ) ) ;
      vbxy =  0.25 * (*vd  + *(vd+M) + *(vd-1) + *(vd+M-1) );   // DR2(V,i,j) + DR2(V,i,jm) + DR2(V,i+1,j) + DR2(V,i+1,jm) 
      nlvuy = hidy * ( vbxy * ( *(ud+1) - *(ud-1) ) ) ;

      *unld++ = -(nluux + nlvuy);
      ud++;
      vd++;
    }

    /* Now do j=M-1  (ny-1) */

    /* Here go the nonlinear terms 1) u \delta_x (u^{bx}) +  v^{bxy} \delta_y  */

    nluux = hidx * ( *ud * ( *(ud+M) - *(ud-M) ) ) ;
    vbxy =  0.25 * (*vd  + *(vd+M) + *(vd-1) + *(vd+M-1) );   // DR2(V,i,j) + DR2(V,i,jm) + DR2(V,i+1,j) + DR2(V,i+1,jm) 
    nlvuy = hidy * ( vbxy * ( *(ud-MN) - *(ud-1) ) ) ;

    *unld++ = -(nluux + nlvuy);
    ud++;
    vd++;

  }

}




void b_calc_v_nonlinear(field2D *BVNL, const field2D *U, const field2D *V)
{

  int M = V->M;
  int MN = M-1;
  int N_BVNL = BVNL->N;

  const double *ud = &(U->data[M]);
  const double *vd = &(V->data[M]);
  double *vnld = BVNL->data;

  double nluvx, nlvvy, ubxy;
  const double hidx = 0.5*idx;
  const double hidy = 0.5*idy;

  int i,j;

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

    /* First do j=0 */

    /* Here go the nonlinear terms 1) u \delta_x (u^{bx}) +  v^{bxy} \delta_y  */

    nlvvy = hidy * ( *vd * ( *(vd+1) - *(vd+MN) ) ) ;
    ubxy =  0.25 * (*ud  + *(ud-M) + *(ud+1) + *(ud-M+1) );   // DR2(V,i,j) + DR2(V,i,jm) + DR2(V,i+1,j) + DR2(V,i+1,jm) 
    nluvx = hidx * ( ubxy * ( *(vd+M) - *(vd-M) ) ) ;

    *vnld++ = -(nluvx + nlvvy);
    ud++;
    vd++;

       /* Done with j=0 */

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

      /* Here go the nonlinear terms 1) u \delta_x (u^{bx}) +  v^{bxy} \delta_y  */
 
      nlvvy = hidy * ( *vd * ( *(vd+1) - *(vd-1) ) ) ;
      ubxy =  0.25 * (*ud  + *(ud-M) + *(ud+1) + *(ud-M+1) );   // DR2(V,i,j) + DR2(V,i,jm) + DR2(V,i+1,j) + DR2(V,i+1,jm) 
      nluvx = hidx * ( ubxy * ( *(vd+M) - *(vd-M) ) ) ;
  
      *vnld++ = -(nluvx + nlvvy);
      ud++;
      vd++;
    }

    /* Now do j=M-1 (ny-1) */

    /* Here go the nonlinear terms 1) u \delta_x (u^{bx}) +  v^{bxy} \delta_y  */

    nlvvy = hidy * ( *vd * ( *(vd+1) - *(vd-1) ) ) ;
    ubxy =  0.25 * (*ud  + *(ud-M) + *(ud-MN) + *(ud-M-MN) );   // DR2(V,i,j) + DR2(V,i,jm) + DR2(V,i+1,j) + DR2(V,i+1,jm) 
    nluvx = hidx * ( ubxy * ( *(vd+M) - *(vd-M) ) ) ;

    *vnld++ = -(nluvx + nlvvy);
    ud++;
    vd++;

  }

}


/* F1 = output, HBX = h^{bx}, HBY = h^{by}, 

   Calculates the term F1 from the FUNWAVE manual:  F1 = -h[ b1*h vt_{xy} + b2 (h*vt)_{xy} ] 
   where vt is the time derivative of v.     The numerics of this are to be on U points so
   
   F1 = - h^{bx} * [ b1 h^{bx} \delta_y \delta_x (vt) +  b2 * \delta_y \delta_x (h^{by} * vt)

   note that field2D dV/dt should have same size as V 

*/


void b_calc_umom_F1(field2D  *F1, const field2D *HBX, const field2D *V, const field2D *HV) 
{

  int i,j;
  int M = V->M;
  int MN = M-1;
  int NF = F1->N;

  const double *vd = REF2(V, 1);    // &(V->data[M]);    /* nx+1,ny */
  const double *hvd = REF2(HV, 1);  //&(HV->data[M]);  /* nx+1,ny */
  const double *hbxd = REF2(HBX, 1);   //&(HBX->data[M]);  /* nx, ny */

  double *f1d = &(F1->data[0]);        /* nx-2, ny */

  double dvdxdy;
  double dhvdxdy;
  double hbx;


  /* Now go through and calculate  F1 */

  for (i=0;i<NF;i++) {
    /* j = 0 */

    hbx = *hbxd++;
    dvdxdy = idx * idy * ( (*(vd+M)- *vd) - (*(vd+M+MN) - *(vd+MN)) );
    dhvdxdy = idx * idy * ( (*(hvd+M)- *hvd) - (*(hvd+M+MN) - *(hvd+MN)) );
    *f1d++  = -hbx * (b1 * hbx *dvdxdy + b2 * dhvdxdy);
    vd++;
    hvd++;

    for (j=1;j<M;j++) {
      hbx = *hbxd++;
      dvdxdy = idx * idy * ( (*(vd+M)- *vd) - (*(vd+M-1) - *(vd-1)) );
      dhvdxdy = idx * idy * ( (*(hvd+M)- *hvd) - (*(hvd+M-1) - *(hvd-1)) );
      *f1d++  = -hbx * (b1 * hbx *dvdxdy + b2 * dhvdxdy);
      vd++;
      hvd++;
    }

  }


}


/* G1 = output, HBX = h^{bx}, HBY = h^{by}, 

   Calculates the term G1 from the FUNWAVE manual:  G1 = -h[ b1*h ut_{xy} + b2 (h*ut)_{xy} ] 
   where ut is the time derivative of u.     The numerics of this are to be on V points so
   
   G1 = - h^{by} * [ b1 h^{by} \delta_y \delta_x (u) +  b2 * \delta_y \delta_x (h^{bx} * u)

   note that field2D U  has size of U

*/


void b_calc_vmom_G1(field2D  *G1, const field2D *HBY, const field2D *U, const field2D *HU )
{


  int i,j;
  int M = U->M;
  int MN = M-1;
  int NG = G1->N;

  const double *ud = &(U->data[M]);
  const double *hbyd = &(HBY->data[M]);
  const double *hud =&(HU->data[M]);

  double *g1d = &(G1->data[0]);

  double dudxdy;
  double dhudxdy;
  double hby;


  /* Now go through and calculate  G1t */

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

    for (j=0;j<M-1;j++) {
      hby = *hbyd++;
      dudxdy = idx * idy * (  (*(ud+1) - *(ud-M+1)) - (*ud - *(ud-M)));
      dhudxdy = idx * idy * (  (*(hud+1) - *(hud-M+1))- (*hud - *(hud-M)) );
      *g1d++  = -hby * (b1 * hby *dudxdy + b2 * dhudxdy);
      ud++;
      hud++;
    }

    /* now j= ny-1 */

    hby = *hbyd++;
    dudxdy = idx * idy * ( (*(ud-MN) - *(ud-M-MN)) - (*ud - *(ud-M))   );
    dhudxdy = idx * idy * ( (*(hud-MN) - *(hud-M-MN)) - (*hud - *(hud-M))   );
    *g1d++  = -hby * (b1 * hby *dudxdy + b2 * dhudxdy);
    ud++;
    hud++;

  }


}


/* Wrapper to calculate F1t and G1t individually */

void b_calc_uvmom_F1_G1(field2D  *F1, field2D *G1, const field2D *HBX, const field2D *HBY, const field2D *U, const field2D *V, const field2D *HU, const field2D *HV)
{

   b_calc_umom_F1(F1, HBX, V, HV);
   b_calc_vmom_G1(G1, HBY, U, HU);
}