shallowwater2d.h

Go to the documentation of this file.

#include <basic.h>
#include <math.h>
#include <matops.h>

#define _SHALLOW_WATER_2D_ "shallow-water-2d"

/* define ndims and nvars for this model */
#undef _MODEL_NDIMS_
#undef _MODEL_NVARS_

#define _MODEL_NDIMS_ 2

#define _MODEL_NVARS_ 3

/* choices for upwinding schemes */
#define _ROE_     "roe"

#define _LLF_     "llf-char"

/* grid direction */
#undef _XDIR_
#undef _YDIR_

#define _XDIR_ 0
#define _YDIR_ 1

#define _ShallowWater2DGetFlowVar_(u,h,uvel,vvel) \
  { \
    h    = u[0]; \
    uvel = u[1] / h; \
    vvel = u[2] / h; \
  }

#define _ShallowWater2DSetFlux_(f,h,uvel,vvel,g,dir) \
  { \
    if (dir == _XDIR_) { \
      f[0] = (h) * (uvel); \
      f[1] = (h) * (uvel) * (uvel) + 0.5 * (g) * (h) * (h); \
      f[2] = (h) * (uvel) * (vvel); \
    } else if (dir == _YDIR_) { \
      f[0] = (h) * (vvel); \
      f[1] = (h) * (uvel) * (vvel); \
      f[2] = (h) * (vvel) * (vvel) + 0.5 * (g) * (h) * (h); \
    } \
  }

#define _ShallowWater2DRoeAverage_(uavg,uL,uR,p) \
  { \
    double h , uvel , vvel ; \
    double hL, uvelL, vvelL; \
    double hR, uvelR, vvelR; \
    _ShallowWater2DGetFlowVar_(uL,hL,uvelL,vvelL); \
    _ShallowWater2DGetFlowVar_(uR,hR,uvelR,vvelR); \
    h    = 0.5 * (hL    + hR   ); \
    uvel = (sqrt(hL)*uvelL + sqrt(hR)*uvelR) / (sqrt(hL) + sqrt(hR)); \
    vvel = (sqrt(hL)*vvelL + sqrt(hR)*vvelR) / (sqrt(hL) + sqrt(hR)); \
    uavg[0] = h; \
    uavg[1] = h*uvel; \
    uavg[2] = h*vvel; \
  }

#define _ShallowWater2DEigenvalues_(u,D,p,dir) \
  { \
    double h,uvel,vvel,c; \
    _ShallowWater2DGetFlowVar_(u,h,uvel,vvel); \
    c = sqrt(p->g*h); \
    if (dir == _XDIR_) { \
      D[0*_MODEL_NVARS_+0] = uvel-c; D[0*_MODEL_NVARS_+1] = 0;      D[0*_MODEL_NVARS_+2] = 0;      \
      D[1*_MODEL_NVARS_+0] = 0;      D[1*_MODEL_NVARS_+1] = uvel;   D[1*_MODEL_NVARS_+2] = 0;      \
      D[2*_MODEL_NVARS_+0] = 0;      D[2*_MODEL_NVARS_+1] = 0;      D[2*_MODEL_NVARS_+2] = uvel+c; \
    } else if (dir == _YDIR_) { \
      D[0*_MODEL_NVARS_+0] = vvel-c; D[0*_MODEL_NVARS_+1] = 0;      D[0*_MODEL_NVARS_+2] = 0;      \
      D[1*_MODEL_NVARS_+0] = 0;      D[1*_MODEL_NVARS_+1] = vvel;   D[1*_MODEL_NVARS_+2] = 0;      \
      D[2*_MODEL_NVARS_+0] = 0;      D[2*_MODEL_NVARS_+1] = 0;      D[2*_MODEL_NVARS_+2] = vvel+c; \
    } \
  }

#define _ShallowWater2DLeftEigenvectors_(u,L,p,dir) \
  { \
    double h,uvel,vvel,c; \
    _ShallowWater2DGetFlowVar_(u,h,uvel,vvel); \
    c    = sqrt(p->g*h); \
    \
    if (dir == _XDIR_) { \
      L[0*_MODEL_NVARS_+0] = 0.5 + uvel/(2*c); \
      L[0*_MODEL_NVARS_+1] = -1.0/(2*c); \
      L[0*_MODEL_NVARS_+2] = 0.0; \
      \
      L[1*_MODEL_NVARS_+0] = vvel; \
      L[1*_MODEL_NVARS_+1] = 0.0; \
      L[1*_MODEL_NVARS_+2] = -1.0; \
      \
      L[2*_MODEL_NVARS_+0] = 0.5 - uvel/(2*c); \
      L[2*_MODEL_NVARS_+1] = 1.0/(2*c); \
      L[2*_MODEL_NVARS_+2] = 0.0; \
      \
    } else if (dir == _YDIR_) { \
      L[0*_MODEL_NVARS_+0] = 0.5 + vvel/(2*c); \
      L[0*_MODEL_NVARS_+1] = 0.0; \
      L[0*_MODEL_NVARS_+2] = -1.0/(2*c); \
      \
      L[1*_MODEL_NVARS_+0] = -uvel; \
      L[1*_MODEL_NVARS_+1] = 1.0; \
      L[1*_MODEL_NVARS_+2] = 0.0; \
      \
      L[2*_MODEL_NVARS_+0] = 0.5 - vvel/(2*c); \
      L[2*_MODEL_NVARS_+1] = 0.0; \
      L[2*_MODEL_NVARS_+2] = 1.0/(2*c); \
    } \
  }

#define _ShallowWater2DRightEigenvectors_(u,R,p,dir) \
  { \
    double h,uvel,vvel,c; \
    _ShallowWater2DGetFlowVar_(u,h,uvel,vvel); \
    c    = sqrt(p->g*h); \
    \
    if (dir == _XDIR_) { \
      R[0*_MODEL_NVARS_+0] = 1.0; \
      R[1*_MODEL_NVARS_+0] = uvel-c; \
      R[2*_MODEL_NVARS_+0] = vvel; \
      \
      R[0*_MODEL_NVARS_+1] = 0.0; \
      R[1*_MODEL_NVARS_+1] = 0.0; \
      R[2*_MODEL_NVARS_+1] = -1.0;\
      \
      R[0*_MODEL_NVARS_+2] = 1.0; \
      R[1*_MODEL_NVARS_+2] = uvel+c; \
      R[2*_MODEL_NVARS_+2] = vvel; \
      \
    } else if (dir == _YDIR_) { \
      R[0*_MODEL_NVARS_+0] = 1.0; \
      R[1*_MODEL_NVARS_+0] = uvel; \
      R[2*_MODEL_NVARS_+0] = vvel-c; \
      \
      R[0*_MODEL_NVARS_+1] = 0.0; \
      R[1*_MODEL_NVARS_+1] = 1.0; \
      R[2*_MODEL_NVARS_+1] = 0.0;\
      \
      R[0*_MODEL_NVARS_+2] = 1.0; \
      R[1*_MODEL_NVARS_+2] = uvel; \
      R[2*_MODEL_NVARS_+2] = vvel+c; \
    } \
  }

typedef struct shallowwater2d_parameters {
  int     bt_type,   
          topo_flag; 
  double  g,         
          *b,        
          fhat,      
          beta,      
          D;         
  char    upw_choice[_MAX_STRING_SIZE_]; 
  int (*SourceUpwind)(double*,double*,double*,double*,int,void*,double);
} ShallowWater2D;

int    ShallowWater2DInitialize (void*,void*);
int    ShallowWater2DCleanup    (void*);