funwavec_timestep.c

波浪数值模拟

}



field2D *UUT1, *UUT2, *UUT3;
field2D *VVT1, *VVT2,  *VVT3;
field2D *ETAT1, *ETAT2, *ETAT3;

field2D *UUtmp, *VVtmp;
field2D *UUTtmp, *VVTtmp;  /* preliminary dU/dt and dV/dt estimates for use in Wei & Kirby calculations.  Estimate FIRST! */

field2D_bousinessq_var *BVARS;
field2D_time_derivative_var *DTVARS;  /* collection of all the time derivative variables */

dynamics_t Bousinessq_Dynamics;

void funwaveC_timestep_AB3(timing *T)
{

   double dt;

   dt = T->dt;

   preliminary_estimate_dUdt_dVdt(UUTtmp, VVTtmp, DTVARS, dt );                  /* Wei&Kirby need dU/dt at this time step for nonlinear terms Ft & Gt */
   preliminary_calculations(BVARS, Depth_Vars, UUTtmp, VVTtmp, DTVARS->ETAT2 );  /* here calc uxx, huxx, for rhs eta & u,v eq etc */

   calculate_eta_rhs( DTVARS->ETAT1, BVARS, Depth_Vars, Bousinessq_Dynamics );                         /* Calculate  d(eta)/dt = rhs*/
   calculate_UV_rhs( DTVARS->UT1, DTVARS->VT1, BVARS, Depth_Vars, Bousinessq_Dynamics );               /* Calculate d([u,v])/dt = rhs  */
   //   calculate_uvmom_rhs_FG_terms( DTVARS->F1_1, DTVARS->G1_1, BVARS, Depth_Vars ); /* calculate F1 and G1 with U & V (not dU/dt and dV/dt) */

   if (tracer_flag == ON)
     calc_passive_tracer_rhs( DTVARS->QT1, UU, VV, QQ, Depth_Vars );

   update_funwaveC_AB3_FG(UUtmp,VVtmp, ETA, DTVARS, dt);
   tridiagonal_solve_UV(UU, VV, UUtmp, VVtmp);


   if(tracer_flag == ON)
     update_passive_tracerAB3(QQ, DTVARS, Depth_Vars,dt);

   if(floats_flag == ON)
     update_floats_AB2(UU, VV, dt) ;

   cycle_time_derivative_var(DTVARS);   /* this makes UT3 = UT2, UT2 = UT1, UT1 = UT3, etc in DTVARS */

   T->time_index++;
   T->current_time += dt;
}


void funwaveC_timestep_AB2(timing *T)
{

   double dt;

   dt = T->dt;

   preliminary_estimate_dUdt_dVdt_AB2(UUTtmp, VVTtmp, DTVARS,  dt );            /* Wei&Kirby need dU/dt at this time step for nonlinear terms Ft & Gt */
   preliminary_calculations(BVARS, Depth_Vars, UUTtmp, VVTtmp, DTVARS->ETAT3 );        /* here calculate uxx, huxx, for use in rhs eta & (u,v) eq etc */

   calculate_eta_rhs( DTVARS->ETAT2, BVARS, Depth_Vars ,Bousinessq_Dynamics);            /* Calculate  d(eta)/dt = rhs*/
   calculate_UV_rhs( DTVARS->UT2, DTVARS->VT2, BVARS, Depth_Vars, Bousinessq_Dynamics );            /* Calculate d([u,v])/dt = rhs  */
   //   calculate_uvmom_rhs_FG_terms( DTVARS->F1_2, DTVARS->G1_2, BVARS, Depth_Vars ); /* calculate F1 and G1 with U & V (not dU/dt and dV/dt) */

   if (tracer_flag == ON) 
     calc_passive_tracer_rhs(DTVARS->QT2, UU, VV, QQ, Depth_Vars );


   update_funwaveC_AB2(UUtmp,VVtmp, ETA, DTVARS, dt);
   tridiagonal_solve_UV(UU, VV, UUtmp, VVtmp);

   if(tracer_flag == ON)
     update_passive_tracerAB2(QQ, DTVARS->QT2, DTVARS->QT3, Depth_Vars, dt);

   if(floats_flag == ON)
     update_floats_AB2(UU, VV, dt) ;


   T->time_index += 1;
   T->current_time += dt;
}


/*
   This function is called to do the initial time step only
   We also set UT,... to the proper quantity
   
*/


void funwaveC_timestep_euler(timing *T)
{

   double dt;

   dt = T->dt;

   preliminary_calculations(BVARS, Depth_Vars, UUTtmp, VVTtmp, DTVARS->ETAT1 );    /* here calculate uxx, huxx, for use in rhs eta & (u,v) eq etc */

   calculate_eta_rhs( DTVARS->ETAT3, BVARS, Depth_Vars, Bousinessq_Dynamics );            /* Calculate  */

   calculate_UV_rhs( DTVARS->UT3, DTVARS->VT3, BVARS, Depth_Vars, Bousinessq_Dynamics );            /* Calculate  */

   if (tracer_flag == ON) 
     calc_passive_tracer_rhs(DTVARS->QT3, UU, VV, QQ, Depth_Vars);

   update_funwaveC_euler(UUtmp, VVtmp, ETA, DTVARS->UT3, DTVARS->VT3, DTVARS->ETAT3,  dt); /* start with UT2 not UT1  */
   tridiagonal_solve_UV(UU, VV, UUtmp, VVtmp);

   if (tracer_flag == ON)
     update_passive_tracer_euler(QQ, DTVARS->QT3, Depth_Vars, dt);

   if(floats_flag == ON)
     update_floats_AB2(UU, VV, dt);


   /* update the time stuff */

   T->time_index=1;
   T->current_time += dt;
}






void  init_funwaveC_timestep(int nx, int ny, dynamics_t Dynamics)
{

   b_init_constants(Dynamics);           /* initialze the bousinessq constants */
   Bousinessq_Dynamics = Dynamics;

   UUT1 = allocate_field2D_grid(nx,ny, UGRID);
   UUT2 = allocate_field2D_grid(nx,ny, UGRID);
   UUT3 = allocate_field2D_grid(nx,ny, UGRID);

   VVT1 = allocate_field2D_grid(nx+1,ny, VGRID);
   VVT2 = allocate_field2D_grid(nx+1,ny, VGRID);
   VVT3 = allocate_field2D_grid(nx+1,ny, VGRID);

   ETAT1 = allocate_field2D_grid(nx-1,ny, ETAGRID);
   ETAT2 = allocate_field2D_grid(nx-1,ny, ETAGRID);
   ETAT3 = allocate_field2D_grid(nx-1,ny, ETAGRID);

   UUT2 = allocate_field2D_grid(nx,ny, UGRID);
   UUT3 = allocate_field2D_grid(nx,ny, UGRID);

   field2D_set_to_value(UUT1,0.0);
   field2D_set_to_value(UUT2,0.0);
   field2D_set_to_value(UUT3,0.0);

   field2D_set_to_value(VVT1,0.0);
   field2D_set_to_value(VVT2,0.0); 
   field2D_set_to_value(VVT3,0.0); 

   field2D_set_to_value(ETAT1,0.0);
   field2D_set_to_value(ETAT2,0.0); 
   field2D_set_to_value(ETAT3,0.0); 


   /* Init all the other dynamic variables */

   BVARS = (field2D_bousinessq_var * ) g_malloc(sizeof(field2D_bousinessq_var));

   BVARS->U = UU;   /* collect the global variables */
   BVARS->V = VV;
   BVARS->ETA = ETA;  
   
   BVARS->ETABX = allocate_field2D_grid(nx,ny, UGRID);  // ?? dimensions
   BVARS->ETABY = allocate_field2D_grid(nx-1,ny, VGRID);  // ?? dimensions
   BVARS->ETABXY = allocate_field2D_grid(nx,ny, VORTGRID);  // ?? dimensions

   BVARS->UXX = allocate_field2D_grid(nx,ny, UGRID);  // this is on a U grid w/ bc Uxx=0 at boundaries
   BVARS->UXY = allocate_field2D_grid(nx-1,ny, VGRID);  // ?? dimensions 
   BVARS->VYY = allocate_field2D_grid(nx-1,ny, VGRID);  // ?? dimensions 
   BVARS->VXY = allocate_field2D_grid(nx,ny, UGRID);  // ?? dimensions 

   BVARS->HUXX = allocate_field2D_grid(nx,ny, UGRID);  // this is on a U grid w/ bc Uxx=0 at boundaries
   BVARS->HUXY = allocate_field2D_grid(nx-1,ny, VGRID);  // ?? dimensions 
   BVARS->HVYY = allocate_field2D_grid(nx-1,ny, VGRID);  // ?? dimensions 
   BVARS->HVXY = allocate_field2D_grid(nx,ny, UGRID);  // ?? dimensions 

   BVARS->HU = allocate_field2D_grid(nx,ny, UGRID);  // this is on a U grid w/ bc Uxx=0 at boundaries
   BVARS->HV = allocate_field2D_grid(nx+1,ny, VGRID);  // ?? dimensions 

   BVARS->HUT = allocate_field2D_grid(nx,ny, UGRID);  // this is on a U grid w/ bc Uxx=0 at boundaries
   BVARS->HVT = allocate_field2D_grid(nx+1,ny, VGRID);  // ?? dimensions 

   /* allocate the UUtmp = h (b1 h u_xx + b2 (hu)_xx) variables */

   UUtmp = allocate_field2D_grid(nx,ny, UGRID);
   VVtmp = allocate_field2D_grid(nx+1,ny, VGRID);

   /* set up the initial condition for UUtmp and VVtmp */

   b_calc_UVtmp_from_UV(UUtmp, VVtmp, Depth_Vars->HBX, Depth_Vars->HBY, UU, BVARS->UXX, BVARS->HUXX, VV, BVARS->VYY, BVARS->HVYY);

   UUTtmp = allocate_field2D_grid(nx,ny, UGRID);
   VVTtmp = allocate_field2D_grid(nx+1,ny, VGRID);
   field2D_set_to_value(UUTtmp,0.0);
   field2D_set_to_value(VVTtmp,0.0);



   /* Now set up the time derivative Variables and the F1 and G1 variables */

   DTVARS = (field2D_time_derivative_var * ) g_malloc(sizeof(field2D_time_derivative_var));

   DTVARS->UT1 = UUT1;
   DTVARS->UT2 = UUT2;
   DTVARS->UT3 = UUT3;

   DTVARS->VT1 = VVT1;
   DTVARS->VT2 = VVT2;
   DTVARS->VT3 = VVT3;

   DTVARS->ETAT1 = ETAT1;
   DTVARS->ETAT2 = ETAT2;
   DTVARS->ETAT3 = ETAT3;

   DTVARS->U1 = allocate_field2D_grid(nx,ny,UGRID);
   DTVARS->U2 = allocate_field2D_grid(nx,ny,UGRID);
   DTVARS->U3 = allocate_field2D_grid(nx,ny,UGRID);
   DTVARS->V1 = allocate_field2D_grid(nx+1,ny,VGRID);
   DTVARS->V2 = allocate_field2D_grid(nx+1,ny,VGRID);
   DTVARS->V3 = allocate_field2D_grid(nx+1,ny,VGRID);

   field2D_set_to_value(DTVARS->U1,0.0);
   field2D_set_to_value(DTVARS->U2,0.0);
   field2D_set_to_value(DTVARS->U3,0.0);
   field2D_set_to_value(DTVARS->V1,0.0);
   field2D_set_to_value(DTVARS->V2,0.0);
   field2D_set_to_value(DTVARS->V3,0.0);

/*    DTVARS->F1_1 = allocate_field2D_grid(nx-2,ny,UTGRID); */
/*    DTVARS->F1_2 = allocate_field2D_grid(nx-2,ny,UTGRID); */
/*    DTVARS->F1_3 = allocate_field2D_grid(nx-2,ny,UTGRID); */

/*    DTVARS->G1_1 = allocate_field2D_grid(nx-1,ny,VTGRID); */
/*    DTVARS->G1_2 = allocate_field2D_grid(nx-1,ny,VTGRID); */
/*    DTVARS->G1_3 = allocate_field2D_grid(nx-1,ny,VTGRID); */
 
/*    field2D_set_to_value(DTVARS->F1_1,0.0); */
/*    field2D_set_to_value(DTVARS->F1_2,0.0); */
/*    field2D_set_to_value(DTVARS->F1_3,0.0); */

/*    field2D_set_to_value(DTVARS->G1_1,0.0); */
/*    field2D_set_to_value(DTVARS->G1_2,0.0); */
/*    field2D_set_to_value(DTVARS->G1_3,0.0); */

   /* initialize the tridiagonal solver */

   init_tridiagonal_solver(Depth_Vars->HBX, Depth_Vars->HBY, Dynamics);


}



void deallocate_funwaveC_timestep()
{

  deallocate_field2D(UUT1);
  deallocate_field2D(UUT2);
  deallocate_field2D(UUT3);
  deallocate_field2D(VVT1);
  deallocate_field2D(VVT2); 
  deallocate_field2D(VVT3); 
  deallocate_field2D(ETAT1);
  deallocate_field2D(ETAT2); 
  deallocate_field2D(ETAT3); 

  deallocate_tridiagonal_solver();


}



/* void calculate_UV_rhs(field2D *UT1, field2D *VT1, const field2D_bousinessq_var *BVARS, const field2D_depth_var *Depth_Vars)  // need more here */
/* { */

/*   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 *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); */
/*     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 *\/ */

/*   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 *\/ */
/*   b_calc_pressure_gradient(PGU, PGV, BVARS->ETA);  /\* -g d(eta)/dx & -g d(eta)/dy *\/ */
/*   b_calc_uvmom_F1t_G1t(F1t, G1t, HBX, HBY, UT, VT, HUT, HVT);       /\* calculates the Ft and Gt terms in funwave *\/ */
/*   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 *\/ */
/*   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 terms *\/ */

/*   field2D_copy(UT1, TAUBX); */
/*   field2D_self_add(UT1, F2); */
/*   field2D_self_add(UT1, Ft); */
/*   field2D_self_add(UT1, F1t); */
/*   field2D_self_add(UT1, PGU); */
/* #ifndef LINEAR_DYNAMICS   /\* only add nonlinear terms if the dynamcis arent' linear *\/ */
/*   field2D_self_add(UT1, BUNL); */
/* #endif /\* LINEAR_DYNAMICS *\/ */
/*   field2D_self_add(UT1, DISS_U); */

/*   /\* Add up all the dV/dt terms *\/ */

/*   field2D_copy(VT1, TAUBY);     /\* VT has N+2 dimensions to TAUBY *\/ */
/*   field2D_self_add(VT1, G2); */
/*   field2D_self_add(VT1, G1t); */
/*   field2D_self_add(VT1, Gt); */
/*   field2D_self_add(VT1, PGV); */
/* #ifndef LINEAR_DYNAMICS      /\* only add nonlinear terms if the dynamcis arent' linear *\/ */
/*   field2D_self_add(VT1, BVNL); */
/* #endif /\* LINEAR_DYNAMICS *\/ */
/*   field2D_self_add(VT1, 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); */

/* } */