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



/* breaking.c   Falk Feddersen */       

#include <math.h>
#include "breaking.h"


field2D *BR_H;     /* just a pointer to the global H   - on a V grid */
field2D *BR_HBX;   /* just a pointer to the global HBX - on a U grid */
field2D *BR_HBY;   /* just a pointer to the global HBY - on a V grid */
field2D *BR_HBXY;  /* just a pointer to the allocated here HBXY - on a vorticity grid */

/* These are the total water depth variables.  Note that this should probably all be done
   in depth.c or something like that at each time step because other things like to have
   (h+\eta) on the various grids such as the d(eta)/dt  equation - and probably others */

field2D *BR_DEPTH;     /* the total depth (h+\eta) on eta grid */
field2D *BR_DEPTHBX;   /* the total global HBX - on a U grid */
field2D *BR_DEPTHBY;   /* just a pointer to the global HBY - on a V grid */
field2D *BR_DEPTHBXY;  /* just a pointer to the allocated here HBXY - on a vorticity grid */

field2D *BBB;         /* this is the breaking B parameter in Kennedy et al. */
field2D *BR_NU;       /* the breaking eddy viscosity which is on an eta grid */
field2D *BR_NUBXY;    /* the breaking eddy viscosity which is on an vort grid */

/* These are the variables from teh breaking_info_t *BR for initializtion */

switch_t breaking_flag; 
breaking_model_t  br_model;   /* The type of breaking model:  ZELT_BR, FUNWAVE_BR, FALK(other formulation not yet derived) */
switch_t active_wave_breaking_flag;  /* ON if wave breaking has occured since last diagnostic output */
double delta_b;
double delta_b2;  /* this = delta_b * delta_b */
double c_eta_tI;  
double c_eta_tF;
double cT;


void  init_breaking(breaking_info_t *BR, field2D_depth_var *DD)
{

  const int M = H->M;
  //  const int NETA = ETA->N;
  const int NHBX = DD->HBX->N;
  //  const int NH = DD->H->N;
  //  const int NV = V->N;

  /* First set up the breaking constants */

  breaking_flag = BR->breaking_flag;
  br_model = ZELT_BR;
  active_wave_breaking_flag = OFF;
  delta_b = BR->delta_b;
  delta_b2 = SQR( delta_b );
  c_eta_tI = BR->c_eta_tI;
  c_eta_tF = BR->c_eta_tF;
  cT = BR->cT;

  /* if breaking is turned off then return and don't allocate any more memory */

  if (breaking_flag == OFF) {
    return;
  }


  /* Now set up the points for breaking */

  BR_H = DD->H;
  BR_HBY = DD->HBY;
  BR_HBX = DD->HBX;
  BR_HBXY = DD->HBXY; 

  BR_DEPTH = DD->DEPTH;
  BR_DEPTHBY = DD->DEPTHBY;
  BR_DEPTHBX = DD->DEPTHBX;
  BR_DEPTHBXY = DD->DEPTHBXY; 

  /* allocate the eddy viscosity memory */

  BR_NU = allocate_field2D_grid(NHBX-1,M, ETAGRID);
  BR_NUBXY = allocate_field2D_grid(NHBX,M, VORTGRID);
  BBB = allocate_field2D_grid(NHBX-1,M, ETAGRID);

  field2D_set_to_value(BR_NU, 0.0);
  field2D_set_to_value(BR_NUBXY, 0.0); /* NOTE: at i=0 and N-1 this should always be zero. */

}

void  deallocate_breaking()
{

  deallocate_field2D(BR_NU);
  deallocate_field2D(BR_NUBXY);
  deallocate_field2D(BBB);

}

 /* This is for output purposes - returns a pointer to the eddy viscosity */

field2D*  breaking_eddy_viscosity()
{
  if (breaking_flag == OFF)
    return NULL;

  return BR_NU;
}

/* prints out breaking information */
void  breaking_report()
{
  if (breaking_flag == OFF)
    return;

  fprintf(stderr,"******* Wave Breaking Information **************\n");
  switch (br_model) {
  case ZELT_BR: fprintf(stderr,"Zelt type breaking used \n ");
    break;
  default: fprintf(stderr,"Unknown breaking model\n");
  }
  if (active_wave_breaking_flag == ON) {
    fprintf(stderr,"Wave breaking occurrred in this diagnostic period \n");
  }
  active_wave_breaking_flag = OFF;
  fprintf(stderr,"*************************************************\n");
}


    
void   update_breaking_parameter_zelt(field2D *BBB, const field2D *ETAT)
{

  double *bbbd = BBB->data;
  const double *etatd = ETAT->data;
  const int NUM_ETAT = ETAT->num;
  const double *depthbrd = BR_DEPTH->data;

  const double c_etat = 0.25;
  double etatc;              /* this is \eta_t^* = c_etat * sqrt(g * depth) */
  double etat;               /* local version of d(etat)/dt */
  int i;

  for (i=0;i<NUM_ETAT;i++) {
    etatc = c_etat * sqrt(gravity * (*depthbrd++) );
    etat = *etatd++;

    if (etat < etatc) {
      *bbbd++ = 0.0;
    }
    else { 
      if (etat < 2.0*etatc) {
	*bbbd++ = etat/etatc - 1.0;
	active_wave_breaking_flag = ON;
      }
      else {
	*bbbd++ = 1.0;
	active_wave_breaking_flag = ON;
      }
    }
  }


}


void  update_breaking_eddy_viscosity(const field2D *ETAT)
{

  /* First update BBB - B in funwave manual depending on the breaking style */

  switch(br_model) {
  case ZELT_BR:   update_breaking_parameter_zelt(BBB, ETAT);
    break;
  case FUNWAVE_BR: funwaveC_error_message("** FUNWAVE style breaking not implemented yet");
    break;
  default: funwaveC_error_message("** no other style breaking not implemented yet");
  }

  field2D_multiply(BR_NU, BR_DEPTH, ETAT);    /* BR_NU = (h+\eta)*d(eta)/dt */
  field2D_self_multiply(BR_NU, BBB);       /* BR_NU = B * (h+\eta)*d(eta)/dt -  BBB = B in Kennedy et al. */
  field2D_mult_const(BR_NU, delta_b2);     /* BR_NU = (\delta_b)^2 * B * (h+\eta)*d(eta)/dt : ** move into update_breaking_parameter()*/

  field2D_barxy_eta_to_vort(BR_NUBXY, BR_NU);  /* create the breaking eddy viscosity on the VORT grid */

}

/* This function gets called to add to the RHS of the momentum equations the eddy-viscosity type
   breaking terms.   Also may need the depth terms.   Note that  FUNWAVE does (1/h) ( div [ (nu*h) grad(u) ] )
   or something like that - it's not exactly right.   I'm gonna try  just doing  div ( nu grad(u) ) for u etc. 

  dU/dt +=  .....  & dV/dt += .....
*/


void  rhs_uvmom_breaking_add(field2D *UTBR, field2D *VTBR, const field2D *U, const field2D *V, const field2D *ETAT)
{

  const int M = U->M;
  const int NETA = ETA->N;
  const int NU = U->N;
  const int NV = V->N;

  static field2D *DUDX=NULL;
  static field2D *DUDY=NULL;
  static field2D *DVDX=NULL;
  static field2D *DVDY=NULL;
  static field2D *TMPU=NULL, *TMPV=NULL;
  static field2D *BR_U=NULL,  *BR_V=NULL;  


  /* If breaking dynamics are turned off then simply return */

  if (breaking_flag == OFF) {
    return;
  }


  if ( NOT_ALLOCATED(DUDX) ) {   /* what are the sizes of this??  */
    DUDX = allocate_field2D_grid(NETA, M, ETAGRID);
    DUDY = allocate_field2D_grid(NU, M, VORTGRID);
    DVDX = allocate_field2D_grid(NV-1, M, VORTGRID);
    DVDY = allocate_field2D_grid(NETA, M, ETAGRID);
    BR_U = allocate_field2D_grid(NU-2, M, UTGRID);
    TMPU = allocate_field2D_grid(NU-2, M, UTGRID);
    BR_V = allocate_field2D_grid(NV-2, M, VTGRID);
    TMPV = allocate_field2D_grid(NV-2, M, VTGRID);
  }


  update_breaking_eddy_viscosity(ETAT);


  /* BR_U = \delta_x ( h*nu \delta_x U ) */

  field2D_diffx(DUDX, U, idx);            /* DUDX is on an eta point */
  field2D_self_multiply_two(DUDX, BR_DEPTH, BR_NU);  /* multiply it by the depth (on eta point) */
  //  field2D_self_multiply(DUDX, BR_NU);     /* multiply it by the eddy viscosity */
  field2D_diffx(BR_U, DUDX, idx);         /* BR_U is now on U points again (NX-2,M) */  

  /* += \delta_y ( h*nu \delta_y U ) */

  field2D_diffy_u_to_vort(DUDY, U, idy);  /* DUDY (N,M) */  
  field2D_self_multiply_two(DUDY, BR_DEPTHBXY, BR_NUBXY);    /* multiply it by the depth (on eta point) */
  //  field2D_self_multiply(DUDY, BR_NUBXY);    /* multiply it by the eddy viscosity */
  field2D_diffy_vort_to_u(TMPU, DUDY, idy);     /* BR_U is now on U points again  */  
  field2D_self_add(BR_U, TMPU);                /* BR_U +=  d/dy ( h*nu dU/dy )   */  

  /* BR_U *= 1/HBX */

  field2D_self_divide(BR_U, BR_DEPTHBX);  /* BR_U = (1/HBX) (d/dx( h*nu du/dx) + d/dy( h*nu du/dy))  */  

  /* This is the V momentum parts */

  /* BR_V = \delta_x ( h*nu \delta_x V ) */

  field2D_diffx(DVDX, V, idx);
  field2D_self_multiply_two(DVDX, BR_DEPTHBXY, BR_NUBXY);   /* multiply it by the depth (on eta point) */
  //  field2D_self_multiply(DVDX, BR_NUBXY);    /* multiply it by the eddy viscosity */
  field2D_diffx(BR_V, DVDX, idx);           /* BR_V is on V points again  */  

  field2D_diffy_v_to_eta(DVDY, V, idy);        /* DVDY is on an eta point */
  field2D_self_multiply_two(DVDY, BR_DEPTH, BR_NU);       /* multiply by the depth (on eta point) */
  //  field2D_self_multiply(DVDY, BR_NU);          /* multiply it by the eddy viscosity */
  field2D_diffy_eta_to_v(TMPV, DVDY, idy);     /* BR_V is now on U points again  */  
  field2D_self_add(BR_V, TMPV);                /* BR_V *= (1/HBX)  */  

  /* BR_U *= 1/HBY */

  field2D_self_divide(BR_V, BR_DEPTHBY);  /* BR_U *= (1/HBX) (d/dx( h*nu du/dx) + d/dy( h*nu du/dy))  */  

  
 /* Now add to the RHS of the momentum equations */

  field2D_self_add(UTBR, BR_U);
  field2D_self_add(VTBR, BR_V);



}