updateab.c

来自「波浪数值模拟」· C语言 代码 · 共 290 行

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



/* updateAB.c   Falk Feddersen - handles the updating of the time stepes, ie U^(n+1) += .... */       
     

#include "updateAB.h"


/* ----------------------------------- */


/* 3rd Order Adams-Bashforth Time Stepping routine for U, V, and ETA  but also includes the
   F1 and G1 terms at the appropriate time level which are also added onto the dU/dt and dV/dt equations */


void update_funwaveC_AB3_FG(field2D *U, field2D *V, field2D *ETA, field2D_time_derivative_var *DTVARS, const double dt)
{


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

  const double *ut1d = &(DTVARS->UT1->data[M]);
  const double *ut2d = &(DTVARS->UT2->data[M]);
  const double *ut3d = &(DTVARS->UT3->data[M]);
  const double *vt1d = &(DTVARS->VT1->data[M]);
  const double *vt2d = &(DTVARS->VT2->data[M]);
  const double *vt3d = &(DTVARS->VT3->data[M]);
  const double *etat1d = DTVARS->ETAT1->data;
  const double *etat2d = DTVARS->ETAT2->data;
  const double *etat3d = DTVARS->ETAT3->data;
/*   const double *f1d = DTVARS->F1_1->data; */
/*   const double *f2d = DTVARS->F1_2->data; */
/*   const double *f3d = DTVARS->F1_3->data; */
/*   const double *g1d = DTVARS->G1_1->data; */
/*   const double *g2d = DTVARS->G1_2->data; */
/*   const double *g3d = DTVARS->G1_3->data; */

  double *ud = &(U->data[M]);
  double *vd = &(V->data[M]);
  double *etad = &(ETA->data[0]);

  double *vd0, *vdN;

  int i,j;

  const int numUT = DTVARS->UT1->num - 2*M;
  const int numVT = DTVARS->VT1->num -2*M;
  const int numETAT = DTVARS->ETAT1->num;

  const double hdt = dt/12.0;

  /* loop over all the U  */

  for (i=0;i<numUT;i++) {
      *ud++ += hdt * ( 23.0 *(*ut1d++) - 16.0 * (*ut2d++) + 5.0* (*ut3d++) );
      //      *ud++ += 2.0* (*f1d++) -3.0* (*f2d++) + *f3d++ ;             /* then add on the F1 terms */
  }

  /* loop over all the V  */

  for (i=0;i<numVT;i++) {
      *vd++  += hdt * ( 23.0 * (*vt1d++) - 16.0 * (*vt2d++) + 5.0* (*vt3d++) );
      //      *vd++ += 2.0* (*g1d++) -3.0* (*g2d++) + *g3d++  ;            /* then add on the F1 terms */
  }

  /* loop over all the ETA  */

  for (i=0;i<numETAT;i++) {
      *etad++ += hdt * ( 23.0 * (*etat1d++) - 16.0 * (*etat2d++) + 5.0* (*etat3d++) );
  }


  /* next update the SLIP boundary conditions on V */

  vd0 = &(V->data[0]);
  vdN = &(V->data[(NV-1)*M]);

  for (j=0;j<M;j++) {
    vd0[j] = vd0[j+M];   // shore boundary condition
    vdN[j] = vdN[j-M];   // far slip BC
  }


}

void update_funwaveC_AB3(field2D *U, field2D *V, field2D *ETA, const field2D *UT1, const field2D *UT2, const field2D *UT3,
                         const field2D *VT1, const field2D *VT2, const field2D *VT3, 
                         const field2D *ETAT1, const field2D *ETAT2, const field2D *ETAT3, const double dt)
{

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

  const double *ut1d = &(UT1->data[M]);
  const double *ut2d = &(UT2->data[M]);
  const double *ut3d = &(UT3->data[M]);
  const double *vt1d = &(VT1->data[M]);
  const double *vt2d = &(VT2->data[M]);
  const double *vt3d = &(VT3->data[M]);
  const double *etat1d = ETAT1->data;
  const double *etat2d = ETAT2->data;
  const double *etat3d = ETAT3->data;

  double *ud = &(U->data[M]);
  double *vd = &(V->data[M]);
  double *etad = &(ETA->data[0]);  

  double *vd0, *vdN;

  int i,j;  

  const int numUT = UT1->num - 2*M;
  const int numVT = VT1->num -2*M;
  const int numETAT = ETAT1->num;

  const double hdt = dt/12.0;

  /* loop over all the U  */

  for (i=0;i<numUT;i++) {
      *ud++ += hdt * ( 23.0 *(*ut1d++) - 16.0 * (*ut2d++) + 5.0* (*ut3d++) );
  }

  /* loop over all the V  */

  for (i=0;i<numVT;i++) {
      *vd++ += hdt * ( 23.0 * (*vt1d++) - 16.0 * (*vt2d++) + 5.0* (*vt3d++) );
  }

  /* loop over all the ETA  */

  for (i=0;i<numETAT;i++) {
      *etad++ += hdt * ( 23.0 * (*etat1d++) - 16.0 * (*etat2d++) + 5.0* (*etat3d++) );
  }


  /* next update the SLIP boundary conditions on V */

  vd0 = &(V->data[0]);
  vdN = &(V->data[(NV-1)*M]);

  for (j=0;j<M;j++) {
    vd0[j] = vd0[j+M];   // shore boundary condition
    vdN[j] = vdN[j-M];   // far slip BC
  }


}



void update_funwaveC_AB2(field2D *U, field2D *V, field2D *ETA, field2D_time_derivative_var *DTVARS, const double dt)
{

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

  const double *ut2d = &(DTVARS->UT2->data[M]);
  const double *ut3d = &(DTVARS->UT3->data[M]);
  const double *vt2d = &(DTVARS->VT2->data[M]);
  const double *vt3d = &(DTVARS->VT3->data[M]);
  const double *etat2d = DTVARS->ETAT2->data;
  const double *etat3d = DTVARS->ETAT3->data;

  double *ud = &(U->data[M]);
  double *vd = &(V->data[M]);
  double *etad = &(ETA->data[0]);  

  double *vd0, *vdN;

  int i,j;  

  const int numUT = DTVARS->UT2->num -2*M;
  const int numVT = DTVARS->VT2->num -2*M;
  const int numETAT = DTVARS->ETAT2->num;

  const double hdt = 0.5*dt;

  /* loop over all the U  */

  for (i=0;i<numUT;i++) {
      *ud++ += hdt*(3.0* (*ut2d++) - *ut3d++);
  }

  /* loop over all the V  */

  for (i=0;i<numVT;i++) {
      *vd++ += hdt*(3.0* (*vt2d++) - *vt3d++);
  }

  /* loop over all the ETA  */

  for (i=0;i<numETAT;i++) {
      *etad++ += hdt*(3.0* (*etat2d++) - *etat3d++);
  }


  /* next update the SLIP boundary conditions on V */

  vd0 = &(V->data[0]);
  vdN = &(V->data[(NV-1)*M]);

  for (j=0;j<M;j++) {
    vd0[j] = vd0[j+M];   // shore boundary condition
    vdN[j] = vdN[j-M];   // far slip BC
  }


}



void update_funwaveC_euler(field2D *U, field2D *V, field2D *ETA, const field2D *UT, const field2D *VT, 
                           const field2D *ETAT, const double dt)
{

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

  const double *utd = &(UT->data[M]);
  const double *vtd = &(VT->data[M]);
  const double *etatd = ETAT->data;


  double *ud = &(U->data[M]);
  double *vd = &(V->data[M]);
  double *etad = &(ETA->data[0]);  

  double *vd0, *vdN;

  int i,j;  

  const int numUT = UT->num -2*M;
  const int numVT = VT->num -2*M;
  const int numETAT = ETAT->num;


  /* loop over all the U  */

  for (i=0;i<numUT;i++) {
      *ud++ += dt * (*utd++);
  }

  /* loop over all the V  */

  for (i=0;i<numVT;i++) {
      *vd++ += dt * (*vtd++);
  }

  /* loop over all the ETA  */

  for (i=0;i<numETAT;i++) {
      *etad++ += dt* (*etatd++);
  }


  /* next update the SLIP boundary conditions on V */

  vd0 = &(V->data[0]);
  vdN = &(V->data[(NV-1)*M]);

  for (j=0;j<M;j++) {
    vd0[j] = vd0[j+M];   // shore boundary condition
    vdN[j] = vdN[j-M];   // far slip BC
  }



}