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



/* lateral_mixing.c   Falk Feddersen */       

/* Includes both lateral friction and bottom stress routines here 


      Functions to calculate various newtonian and biharmonic friction
      for the shallow water equations

      This is done using pointer arithmetic to make things faster

*/

#include <math.h>

#include "lateral_mixing.h"


field2D *LUU, *LVV;   /* this will crash without allocation */


double nu_bi;      // this is bi-harmonic  m^4/s  
double nu_newt;    // Newtonian friction term m^2/s

mix_t lateral_mixing_type;



void laplacian0_all(field2D *LU, field2D *LV, const field2D *U, const field2D *V);


void newtonian0_u_all(field2D *NU, const field2D *U);
void newtonian0_v_all(field2D *NV, const field2D *V);

//void newtonian1_u_all(field2D *NU)
//void newtonian1_v_all(field2D *NV);

//void newtonian2_u_all(field2D *NU);
//void newtonian2_v_all(field2D *NV);

void bi_laplacian0_all(field2D *LXX, const field2D *XX);

//void biharmonic0_u_all(field2D *BU, const field2D *LU);
//void biharmonic0_v_all(field2D *BV, const field2D *LV);

void mult_laplacian0_depth(field2D *LU, field2D *LV);

//void biharmonic1_u_all(field2D *BU, field2D *LU);
//void biharmonic1_v_all(field2D *BV, field2D *LV);



/* here the BC are  lap(U) = 0 at boundaries 
   and if slip then  lap(V)_x = 0 
   no-slip then  lap(V) = 0   at BC   */



void laplacian0_all(field2D *LU, field2D *LV, const field2D *U, const field2D *V)
{

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

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

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

  ud = &(U->data[M]);
  udup = &(U->data[2*M]);
  uddn = U->data;

  vd = &(V->data[M]);
  vdup = &(V->data[2*M]);
  vddn = V->data;

  for (i=1;i<nx-1;i++) {
    utmp0 = (*ud) * 2.0;  // here we do j=0
    utmp1 =  *udup++ - utmp0 + *uddn++;
    utmp2 = *(ud+1) - utmp0 + *(ud+MN);
    *lud++ = (idx2*utmp1 + idy2*utmp2);
    vtmp0 = (*vd) * 2.0;  // here we do j=0
    vtmp1 =  *vdup++ - vtmp0 + *vddn++;
    vtmp2 = *(vd+1) - vtmp0 + *(vd+MN);
    *lvd++ = (idx2*vtmp1 + idy2*vtmp2);
    vd++;  
    ud++; 
    for (j=1;j<ny-1;j++) {
      utmp0 = (*ud) * 2.0;
      utmp1 =  *udup++ - utmp0 + *uddn++;
      utmp2 = *(ud+1) - utmp0 + *(ud-1);
      *lud++ = (idx2*utmp1 + idy2*utmp2);
      ud++; 
      vtmp0 = (*vd) * 2.0;
      vtmp1 =  *vdup++ - vtmp0 + *vddn++;
      vtmp2 = *(vd+1) - vtmp0 + *(vd-1);
      *lvd++ = (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++ = (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++ = (idx2*vtmp1 + idy2*vtmp2);
    vd++; 
  }

  /* lastly we must do i=nx for the V laplacian! */

  vtmp0 = (*vd) * 2.0;  // here we do j=0
  vtmp1 =  *vdup++ - vtmp0 + *vddn++;
  vtmp2 = *(vd+1) - vtmp0 + *(vd+MN);
  *lvd++ = (idx2*vtmp1 + idy2*vtmp2);
  vd++;  
  for (j=1;j<ny-1;j++) {
    vtmp0 = (*vd) * 2.0;
    vtmp1 =  *vdup++ - vtmp0 + *vddn++;
    vtmp2 = *(vd+1) - vtmp0 + *(vd-1);
    *lvd++ = (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++; 


  /* finally we have to apply the BC */

  lvd0 = LV->data;
  lvd1 = LV->data+M;
  lvd2 = lvd - M;

  for (j=0;j<ny;j++) {
#ifdef SLIP
    *lvd0++ = *lvd1++;  /* here lap(V)_x = 0 */
    *lvd++ = *lvd2++;
#endif  /* SLIP */
#ifdef NOSLIP
    *lvd0++ = - *lvd1++;  /* here lap(V)=0 */
    *lvd++ = - *lvd2++;
#endif  /* NOSLIP */
  }


}





void newtonian0_u_all(field2D *NU, const field2D *U)
{

  int i,j,im,M,MN;
  double tmp0, tmp1, tmp2;
  double *nud;
  double *ud, *udup, *uddn;

  M = NU->M;
  MN = M-1;

  nud = NU->data;

  ud = &(UU->data[M]);
  udup = &(UU->data[2*M]);
  uddn = UU->data;

  for (i=1;i<nx-1;i++) {
    im = i-1;

    tmp0 = (*ud) * 2.0;  // here we do j=0
    tmp1 =  *udup++ - tmp0 + *uddn++;
    tmp2 = *(ud+1) - tmp0 + *(ud+MN);
    *nud++ = nu_newt*(idx2*tmp1 + idy2*tmp2);
    ud++;  

    for (j=1;j<ny-1;j++) {
      tmp0 = (*ud) * 2.0;
      tmp1 =  *udup++ - tmp0 + *uddn++;
      tmp2 = *(ud+1) - tmp0 + *(ud-1);
      *nud++ = nu_newt*(idx2*tmp1 + idy2*tmp2);
      ud++;  
    }

    tmp0 = (*ud) * 2.0;  // here we do j=ny-1
    tmp1 =  *udup++ - tmp0 + *uddn++;
    tmp2 = *(ud-MN) - tmp0 + *(ud-1);
    *nud++ = nu_newt*(idx2*tmp1 + idy2*tmp2);
    ud++; 
  }
}



void newtonian0_v_all(field2D *NV, const field2D *V)
{

  int i,j,im,M,MN;
  double tmp0, tmp1, tmp2;
  double *nvd;
  double *vd, *vdup, *vddn;
  //  double h;

  M = NV->M;
  MN = M-1;

  nvd = NV->data;

  vd = &(VV->data[M]);
  vdup = &(VV->data[2*M]);
  vddn = VV->data;

  for (i=1;i<nx;i++) {
    im = i-1;

    tmp0 = (*vd) * 2.0;  // here we do j=0
    tmp1 =  *vdup++ - tmp0 + *vddn++;
    tmp2 = *(vd+1) - tmp0 + *(vd+MN);
    *nvd++ = nu_newt*(idx2*tmp1 + idy2*tmp2);
    vd++;  

    for (j=1;j<ny-1;j++) {
      tmp0 = (*vd) * 2.0;
      tmp1 =  *vdup++ - tmp0 + *vddn++;
      tmp2 = *(vd+1) - tmp0 + *(vd-1);
      *nvd++ = nu_newt*(idx2*tmp1 + idy2*tmp2);
      vd++;  
    }

    tmp0 = (*vd) * 2.0;  // here we do j=ny-1
    tmp1 =  *vdup++ - tmp0 + *vddn++;
    tmp2 = *(vd-MN) - tmp0 + *(vd-1);
    *nvd++ = nu_newt*(idx2*tmp1 + idy2*tmp2);
    vd++;  
  }
}


     
/* **** BIHARMONIC LAPLACIAN FUNCTIONS */

/* this takes the field XX and returns with the laplacian LXX
   which is of x dimension two less than XX  */

void bi_laplacian0_all(field2D *LXX, const field2D *XX)
{

  int i,j,im,M,MN,N;
  double utmp0,utmp1,utmp2;
  double *lxd, *lvd, *lvd0, *lvd1, *lvd2;
  double *xd, *xdup, *xddn;

  M = LXX->M;
  N = XX->N;
  MN = M-1;

  lxd = LXX->data;

  xd = &(XX->data[M]);
  xdup = &(XX->data[2*M]);
  xddn = XX->data;

  for (i=1;i<N-1;i++) {
    utmp0 = (*xd) * 2.0;  // here we do j=0
    utmp1 =  *xdup++ - utmp0 + *xddn++;
    utmp2 = *(xd+1) - utmp0 + *(xd+MN);
    *lxd++ = (idx2*utmp1 + idy2*utmp2);
    xd++; 

    for (j=1;j<ny-1;j++) {
      utmp0 = (*xd) * 2.0;
      utmp1 =  *xdup++ - utmp0 + *xddn++;
      utmp2 = *(xd+1) - utmp0 + *(xd-1);
      *lxd++ = (idx2*utmp1 + idy2*utmp2);
      xd++; 
    }

    utmp0 = (*xd) * 2.0;  // here we do j=ny-1
    utmp1 =  *xdup++ - utmp0 + *xddn++;
    utmp2 = *(xd-MN) - utmp0 + *(xd-1);
    *lxd++ = (idx2*utmp1 + idy2*utmp2);
    xd++;  
  }

}


void bi_dyy_all(field2D *LXX, const field2D *XX)
{

  int i,j,im,M,MN,N;
  double utmp0,utmp1,utmp2;
  double *lxd, *lvd, *lvd0, *lvd1, *lvd2;
  double *xd, *xdup, *xddn;

  M = LXX->M;
  N = XX->N;
  MN = M-1;

  lxd = LXX->data;

  xd = &(XX->data[M]);
  xdup = &(XX->data[2*M]);
  xddn = XX->data;

  for (i=1;i<N-1;i++) {
    utmp0 = (*xd) * 2.0;  // here we do j=0
    utmp2 = *(xd+1) - utmp0 + *(xd+MN);
    *lxd++ = idy2*utmp2;
    xd++; 

    for (j=1;j<M-1;j++) {
      utmp0 = (*xd) * 2.0;
      utmp2 = *(xd+1) - utmp0 + *(xd-1);
      *lxd++ = idy2*utmp2;
      xd++; 
    }

    utmp0 = (*xd) * 2.0;  // here we do j=ny-1
    utmp2 = *(xd-MN) - utmp0 + *(xd-1);
    *lxd++ = idy2*utmp2;
    xd++;  
  }

}

void bi_dyy_all_add(field2D *LXX, const field2D *XX)
{

  int i,j,im,M,MN,N;
  double utmp0,utmp1,utmp2;
  double *lxd, *lvd, *lvd0, *lvd1, *lvd2;
  double *xd, *xdup, *xddn;

  M = LXX->M;
  N = XX->N;
  MN = M-1;

  lxd = LXX->data;

  xd = &(XX->data[M]);
  xdup = &(XX->data[2*M]);
  xddn = XX->data;

  for (i=1;i<N-1;i++) {
    utmp0 = (*xd) * 2.0;  // here we do j=0
    utmp2 = *(xd+1) - utmp0 + *(xd+MN);
    *lxd++ += idy2*utmp2;
    xd++; 

    for (j=1;j<M-1;j++) {
      utmp0 = (*xd) * 2.0;
      utmp2 = *(xd+1) - utmp0 + *(xd-1);
      *lxd++ += idy2*utmp2;
      xd++; 
    }

    utmp0 = (*xd) * 2.0;  // here we do j=ny-1
    utmp2 = *(xd-MN) - utmp0 + *(xd-1);
    *lxd++ += idy2*utmp2;
    xd++;  
  }

}


void bi_dxx_all(field2D *LXX, const field2D *XX)
{

  int i,j,im,M,MN,N;
  double utmp0,utmp1,utmp2;
  double *lxd, *lvd, *lvd0, *lvd1, *lvd2;
  double *xd, *xdup, *xddn;

  M = LXX->M;
  N = XX->N;
  MN = M-1;

  lxd = LXX->data;

  xd = &(XX->data[M]);
  xdup = &(XX->data[2*M]);
  xddn = XX->data;

  for (i=1;i<N-1;i++) {
    for (j=0;j<M;j++) {
      utmp0 = (*xd) * 2.0;
      utmp1 =  *xdup++ - utmp0 + *xddn++;
      *lxd++ = idx2*utmp1;
      xd++; 
    }
  }

}


void bi_dxx_all_add(field2D *LXX, const field2D *XX)
{

  int i,j,im,M,MN,N;
  double utmp0,utmp1;
  double *lxd, *lvd, *lvd0, *lvd1, *lvd2;
  double *xd, *xdup, *xddn;

  M = LXX->M;
  N = XX->N;
  MN = M-1;

  lxd = LXX->data;

  xd = &(XX->data[M]);
  xdup = &(XX->data[2*M]);
  xddn = XX->data;