tracer.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.
 *
 * For a copy of the GNU General Public License see LICENSE or go to
 *     http://www.gnu.org/licenses/gpl.html
 */

/* ---  passive.c ----  */

#include "bdefs.h"
#include "tracer.h"
#include "breaking.h"


switch_t tracer_flag;
tracer_source_t source_flag;

field2D *QQ;
field2D *QQT1, *QQT2, *QQT3;
field2D *QDISS, *QADV;
//field2D *LQQ;     // laplacian of QQ what size should it be?  

/* Variables for breaking-wave induced tracer eddy viscosity */

switch_t breaking_tracer_flag;

field2D *BRNU_Q=NULL;
field2D *BRNU_QBX=NULL;
field2D *BRNU_QBY=NULL;

/* for the tracer advection  equation  */

double *uhq_dn;

/* Tracer decay and diffusivities */

double half_life;  // this is the half_life
double lambda;     // this is the decay scale
double kappa;      // this is the diffusivity
//double kappa_bi;   // this would be the biharmonic diffusivity

int i_source;
int j_source;
double source_start_time;
double source_end_time;
double point_source_flux;

void passive_tracer_report(tracerinfo *T)
{

  fprintf(stdout,"------------- * Passive Tracer Report * --------------------------\n");

  /* Then model dimension information */
  /* - Passive Tracer - */

  if (T->tracer_flag == ON) {
    fprintf(stdout,"Passive Tracer: tracer diffusivity kappa=%g (m^2/s)  half_life=%g (s)\n",T->kappa,T->half_life);
    fprintf(stdout,"Passive Tracer: i=%d j=%d src_flux = %g (A/(m^2 s)  start_time=%g (sec), end_time=%g (sec)\n",T->i_source,T->j_source,
            T->point_source_flux, T->source_start_time, T->source_end_time );
  }
  else
    fprintf(stdout,"Passive Tracer:\t module not used in this run\n");

  fprintf(stdout,"------------------------------------------------------------------\n");
}

void init_passive_tracer_variables(int nx,int ny, tracerinfo *T)
{


   tracer_flag = T->tracer_flag;

   if (tracer_flag == OFF)
     return;

   QQ = allocate_field2D_grid(nx+1,ny, ETAGRID);
   QQT1 = allocate_field2D_grid(nx-1,ny, ETAGRID);
   QQT2 = allocate_field2D_grid(nx-1,ny, ETAGRID);
   QQT3 = allocate_field2D_grid(nx-1,ny, ETAGRID);
   //   QSRC = allocate_field2D_grid(nx-1,ny, ETAGRID);
   QDISS = allocate_field2D_grid(nx-1,ny, ETAGRID);
   QADV = allocate_field2D_grid(nx-1,ny, ETAGRID);
   //   LQQ = allocate_field2D(nx+1,ny);

   field2D_set_to_value(QQ,0);
   field2D_set_to_value(QQT1,0);
   field2D_set_to_value(QQT2,0);
   field2D_set_to_value(QQT3,0);
   //   field2D_set_to_value(QSRC,0);
   field2D_set_to_value(QDISS,0);
   field2D_set_to_value(QADV,0);
   //   field2D_set_to_value(LQQ,0);

   DTVARS->QT1 = QQT1;
   DTVARS->QT2 = QQT2;
   DTVARS->QT3 = QQT3;

   uhq_dn = (double * ) g_malloc(ny * sizeof(double));

   /* Here is the source stuff */


   if (T->source_flag != TRACER_POINT_SOURCE)
     funwaveC_error_message("** funwaveC, tracer.c:init_passive_tracer_variables() - source_flag != TRACER_POINT_SOURCE\n");

   /* Check i and j values */

   if ((T->i_source < 1)||(T->i_source > nx)) 
     funwaveC_error_message("** tracer.c: init_passive_tracer_variables() - i is out of bounds\n");



   if ((T->j_source < 0)||(T->j_source > ny-1)) 
     funwaveC_error_message("** tracer.c: init_passive_tracer_variables() - j is out of bounds\n");

   i_source = T->i_source;
   j_source = T->j_source;
   source_start_time = T->source_start_time;
   source_end_time = T->source_end_time;
   point_source_flux = T->point_source_flux;
   

   kappa = T->kappa;
   half_life = T->half_life;

   if (half_life != 0.0)
     lambda = 1.0/half_life;
   else
     lambda = 0.0;


   /* Now deal with the breaking-wave induced eddy viscosity */

   BRNU_Q=breaking_eddy_viscosity();

   if (BRNU_Q == NULL) {
     breaking_tracer_flag = OFF;
   }
   else {
     breaking_tracer_flag = ON;
   }

}


void deallocate_passive_tracer_variables()
{
  if (tracer_flag == ON) {
    deallocate_field2D(QQ);
    deallocate_field2D(QQT1);
    deallocate_field2D(QQT2);
    deallocate_field2D(QQT3);

    g_free(uhq_dn);
  }
}


/* ----------  NUMERICAL MODEL PASSIVE TRACER DYNAMICS ROUTINUES ------------------- */

void  calc_passive_tracer_source_add(field2D *QT)
{

  if ( (global_time > source_start_time)&& (global_time < source_end_time) ) {
    DR2(QT, i_source, j_source ) += point_source_flux;
  }

}

void calc_passive_tracer_dissipation(field2D *QD, const field2D *QQ, const field2D *HBX, const field2D *HBY, double kappa);
void  calc_passive_tracer_advection(field2D *QADV, const field2D *QQ, const field2D *U, const field2D *V, const field2D *HBX, const field2D *HBY);
void  calc_passive_tracer_source_add(field2D *QT);  


/* This function computes the rhs terms of the tracer equation in ETA coordinates
    
    h q_{n+1}(i,j) = hq_n(i,j) + dt*qrhs(i,j)  - etc. 

    This does not worry explicitly about the type of time stepping 

  */


void   calc_passive_tracer_rhs(field2D *QT, const field2D *U, const field2D *V, const field2D *Q, const field2D_depth_var *Depth_Var)
{
  int M = V->M;

  const double *qd = &(Q->data[M]);
  const double *qdissd;
  const double *qadvd; 

  double *qtd = QT->data;
  double *hd = REF2(H, 1);
  double qt;
  
  int i, j, max;
  double qq, hh;

  field2D *HBX = Depth_Var->HBX;
  field2D *HBY = Depth_Var->HBY;
  field2D *H = Depth_Var->H;

  calc_passive_tracer_dissipation(QDISS, QQ, HBX, HBY, kappa);
  calc_passive_tracer_advection(QADV, QQ, U, V, HBX, HBY);

  qadvd = &(QADV->data[0]);
  qdissd = &(QDISS->data[0]);

  max = (nx-1)*ny;

  for (i=1;i<nx;i++) {
    for (j=0;j<ny;j++) {
      hh = *hd++;
      qq = *qd++;
      qt = -lambda*qq*hh  +  (*qdissd++) + (*qadvd++); 
      *qtd++ = qt;
    }
  }

  calc_passive_tracer_source_add(QT);  
}



/* This function calculates the  -  (d/dx(huq) + d/dy(hvq)) terms */

void  calc_passive_tracer_advection(field2D *QADV, const field2D *QQ, const field2D *U, const field2D *V, const field2D *HBX, const field2D *HBY)
{

  int M = V->M;
  int NADV = QADV->N;

  const double *udup = REF2(U, 1);
  const double *vd = REF2(V, 1);
  const double *qd =  REF2(QQ, 1);
  const double *qdup = REF2(QQ, 2);
  const double *hbxd = REF2(HBX,1);
  const double *hbyd = REF2(HBY, 1);

  double *qadvd =  REF2(QADV, 0);
  
  double *uhq_dnd = uhq_dn;
  
  int i,j,jp,jm;
  double hbx, hby;
  double qq, qqm, qqp, vm, uhq_up;
  double vhq_up, vhq_dn;    


  // nonlinear advective terms terms 1) d/dx(huq) 2) d/dy(huq)  
  // Numerically: d_x( u h^x q^x)   and d_y(v h^y q^y)

  /* First zero out the beach boundary x advective flux */

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


  for (i=0;i<NADV-1;i++) {

     // do j = 0 then j=1..ny-2, then j=ny-1

    uhq_dnd = uhq_dn;  /* reset the down uhq pointer */

    for (j=0;j<M;j++) {

      hbx = *hbxd++;
      qq = *qd++;


     // this is the d/dx(hUq) term

      uhq_up = (*udup++) * hbx  * (*qdup++ + qq);
      *qadvd++  = -0.5 * idx * (uhq_up - *uhq_dnd);
      *uhq_dnd++ = uhq_up;
    }
  }

   /* Here at the offshore boundary the advective flux is also zero */

  uhq_dnd = uhq_dn;
  for (j=0;j<M;j++) {
    *qadvd++  = 0.5 * idx * ( *uhq_dnd++);
  }

  /* Done with d/dx(huq) */
  /* Now onto d/dy (hvq) */

  qadvd = REF2(QADV, 0);
  qd = REF2(QQ, 1);  

  /* Now do the y-advective nonlinear term:  d/dy(hVq) */

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

    /* First do j=0 down part */

    hby = *(hbyd+(M-1));

    qqm = *(qd+ (M-1));
    qq = *qd++;
    vm = *(vd + (M-1));

    vhq_dn = -0.5 * idy * hby * vm * (qq + qqm);

    for (j=0;j<M-1;j++) { /* then do rest of j=0 to j=M-2 */
      hby = *hbyd++;
      qqp = *qd++;

      vhq_up = -0.5 * idy * hby * (*vd++) * (qq + qqp);
      *qadvd++ +=  (vhq_up-vhq_dn);
      vhq_dn = vhq_up;
      qq = qqp;

    }

    /* do j=M-1 only up part needed */

    hby = *hbyd++;
    vhq_up = -0.5 * idy * hby * (*vd++) * (qq + *(qd-M));
    *qadvd++ +=  (vhq_up-vhq_dn);
  }

  /* Here done with the d/dy(hvq) term */

}



void calc_passive_tracer_dissipation(field2D *QDISS, const field2D *QQ, const field2D  *HBX, const field2D *HBY, double kappa)
{


  const int M = QQ->M;
  const int NQ = QQ->N;

  static field2D *DQDX=NULL;
  static field2D *DQDY=NULL;
  static field2D *TMPQ=NULL;
  static field2D *BRNU_QBX=NULL;
  static field2D *BRNU_QBY=NULL;

  /* Note no eddy viscosity here = I should use the breaking eddy viscosity I suppose */


  if ( NOT_ALLOCATED(DQDX) ) {   /* what are the sizes of this??  */
    DQDX = allocate_field2D_grid(NQ-1, M, UTGRID);
    DQDY = allocate_field2D_grid(NQ, M, VTGRID);
    TMPQ = allocate_field2D_grid(NQ, M, ETAGRID);
  }

  /* QDISS = \delta_x ( h^x \delta_x Q ) */

  field2D_diffx(DQDX, QQ, idx);            /* DQDX is on an U point */
  field2D_self_multiply(DQDX, HBX);       /* multiply it by the depth (on U point) */

  /* += \delta_y ( h^y \delta_y Q ) */

  field2D_diffy_eta_to_v(DQDY, QQ, idy);          /* DQDY (N+1,M) */  
  field2D_self_multiply(DQDY,HBY);               /* multiply it by the depth (on eta point) */

  if (breaking_tracer_flag == ON) {
    if ( NOT_ALLOCATED(BRNU_QBX) ) {
      BRNU_QBX = allocate_field2D_grid(NQ-1, M, UTGRID);
      BRNU_QBY = allocate_field2D_grid(NQ-2, M, VGRID);
    }

    field2D_barx(BRNU_QBX, BRNU_Q);     /* set up eddy viscosity on x points */
    field2D_add_const(BRNU_QBX, kappa);  /* add the background eddy viscosity */
    field2D_bary_eta_to_v(BRNU_QBY, BRNU_Q);     /* same thing for y points */
    field2D_add_const(BRNU_QBY, kappa);

    field2D_self_multiply(DQDX, BRNU_QBX);
    field2D_self_multiply(DQDY, BRNU_QBY);