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



/* funwaveC_timestep.c   Falk Feddersen */       
     

#include "bdefs.h"
#include "timing.h"
#include "depth.h"
#include "funwaveC_timestep.h"
#include "bousinessq_dynamics.h"
#include "bousinessq_operations.h"
#include "tracer.h"
#include "floats.h"
#include "sponge.h"
#include "breaking.h"
#include "eta_source.h"
#include "bottom_stress.h"
#include "lateral_mixing.h"
#include "forcing.h"
#include "updateAB.h"
#include "tridiagonal.h"

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


void   preliminary_estimate_dUdt_dVdt(field2D *UUTtmp, field2D *VVTtmp, const field2D_time_derivative_var *D, double dt )
{

  int i;
  const double hdt = 0.5*dt;
  const int numUT = UUTtmp->num;
  const int numVT = VVTtmp->num;

  const double *u1d = REF2(D->U1,0);
  const double *u2d = REF2(D->U2,0);
  const double *u3d = REF2(D->U3,0);
  const double *v1d = REF2(D->V1,0);
  const double *v2d = REF2(D->V2,0);
  const double *v3d = REF2(D->V3,0);

  double *uuttmpd = REF2(UUTtmp,0);
  double *vvttmpd = REF2(VVTtmp,0);

  for (i=0;i<numUT;i++) {
      *uuttmpd++ = hdt * ( 3.0*(*u1d++) -4.0* (*u2d++) + *u3d++ ) ;     
  }

  for (i=0;i<numVT;i++) {
      *vvttmpd++ = hdt * ( 3.0*(*v1d++) -4.0* (*v2d++) + *v3d++ ) ;     
  }


}



void   preliminary_estimate_dUdt_dVdt_AB2(field2D *UUTtmp, field2D *VVTtmp, const field2D_time_derivative_var *D, double dt )
{

  int i;

  const int numUT = UUTtmp->num;
  const int numVT = VVTtmp->num;

  const double *u1d = REF2(D->U1,0);
  const double *u2d = REF2(D->U2,0);
  const double *v1d = REF2(D->V1,0);
  const double *v2d = REF2(D->V2,0);


  double *uuttmpd = REF2(UUTtmp,0);
  double *vvttmpd = REF2(VVTtmp,0);

  for (i=0;i<numUT;i++) {
      *uuttmpd++ = dt * ( *u1d++ - *u2d++);
  }

  for (i=0;i<numVT;i++) {
      *vvttmpd++ = dt * ( *v1d++ - *v2d++);
  }


}


void preliminary_calculations(field2D_bousinessq_var *BVARS, field2D_depth_var *Depth_Vars, field2D *UT, field2D *VT, field2D *ETAT)
{

  b_calc_uxx(BVARS->UXX, BVARS->U);
  b_calc_huxx(BVARS->HUXX, BVARS->U, Depth_Vars->HBX);
  b_calc_uxy(BVARS->UXY, BVARS->U);
  b_calc_huxy(BVARS->HUXY, BVARS->U, Depth_Vars->HBX);

  b_calc_vxy(BVARS->VXY, BVARS->V);
  b_calc_hvxy(BVARS->HVXY, BVARS->V, Depth_Vars->HBY);
  b_calc_vyy(BVARS->VYY, BVARS->V);
  b_calc_hvyy(BVARS->HVYY, BVARS->V, Depth_Vars->HBY);

  field2D_bary_eta_to_v(BVARS->ETABY, BVARS->ETA);
  field2D_barx(BVARS->ETABX, BVARS->ETA);
  field2D_bary_u_to_vort(BVARS->ETABXY, BVARS->ETABX);

  /* Here the total depth (h+\eta) on all the grids are calculated */
 
  field2D_add_sizes(Depth_Vars->DEPTH, BVARS->ETA, Depth_Vars->H);
  field2D_add(Depth_Vars->DEPTHBX, BVARS->ETABX, Depth_Vars->HBX);
  field2D_add_sizes(Depth_Vars->DEPTHBY, BVARS->ETABY, Depth_Vars->HBY);
  field2D_add(Depth_Vars->DEPTHBXY, BVARS->ETABXY, Depth_Vars->HBXY);
 
  BVARS->UT = UT;
  BVARS->VT = VT;
  BVARS->ETAT = ETAT;

  field2D_multiply(BVARS->HU, Depth_Vars->HBX, BVARS->U);
  field2D_multiply(BVARS->HV, Depth_Vars->HBY, BVARS->V);

  field2D_multiply(BVARS->HUT, Depth_Vars->HBX, BVARS->UT);
  field2D_multiply(BVARS->HVT, Depth_Vars->HBY, BVARS->VT);

}




void calculate_UV_rhs(field2D *UT1, field2D *VT1, const field2D_bousinessq_var *BVARS, const field2D_depth_var *Depth_Vars, dynamics_t Dynamics)
{

  int M = UT1->M;
  int NU = UT1->N - 2;
  int NV = VT1->N - 2;

  field2D *U, *V;
  field2D *ETA,  *ETABX, *ETABY;
  field2D *HU, *HV, *HUT, *HVT;
  field2D *UT, *VT;
  field2D *HBX, *HBY, *ZA;

  static field2D *BUNL=NULL, *BVNL=NULL;
  static field2D *PGU=NULL, *PGV=NULL;
  static field2D *F1t=NULL, *G1t=NULL;
  static field2D *Ft=NULL, *Gt=NULL;
  static field2D *F2=NULL, *G2=NULL;
  static field2D *FWK=NULL, *GWK=NULL;
  static field2D *TAUBX=NULL, *TAUBY=NULL;
  static field2D *DISS_U=NULL, *DISS_V=NULL;

  if ( NOT_ALLOCATED( PGU ) ) {
    PGU = allocate_field2D_grid( NU, M, UGRID);
    PGV = allocate_field2D_grid( NV, M, VGRID);
    BUNL = allocate_field2D_grid( NU, M, UGRID);
    BVNL = allocate_field2D_grid( NV, M, VGRID);
    FWK = allocate_field2D_grid( NU, M, UTGRID);
    GWK = allocate_field2D_grid( NV, M, VTGRID);
    //    F2 = allocate_field2D_grid( NU, M, UGRID);
    //    G2 = allocate_field2D_grid( NV, M, VGRID);
    //    Ft = allocate_field2D_grid( NU, M, UGRID);
    //    Gt = allocate_field2D_grid( NV, M, VGRID);
    F1t = allocate_field2D_grid( NU, M, UGRID);
    G1t = allocate_field2D_grid( NV, M, VGRID);
    TAUBX = allocate_field2D_grid( NU, M, UTGRID);
    TAUBY = allocate_field2D_grid( NV, M, VTGRID);
    DISS_U = allocate_field2D_grid( NU, M, UTGRID);
    DISS_V = allocate_field2D_grid( NV, M, VTGRID);
  }

  U = BVARS->U;
  V = BVARS->V;
  ETA = BVARS->ETA;
  ETABX = BVARS->ETABX;
  ETABY = BVARS->ETABY;
  HU = BVARS->HU;
  HV = BVARS->HV;
  UT = BVARS->UT;
  VT = BVARS->VT;
  HUT = BVARS->HUT;
  HVT = BVARS->HVT;

  /* Depth related variables */

  HBX = Depth_Vars->HBX;
  HBY = Depth_Vars->HBY;
  ZA =  Depth_Vars->ZA;

  /* Here do the U & V rhs momentum calculations */

  switch (Dynamics) {
  case WEI_KIRBY_DYNAMICS:
    //    b_calc_uvmom_Ft_Gt(Ft, Gt, ETA, UT, VT, HUT, HVT);            /* calculates the Ft and Gt terms in funwave */
    //    b_calc_uvmom_F2_G2(F2, G2, ETA,  ZA, U, V, HU, HV);           /* calculates the F2 and G2 terms in funwave */
    b_calc_uvmom_FG_WK(FWK, GWK, ETA,  ZA, U, V, HU, HV, UT, VT, HUT, HVT );     /* calculates combined (F2,G2) and (Ft,Gt) terms in funwave */
  case NWOGU_DYNAMICS:

  case PEREGRINE_DYNAMICS:
    b_calc_uvmom_F1t_G1t(F1t, G1t, HBX, HBY, UT, VT, HUT, HVT);       /* calculates the Ft and Gt terms in funwave */
  case NSWE_DYNAMICS:
    b_calc_u_nonlinear(BUNL, U, V);                  /*  u*u_x + v*u_y */
    b_calc_v_nonlinear(BVNL, U, V);                  /*  u*v_x + v*v_y */
  case LINEAR_DYNAMICS:
    b_calc_pressure_gradient(PGU, PGV, BVARS->ETA);  /* -g d(eta)/dx & -g d(eta)/dy */
    rhs_uvmom_bottom_stress(TAUBX, TAUBY, U, V, HBX, HBY, ETABX, ETABY );   // note that here c_d is not defined but inside bottom_stress.c
    rhs_uvmom_lateral_mixing(DISS_U, DISS_V, U, V, Depth_Vars);                       /* this is in lateral_friction.c  (and .h) */
  }



  /* Add up all the dU/dt and dV/dt terms */

  field2D_copy(UT1, TAUBX);
  field2D_copy(VT1, TAUBY);     /* VT has N+2 dimensions to TAUBY */

  switch (Dynamics) {

  case WEI_KIRBY_DYNAMICS:
    field2D_self_add_three(UT1, FWK, F1t, BUNL);
    field2D_self_add_three(VT1, GWK, G1t, BVNL);
    break;
  case NWOGU_DYNAMICS:

  case PEREGRINE_DYNAMICS:
    field2D_self_add_two(UT1, F1t, BUNL);
    field2D_self_add_two(VT1, G1t, BVNL);
    break;
  case NSWE_DYNAMICS:
    field2D_self_add(UT1, BUNL);
    field2D_self_add(VT1, BVNL);
    break;


  }

  field2D_self_add_two(UT1, PGU, DISS_U);
  field2D_self_add_two(VT1, PGV, DISS_V);

  /* Now do the external forcing */

  rhs_uvmom_forcing_add(UT1, VT1);     /* add on the external alongshore forcing */

  /* Now do the sponge layer */

  rhs_uvmom_sponge_layer_add(UT1, VT1, U, V);

  /* Now do the wave breaking stuff  */

  rhs_uvmom_breaking_add(UT1, VT1, U, V, BVARS->ETAT);

  /* just to be safe apply boundary conditions on UT1 and VT1 */

  field2D_apply_zero_bc(UT1);
  field2D_apply_nogradient_bc(VT1);

}







void calculate_eta_rhs(field2D *ETAT, const field2D_bousinessq_var *BVARS, const field2D_depth_var *Depth_Vars, dynamics_t Dynamics)
{

  int NETA = ETAT->N;
  int M = ETAT->M;

  field2D *U, *V, *ETABX, *ETABY;
  field2D *UXX, *UXY, *VXY, *VYY;
  field2D *HUXX, *HUXY, *HVXY, *HVYY;

  field2D *HBX, *HBY;


  U = BVARS->U;
  V = BVARS->V;
  ETABX = BVARS->ETABX;
  ETABY = BVARS->ETABY;
  UXX = BVARS->UXX;
  UXY = BVARS->UXY;
  VXY = BVARS->VXY;
  VYY = BVARS->VYY;
  HUXX = BVARS->HUXX;
  HUXY = BVARS->HUXY;
  HVXY = BVARS->HVXY;
  HVYY = BVARS->HVYY;

  HBX = Depth_Vars->HBX;
  HBY = Depth_Vars->HBY;

  switch (Dynamics) {
  case LINEAR_DYNAMICS:  /* use linearized (mean) depth */
    b_calc_eta_div_hu(ETAT, U, V, HBX, HBY);
    break;
  case NSWE_DYNAMICS:  /* use true depth (h+\eta) */
  case PEREGRINE_DYNAMICS:
    b_calc_eta_div_hu(ETAT, U, V, Depth_Vars->DEPTHBX, Depth_Vars->DEPTHBY); 
    break;
  case NWOGU_DYNAMICS:
    b_calc_eta_E(ETAT, U, V, Depth_Vars->DEPTHBX, Depth_Vars->DEPTHBY, HBX, HBY, UXX, VXY, HUXX, HVXY, UXY, VYY, HUXY, HVYY);
    break;
  case WEI_KIRBY_DYNAMICS:
    b_calc_eta_E(ETAT, U, V, Depth_Vars->DEPTHBX, Depth_Vars->DEPTHBY, HBX, HBY, UXX, VXY, HUXX, HVXY, UXY, VYY, HUXY, HVYY);
    b_calc_eta_E2(ETAT, ETABX, ETABY, HBX, HBY, UXX, VXY, HUXX, HVXY, UXY, VYY, HUXY, HVYY);   /* d(eta)/dt += E2 */
  }

  rhs_calc_eta_source_f_add(ETAT);    /*  d(eta)/dt += f(x,y,t);  this is the source func for eta in eta_source.c */



}


void cycle_time_derivative_var(field2D_time_derivative_var *D)
{
  field2D *TMP;

  /* First do the dU/dt pointers */

  TMP = D->UT3;
  D->UT3 = D->UT2;
  D->UT2 = D->UT1;
  D->UT1 = TMP;

  /* Then  the dV/dt pointers */

  TMP = D->VT3;
  D->VT3 = D->VT2;
  D->VT2 = D->VT1;
  D->VT1 = TMP;

  /* Then the d(eta)/dt pointers */

  TMP = D->ETAT3;
  D->ETAT3 = D->ETAT2;
  D->ETAT2 = D->ETAT1;
  D->ETAT1 = TMP;

  /* Then the d(tracer)/dt pointers */

  TMP = D->QT3;
  D->QT3 = D->QT2;
  D->QT2 = D->QT1;
  D->QT1 = TMP;

/*   /\* Then the F1 pointers *\/ */

/*   TMP = D->F1_3; */
/*   D->F1_3 = D->F1_2; */
/*   D->F1_2 = D->F1_1; */
/*   D->F1_1 = TMP; */

/*   /\* Then the G1 pointers *\/ */

/*   TMP = D->G1_3; */
/*   D->G1_3 = D->G1_2; */
/*   D->G1_2 = D->G1_1; */
/*   D->G1_1 = TMP; */

  /* Now do the U & V pointers - this is more tricky */

  TMP = D->U3;
  D->U3 = D->U2;
  D->U2 = D->U1;
  D->U1 = TMP;
  field2D_copy(D->U1, UU);

  /* Now do the U & V pointers - this is more tricky */

  TMP = D->V3;
  D->V3 = D->V2;
  D->V2 = D->V1;
  D->V1 = TMP;
  field2D_copy(D->V1, VV);