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



/* sponge.c   Falk Feddersen */       

#include <math.h>

#include "sponge.h"


double *WS1_X0_U;   /* the Upts x=0 boundary sponge terms for ws1 in funwave manual */
double *WS2_X0_U;   /* the Upts x=0 boundary sponge terms for ws2 in funwave manual */
double *WS3_X0_U;   /* the Upts x=0 boundary sponge terms for ws3 in funwave manual */

double *WS1_X0_V;   /* the Vpts x=0 boundary sponge terms for ws1 in funwave manual */
double *WS2_X0_V;   /* the Vpts x=0 boundary sponge terms for ws2 in funwave manual */
double *WS3_X0_V;   /* the Vpts x=0 boundary sponge terms for ws3 in funwave manual */

double *WS1_XL_U;   /* the Upts x=L boundary sponge terms for ws1 in funwave manual */
double *WS2_XL_U;   /* the Upts x=L boundary sponge terms for ws2 in funwave manual */
double *WS3_XL_U;   /* the Upts x=L boundary sponge terms for ws3 in funwave manual */

double *WS1_XL_V;   /* the Vpts x=L boundary sponge terms for ws1 in funwave manual */
double *WS2_XL_V;   /* the Vpts x=L boundary sponge terms for ws2 in funwave manual */
double *WS3_XL_V;   /* the Vpts x=L boundary sponge terms for ws3 in funwave manual */


switch_t  sponge_layer_flag;   // sponge layer either on/off
int ix0_width;                // width of the sponge layer starting at x=0 (the shore)
int ixL_width;               // width of the sponge layer starting at x=Lx (the far boundary)
double cspg1;                 // factor that goes in front of Wspg1
double cspg2;                 // factor that goes in front of Wspg1
double cspg3;                 // factor that goes in front of Wspg1
double frequency;             // this is the wave frequency
double omega;


void allocate_all_sponge_layer_arrays() 
{

  /* allocate all the arrays */

  WS1_X0_U = (double * ) g_malloc(sizeof(double)*ix0_width);
  WS2_X0_U = (double * ) g_malloc(sizeof(double)*ix0_width);
  WS3_X0_U = (double * ) g_malloc(sizeof(double)*ix0_width);

  WS1_X0_V = (double * ) g_malloc(sizeof(double)*ix0_width);
  WS2_X0_V = (double * ) g_malloc(sizeof(double)*ix0_width);
  WS3_X0_V = (double * ) g_malloc(sizeof(double)*ix0_width);

  WS1_XL_U = (double * ) g_malloc(sizeof(double)*ixL_width);
  WS2_XL_U = (double * ) g_malloc(sizeof(double)*ixL_width);
  WS3_XL_U = (double * ) g_malloc(sizeof(double)*ixL_width);

  WS1_XL_V = (double * ) g_malloc(sizeof(double)*ixL_width);
  WS2_XL_V = (double * ) g_malloc(sizeof(double)*ixL_width);
  WS3_XL_V = (double * ) g_malloc(sizeof(double)*ixL_width);

  //  SP_F_X0 = (double * ) g_malloc(sizeof(double)*ix0_width);    
  //  SP_F_XL = (double * ) g_malloc(sizeof(double)*ixL_width);    
}


void  init_sponge_layer(sponge_layer_info_t *SP, int nx, int ny)
{
  int i;
  double f_u, f_v;
  double x_u, x_v, dx_u, dx_v, tmp;
  double xwidth_x0, xwidth_xL;
  double xs_x0, xl_x0;
  double xs_xL, xl_xL;

  sponge_layer_flag = SP->sponge_layer_flag;

  if (sponge_layer_flag == OFF) /* if the flag is off then don't do anything! */
    return;                        

  ix0_width = SP->ix0_width;
  ixL_width = SP->ixL_width;
  cspg1 = SP->cspg1;
  cspg2 = SP->cspg2;
  cspg3 = SP->cspg3;


  
  omega = 2.0 * pi * SP->frequency;

  allocate_all_sponge_layer_arrays();

  /* set up x0 sponge stuff first */

  xs_x0 = ix0_width*dx;
  xwidth_x0 = ix0_width*dx;
  xl_x0 = 0;

  for (i=0;i<ix0_width;i++) {
    x_u = (i+1)*dx;
    x_v = (i+0.5)*dx;
    dx_u = (x_u - xs_x0);
    dx_v = (x_v - xs_x0);
    tmp = exp(1)-1;
    f_u = ( exp(  dx_u*dx_u / (xwidth_x0 * xwidth_x0) ) - 1.0 )/tmp;
    WS1_X0_U[i] = cspg1*omega*f_u;
    WS2_X0_U[i] = cspg2*omega*f_u;
    WS3_X0_U[i] = cspg3*omega*f_u;

    f_v = ( exp(  dx_v*dx_v / (xwidth_x0 * xwidth_x0) ) - 1.0 )/tmp;

    WS1_X0_V[i] = cspg1*omega*f_v;
    WS2_X0_V[i] = cspg2*omega*f_v;
    WS3_X0_V[i] = cspg3*omega*f_v;
  }

  /* set up xL sponge stuff next */

  xs_xL = ixL_width*dx;
  xwidth_xL = ixL_width*dx;
  xl_xL = 0;


  for (i=0;i<ixL_width;i++) {
    x_u = (i+1)*dx;
    x_v = (i+0.5)*dx;
    dx_u = (x_u - xs_xL);
    dx_v = (x_v - xs_xL);
    tmp = exp(1)-1;
    f_u = ( exp(  dx_u*dx_u / (xwidth_xL * xwidth_xL) ) - 1.0 )/tmp;
    WS1_XL_U[ixL_width-i-1] = cspg1*omega*f_u;  /* reverse the order of application */
    WS2_XL_U[ixL_width-i-1] = cspg2*omega*f_u;
    WS3_XL_U[ixL_width-i-1] = cspg3*omega*f_u;

    f_v = ( exp(  dx_v*dx_v / (xwidth_xL * xwidth_xL) ) - 1.0 )/tmp;
    WS1_XL_V[ixL_width-i-1] = cspg1*omega*f_v;
    WS2_XL_V[ixL_width-i-1] = cspg2*omega*f_v;
    WS3_XL_V[ixL_width-i-1] = cspg3*omega*f_v;
  }

  f_v=1.0;

}

void  deallocate_sponge_layer( )
{
  g_free(WS1_X0_U);
  g_free(WS2_X0_U);
  g_free(WS3_X0_U);
  g_free(WS1_X0_V);
  g_free(WS2_X0_V);
  g_free(WS3_X0_V);

  g_free(WS1_XL_U);
  g_free(WS2_XL_U);
  g_free(WS3_XL_U);
  g_free(WS1_XL_V);
  g_free(WS2_XL_V);
  g_free(WS3_XL_V);

  //  g_free(SP_F_X0);
  //  g_free(SP_F_XL);

    
}


void  rhs_uvmom_sponge_layer_add(field2D *UTSP, field2D *VTSP, const field2D *U, const field2D *V)
     //                                 const field2D *LAP_U, const field2D *LAP_V)
{
  int i,j;
  int ixLstart_u, ixLstart_v;

  const int M = UTSP->M;

  const double *ux0d, *vx0d, *uxLd, *vxLd;

  double *utspx0d, *utspxLd;
  double *vtspx0d, *vtspxLd;

  double ws1u, ws2u, ws3u;
  double ws1v, ws2v, ws3v;

  /* Zero: check to see if sponge layer is on - if not return */

  if (sponge_layer_flag == OFF)
    return;

  /* First the shoreward sponge layer */

  utspx0d = &(UTSP->data[M]);
  vtspx0d = &(VTSP->data[M]);
  ux0d  = &(U->data[M]);
  vx0d  = &(V->data[M]);


  for (i=0;i<ix0_width;i++) {
    ws1u = -WS1_X0_U[i];
    ws1v = -WS1_X0_V[i];

    for (j=0;j<M;j++) {
      *utspx0d++ += ws1u * *ux0d++;
      *vtspx0d++ += ws1v * *vx0d++;
    }
  }

  /* Next the seaward xL sponge layer */

  ixLstart_u = U->N - ixL_width - 1;
  ixLstart_v = V->N - ixL_width - 1;

  utspxLd = &(UTSP->data[M*(ixLstart_u)]);   //   utspxLd = &(UTSP->data[M*(ixLstart_u-1)]);
  vtspxLd = &(VTSP->data[M*(ixLstart_v)]);   //   vtspxLd = &(VTSP->data[M*(ixLstart_v-1)]);
  uxLd  = &(U->data[M*ixLstart_u]);
  vxLd  = &(V->data[M*ixLstart_v]);

  for (i=0;i<ixL_width;i++) {
    ws1u = -WS1_XL_U[i];
    ws1v = -WS1_XL_V[i];

    for (j=0;j<M;j++) {
      *utspxLd++ += ws1u * *uxLd++;
      *vtspxLd++ += ws1v * *vxLd++;
    }
  }


  /* First do the seaward sponge layer */





}




/* report on bottom stress stability parameters by doing a linearlized rayliegh friction analysis */

//#define SMAX(a,b) ( (a)>(b) ) ? (a) : (b) )

void  sponge_layer_stability_report( double dt )
{

    // find the largest S_u, find smallest h  
 
  double S_sp1, S_sp2;
  double max_Ws1_x0, max_Ws1_xL; 

  if ( sponge_layer_flag == OFF )
    return;

  max_Ws1_x0 = WS1_X0_U[0];
  max_Ws1_xL = WS1_XL_U[ixL_width-1];


  S_sp1 = fabs( max_Ws1_x0 *dt);     // < 0.96 (AB2)
  S_sp2 = fabs( max_Ws1_xL *dt);     // < 0.96 (AB2)

  fprintf(stderr,"Sponge Rayleigh Friction (x0,xL) # = (%g , %g) \n ",S_sp1, S_sp2);

  /* note that this limit of 0.25 is arbitary based on some limited experience.
     It is definitely different than what is used in shallow w/ AB2 time stepping (0.96) */

  if (( S_sp1>0.25) || ( S_sp2>0.25)) 
    fprintf(stderr,"** Warning: sponge layer may be Rayleigh Friction unstable, S_sp>0.25! \n");
  else
    fprintf(stderr,"\n");
  
    
}



/* Calculates \lap(U) and \lap(V) from istart to iend over all j.  Requires teh LU and LV be
   of the right dimensions.   Maybe multiply it by the sponge layer viscosity as well */


void sponge_layer_laplacian_UV(field2D *LU, field2D *LV, const field2D *U, const field2D *V, int istart, double *ws2U, double *ws2V)
{

  int i,j,im,M,MN;
  int NLU;
  double utmp0,utmp1,utmp2;
  double vtmp0,vtmp1,vtmp2;
  double *lud, *lvd, *lvd0, *lvd1, *lvd2;
  double *ud, *udup, *uddn;
  double *vd, *vdup, *vddn;

  double nu_sponge_u;
  double nu_sponge_v;

  /* First check to make sure that the size of LU and LV are appropriate given istgart and iend */

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

  lud = LU->data;
  lvd = LV->data;

  ud =  REF2(U, istart);
  udup = REF2(U, istart+1);
  uddn = REF2(U, istart-1);

  vd =  REF2(V, istart);
  vdup = REF2(V, istart+1);
  vddn = REF2(V, istart-1);

  for (i=0;i<NLU;i++) {
    nu_sponge_u = *ws2U++;    // this is the viscosity 
    nu_sponge_v = *ws2V++;    // this is the viscosity 
    utmp0 = (*ud) * 2.0;  // here we do j=0
    utmp1 =  *udup++ - utmp0 + *uddn++;
    utmp2 = *(ud+1) - utmp0 + *(ud+MN);
    *lud++ = nu_sponge_u * (idx2*utmp1 + idy2*utmp2);
    vtmp0 = (*vd) * 2.0;  // here we do j=0
    vtmp1 =  *vdup++ - vtmp0 + *vddn++;
    vtmp2 = *(vd+1) - vtmp0 + *(vd+MN);
    *lvd++ = nu_sponge_v * (idx2*vtmp1 + idy2*vtmp2);
    vd++;  
    ud++; 
    for (j=1;j<M-1;j++) {
      utmp0 = (*ud) * 2.0;
      utmp1 =  *udup++ - utmp0 + *uddn++;
      utmp2 = *(ud+1) - utmp0 + *(ud-1);
      *lud++ = nu_sponge_u * (idx2*utmp1 + idy2*utmp2);
      ud++; 
      vtmp0 = (*vd) * 2.0;
      vtmp1 =  *vdup++ - vtmp0 + *vddn++;
      vtmp2 = *(vd+1) - vtmp0 + *(vd-1);
      *lvd++ = nu_sponge_v * (idx2*vtmp1 + idy2*vtmp2);
      vd++; 
    }
    utmp0 = (*ud) * 2.0;  // here we do j=ny-1
    utmp1 =  *udup++ - utmp0 + *uddn++;
    utmp2 = *(ud-MN) - utmp0 + *(ud-1);
    *lud++ = nu_sponge_u * (idx2*utmp1 + idy2*utmp2);
    ud++;  
    vtmp0 = (*vd) * 2.0;  // here we do j=ny-1
    vtmp1 =  *vdup++ - vtmp0 + *vddn++;
    vtmp2 = *(vd-MN) - vtmp0 + *(vd-1);
    *lvd++ = nu_sponge_v * (idx2*vtmp1 + idy2*vtmp2);
    vd++; 
  }

  /* lastly we must do i=nx for the V laplacian! */
  nu_sponge_v = *ws2V++;    // this is the viscosity 

  vtmp0 = (*vd) * 2.0;  // here we do j=0
  vtmp1 =  *vdup++ - vtmp0 + *vddn++;
  vtmp2 = *(vd+1) - vtmp0 + *(vd+MN);
  *lvd++ = nu_sponge_v * (idx2*vtmp1 + idy2*vtmp2);
  vd++;  
  for (j=1;j<M-1;j++) {
    vtmp0 = (*vd) * 2.0;
    vtmp1 =  *vdup++ - vtmp0 + *vddn++;
    vtmp2 = *(vd+1) - vtmp0 + *(vd-1);
    *lvd++ = nu_sponge_v * (idx2*vtmp1 + idy2*vtmp2);
    vd++;  
    }
  vtmp0 = (*vd) * 2.0;  // here we do j=ny-1
  vtmp1 =  *vdup++ - vtmp0 + *vddn++;
  vtmp2 = *(vd-MN) - vtmp0 + *(vd-1);
  *lvd++ = (idx2*vtmp1 + idy2*vtmp2);
  vd++; 



}