bousinessq_operations.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_operations.c  */

#include "bdefs.h"



/* Calculates field2D operators such as  U_xx, V_yy etc.  */

/* Returns d^2U/dx^2 on a U grid - applies the boundary condition that U_xx = 0 at boundaries */

void b_calc_uxx(field2D *UXX, const field2D *U)
{
  int i,j;
  const int N = U->N;
  const int M = UXX->M;

  const double *ud = &(U->data[M]);
  double *uxxd = &(UXX->data[M]);

  double *uxxd0 = &(UXX->data[0]);
  double *uxxdN = &(UXX->data[(N-1)*M]);

  for (i=1;i<N-1;i++) {
    for (j=0;j<M;j++)  {
      *uxxd++ = idx2 * ( *(ud+M) - 2.0 * (*ud) + *(ud-M) );
      ud++;
    }
  }

  /* This is the application fo the boundary condition - but may be redundant if Uxx = 0 originally */

  for (j=0;j<M;j++) {
    *uxxd0++ = 0.0;
    *uxxdN++ = 0.0;
  }
}


/* Returns d^2(hU)/dy^2 on a U grid */

void b_calc_huxx(field2D *HUXX, const field2D *U, const field2D *HBX)
{
  int i,j;
  const int N = U->N;
  const int M = HUXX->M;

  const double *ud = &(U->data[M]);
  const double *hbxd = &(HBX->data[M]);

  double *huxxd = &(HUXX->data[M]);
  double huup, hu, hudn;

  double *huxxd0 = &(HUXX->data[0]);
  double *huxxdN = &(HUXX->data[(N-1)*M]);

  for (i=1;i<N-1;i++) {
    for (j=0;j<M;j++)  {
      huup = *(ud+M) * (*(hbxd+M));
      hu = *ud * (*hbxd);
      hudn = *(ud-M) * (*(hbxd-M));

      *huxxd++ = idx2 * ( huup - 2.0*hu + hudn) ;
      ud++;
      hbxd++;
    }
  }

  /* This is the application fo the boundary condition - but may be redundant if (hU)xx = 0 originally */

  for (j=0;j<M;j++) {
    *huxxd0++ = 0.0;
    *huxxdN++ = 0.0;
  }

}




/* Returns d^2(V)/dxdy on a U grid */
/* Note that V = (NX+1,NY),  VXY = (NX,NY)  - U points! */

void b_calc_vxy(field2D *VXY, const field2D *V)
{
  int i,j;
  const int M = VXY->M;
  const int MN = M-1;
  const int NVXY = VXY->N;

  const double *vd = &(V->data[0]);
  double *vxyd = &(VXY->data[0]);
  double tmp1, tmp0;

  const double idxdy = idx*idy;

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

    /* j = 0 */

      tmp1 = *(vd+M) - *(vd+M+MN);
      tmp0 = *vd  - *(vd+MN);;
      *vxyd++ = idxdy * (tmp1 - tmp0);
      vd++;

    for (j=1;j<M;j++)  {
      tmp1 = *(vd+M) - *(vd+M-1);
      tmp0 = *vd  - *(vd-1);;
      *vxyd++ = idxdy * (tmp1 - tmp0);
      vd++;
    }
  }
}


/* Returns d^2(hV)/dxdy on a U grid */
/* Note that V = (NX+1,NY),  HVXY = (NX,NY)  - U points! */

void b_calc_hvxy(field2D *HVXY, const field2D *V, const field2D *HBY)
{
  int i,j;

  const int M = HVXY->M;
  const int MN = M-1;
  const int NHVXY = HVXY->N;

  const double *vd = &(V->data[0]);
  const double *hbyd = &(HBY->data[0]);

  double *hvxyd = &(HVXY->data[0]);
  double tmp1, tmp0;

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

    tmp1 = (*(hbyd+M) * (*(vd+M)))  - (*(hbyd+M+MN) * (*(vd+M+MN)));
    tmp0 = (*hbyd * (*vd))  - (*(hbyd+MN)* (*(vd+MN)));
    *hvxyd++ = idx * idy * (tmp1 - tmp0);
    vd++;
    hbyd++;

    for (j=1;j<M;j++)  {
      tmp1 = (*(hbyd+M) *(*(vd+M)))  - (*(hbyd+M-1) * (*(vd+M-1)));
      tmp0 = *hbyd * (*vd)  - *(hbyd-1)* (*(vd-1));
      //      tmp1 = *(vd+M) - *(vd+M-1);
      //      tmp0 = *vd  - *(vd-1);;
      *hvxyd++ = idx * idy * (tmp1 - tmp0);
      vd++;
      hbyd++;
    }
  }
}


/* ----- HERE ARE THE V POINT ROUTINES -------- */

/* Returns d^2V/dy^2 on a V grid - note that VYY is (nx-1,ny) not (nx+1,ny) */

void b_calc_vyy(field2D *VYY, const field2D *V)
{
  int i,j;
  const int M = VYY->M;
  const int MN = M-1;
  const int NVYY = VYY->N;

  const double *vd = &(V->data[M]);
  double *vyyd = &(VYY->data[0]);

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

    /* j=0 */

    *vyyd++ = idy2 * ( *(vd+1) - 2.0 * (*vd) + *(vd+MN) );
    vd++;

    for (j=1;j<M-1;j++)  {
      *vyyd++ = idy2 * ( *(vd+1) - 2.0 * (*vd) + *(vd-1) );
      vd++;
    }

    /* j=M-1 */

    *vyyd++ = idy2 * ( *(vd-MN) - 2.0 * (*vd) + *(vd-1) );
    vd++;

  }
}


/* Returns d^2(hV)/dx^2 on a V grid */

void b_calc_hvyy(field2D *HVYY, const field2D *V, const field2D *HBY)
{
  int i,j;
  int M = HVYY->M;
  int MN = M-1;
  const int NHVYY = HVYY->N;

  const double *vd = &(V->data[M]);
  const double *hbyd = &(HBY->data[M]);

  double *hvyyd = &(HVYY->data[0]);
  double hvup, hv, hvdn;

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


    /* j = 0 */

    hvup = *(vd+1) * (*(hbyd+1));
    hv = *vd * (*hbyd);
    hvdn = *(vd+MN) * (*(hbyd+MN));

    *hvyyd++ = idx2 * ( hvup - 2.0*hv + hvdn) ;
    vd++;
    hbyd++;

    for (j=1;j<M-1;j++)  {
      hvup = *(vd+1) * (*(hbyd+1));
      hv = *vd * (*hbyd);
      hvdn = *(vd-1) * (*(hbyd-1));

      *hvyyd++ = idx2 * ( hvup - 2.0*hv + hvdn) ;
      vd++;
      hbyd++;
    }

    /* j = M-1 */

    hvup = *(vd-MN) * (*(hbyd-MN));
    hv = *vd * (*hbyd);
    hvdn = *(vd-1) * (*(hbyd-1));

    *hvyyd++ = idx2 * ( hvup - 2.0*hv + hvdn) ;
    vd++;
    hbyd++;

  }

}


/* Returns d^2(U)/dxdy on a V grid */
/* Note that U = (NX,NY),  UXY = (NX-1,NY)  - V points! */

void b_calc_uxy(field2D *UXY, const field2D *U)
{
  int i,j;
  const int M = UXY->M;
  const int MN = M-1;
  const int NUXY = UXY->N;

  const double *ud = &(U->data[M]);
  double *uxyd = &(UXY->data[0]);
  double tmp1, tmp0;

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

    for (j=0;j<M-1;j++)  {
      tmp1 = *(ud+1) - *ud;
      tmp0 = *(ud-M+1)  - *(ud-M);
      *uxyd++ = idx * idy * (tmp1 - tmp0);
      ud++;
    }

    /* j = M-1 */

    tmp1 = *(ud-MN) - *ud;
    tmp0 = *(ud-M-MN)  - *(ud-M);
    *uxyd++ = idx * idy * (tmp1 - tmp0);
    ud++;

  }
}

/* Returns d^2(U)/dxdy on a V grid */
/* Note that U = (NX,NY),  UXY = (NX-1,NY)  - V points! */

void b_calc_huxy(field2D *HUXY, const field2D *U, const field2D *HBX)
{
  int i,j;
  const int M = HUXY->M;
  const int MN = M-1;
  const int NHUXY = HUXY->N;

  const double *ud = &(U->data[M]);
  const double *hbxd = &(HBX->data[M]);

  double *huxyd = &(HUXY->data[0]);
  double tmp1, tmp0;

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

    for (j=0;j<M-1;j++)  {
      tmp1 = (*(hbxd+1) * (*(ud+1))) - ( *hbxd * (*ud));
      tmp0 = (*(hbxd-M+1) * (*(ud-M+1)))   -  ( *(hbxd-M) * (*(ud-M)) );
      *huxyd++ = idx * idy * (tmp1 - tmp0);
      ud++;
      hbxd++;
    }

    /* j = M-1 */

    tmp1 = (*(hbxd-MN) * (*(ud-MN))) - (*hbxd * (*ud));
    tmp0 = (*(hbxd-M-MN) * (*(ud-M-MN)))  - (*(hbxd-M) * (*(ud-M)));
    *huxyd++ = idx * idy * (tmp1 - tmp0);
    ud++;
    hbxd++;


  }
}


/* Calculates  UXPVY = delta_x(U) + \delta_y(V) - Can also be used for  \delta_x(hU) etc. 
    -- U is on a U grid, V is on a VGRID,  UXPVY is on an ETA GRID, ie size (NX-1,NY)            */

void b_calc_ux_plus_vy(field2D *UXPVY, const field2D *U, const field2D *V, double idx, double idy)
{

  int NU = U->N;
  int M = U->M;
  int NV = V->N;

  int i,j;

  double *uxd;
  double *vyd;
  double *uxpvyd = UXPVY->data;

  static field2D *UX=NULL;
  static field2D *VY=NULL;

  if ( NOT_ALLOCATED(UX) ) {
     UX = allocate_field2D_grid( NU-1, M, ETAGRID);
  }

  if ( NOT_ALLOCATED(VY) ) {
     VY = allocate_field2D_grid( NV, M, ETAGRID);
  }

  field2D_diffx(UX, U, idx);
  field2D_diffy_v_to_eta(VY, V, idy);

  uxd = UX->data;
  vyd = &(VY->data[M]);

    for (i=0;i<NU-1;i++)
      for (j=0;j<M;j++) 
        *uxpvyd++ = *uxd++ + *vyd++;

    //  field2D_add(UXPVY, UX, VY); /* This will break because UX has NX-1  and VY has NX+1 !! */

}