shallowwater2d.h File Reference

2D Shallow Water Equations More...

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

Go to the source code of this file.

Data Structures
struct	ShallowWater2D
	Structure containing variables and parameters specific to the 2D Shallow Water equations. This structure contains the physical parameters, variables, and function pointers specific to the 2D ShallowWater equations. More...

Macros
#define	_SHALLOW_WATER_2D_   "shallow-water-2d"
#define	_MODEL_NDIMS_   2
#define	_MODEL_NVARS_   3
#define	_ROE_   "roe"
#define	_LLF_   "llf-char"
#define	_XDIR_   0
#define	_YDIR_   1
#define	_ShallowWater2DGetFlowVar_(u, h, uvel, vvel)
#define	_ShallowWater2DSetFlux_(f, h, uvel, vvel, g, dir)
#define	_ShallowWater2DRoeAverage_(uavg, uL, uR, p)
#define	_ShallowWater2DEigenvalues_(u, D, p, dir)
#define	_ShallowWater2DLeftEigenvectors_(u, L, p, dir)
#define	_ShallowWater2DRightEigenvectors_(u, R, p, dir)

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

Detailed Description

2D Shallow Water Equations

Author

Debojyoti Ghosh

2D Shallow Water Equations

\begin{equation} \frac {\partial} {\partial t} \left[\begin{array}{c} h \\ hu \\ hv \end{array} \right] + \frac {\partial} {\partial x} \left[\begin{array}{c} hu \\ hu^2 + \frac{1}{2}gh^2 \\ huv \end{array} \right] + \frac {\partial} {\partial x} \left[\begin{array}{c} hv \\ huv \\ hv^2 + \frac{1}{2}gh^2 \end{array} \right] = \left[\begin{array}{c} 0 \\ -ghb_x \\ -ghb_y \end{array}\right] + \left[\begin{array}{c} 0 \\ fv \\ -fu \end{array}\right] \end{equation}

where \(h\) is the height, \(\left(u,v\right)\) are the velocity components, \(b(x,y)\) is the bottom topography, and \(g\) is the gravitational constant.

The Coriolis parameter \(f\) is defined as

\begin{equation} f = \hat{f} + \beta\left( y - \frac{D}{2} \right) \end{equation}

based on:

• Zhu, Et. al., "Variational Data Assimilation with a Variable Resolution Finite-Element Shallow-Water Equations Model", Monthly Weather Review, 122, 1994, pp. 946–965 http://dx.doi.org/10.1175/1520-0493(1994)122%3C0946:VDAWAV%3E2.0.CO;2
  Eqns. (2.1)-(2.3), (2.4)

For the treatment of the topography gradient source term (well-balanced formulation), refer to:

• Xing, Y., Shu, C.-W., "High order finite difference WENO schemes with the exact conservation property for the shallow water equations", Journal of Computational Physics, 208, 2005, pp. 206-227. http://dx.doi.org/10.1016/j.jcp.2005.02.006

Definition in file shallowwater2d.h.

Macro Definition Documentation

_SHALLOW_WATER_2D_

#define _SHALLOW_WATER_2D_   "shallow-water-2d"

2D Shallow Water equations

Definition at line 43 of file shallowwater2d.h.

_MODEL_NDIMS_

#define _MODEL_NDIMS_   2

Number of spatial dimensions

Definition at line 49 of file shallowwater2d.h.

_MODEL_NVARS_

#define _MODEL_NVARS_   3

Number of variables per grid point

Definition at line 51 of file shallowwater2d.h.

_ROE_

#define _ROE_   "roe"

Roe upwinding scheme

Definition at line 55 of file shallowwater2d.h.

_LLF_

#define _LLF_   "llf-char"

Local Lax-Friedrich upwinding scheme

Definition at line 57 of file shallowwater2d.h.

_XDIR_

#define _XDIR_   0

dimension corresponding to the x spatial dimension

Definition at line 63 of file shallowwater2d.h.

_YDIR_

#define _YDIR_   1

Definition at line 64 of file shallowwater2d.h.

_ShallowWater2DGetFlowVar_

#define _ShallowWater2DGetFlowVar_	(		u,
			h,
			uvel,
			vvel
	)

Value:

{ \
    h    = u[0]; \
    uvel = u[1] / h; \
    vvel = u[2] / h; \
  }

Compute the flow variables (height, velocity) from the conserved solution vector.

\begin{equation} {\bf u} = \left[\begin{array}{c} h \\ hu \\ hv \end{array}\right] \end{equation}

Definition at line 73 of file shallowwater2d.h.

_ShallowWater2DSetFlux_

#define _ShallowWater2DSetFlux_	(		f,
			h,
			uvel,
			vvel,
			g,
			dir
	)

Value:

{ \
    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); \
    } \
  }

Set the flux vector given the flow variables (height, velocity).

\begin{equation} {\bf F}\left({\bf u}\right) = \left\{\begin{array}{cc} \left[\begin{array}{c} hu \\ hu^2 + \frac{1}{2} gh^2 \\ huv \end{array}\right] & {\rm dir} = x \\ \left[\begin{array}{c} hv \\ huv \\ hv^2 + \frac{1}{2} gh^2 \end{array}\right] & {\rm dir} = y \\ \end{array}\right. \end{equation}

Definition at line 91 of file shallowwater2d.h.

_ShallowWater2DRoeAverage_

#define _ShallowWater2DRoeAverage_	(		uavg,
			uL,
			uR,
			p
	)

Value:

{ \
    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; \
  }

Compute the Roe average of two conserved solution vectors.

Definition at line 108 of file shallowwater2d.h.

_ShallowWater2DEigenvalues_

#define _ShallowWater2DEigenvalues_	(		u,
			D,
			p,
			dir
	)

Value:

{ \
    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; \
    } \
  }

Compute the eigenvalues, given the conserved solution vector. The eigenvalues are returned as a 3x3 matrix stored in row-major format. It is a diagonal matrix with the eigenvalues as diagonal elements. Admittedly, it is a wasteful way of storing the eigenvalues.

Definition at line 129 of file shallowwater2d.h.

_ShallowWater2DLeftEigenvectors_

#define _ShallowWater2DLeftEigenvectors_	(		u,
			L,
			p,
			dir
	)

Compute the matrix that has the left-eigenvectors as its rows. Stored in row-major format. Reference:

• Ambrosi, D., "Approximation of Shallow Water Equations by Roe's Riemann Solver", Intl. J. Num. Meth. Fluids, 20 (1995), pp. 157 – 168, Page 160, Eqn. 15
  http://dx.doi.org/10.1002/fld.1650200205

Definition at line 154 of file shallowwater2d.h.

_ShallowWater2DRightEigenvectors_

#define _ShallowWater2DRightEigenvectors_	(		u,
			R,
			p,
			dir
	)

Compute the matrix that has the right-eigenvectors as its columns. Stored in row-major format.
Reference:

• Ambrosi, D., "Approximation of Shallow Water Equations by Roe's Riemann Solver", Intl. J. Num. Meth. Fluids, 20 (1995), pp. 157 – 168, Page 160, Eqn. 14
  http://dx.doi.org/10.1002/fld.1650200205

Definition at line 197 of file shallowwater2d.h.

Function Documentation

ShallowWater2DInitialize()

int ShallowWater2DInitialize	(	void *	s,
		void *	m
	)

Function to initialize the 2D ShallowWater module

Function to initialize the 2D shallow water equations (ShallowWater2D) module: Sets the default parameters, read in and set physics-related parameters, and set the physics-related function pointers in HyPar.

Parameters

s	Solver object of type HyPar
m	Object of type MPIVariables containing MPI-related info

Definition at line 37 of file ShallowWater2DInitialize.c.

{
  HyPar           *solver  = (HyPar*)         s;
  MPIVariables    *mpi     = (MPIVariables*)  m; 
  ShallowWater2D  *physics = (ShallowWater2D*)       solver->physics;
  int             ferr, d;

  static int count = 0;

  if (solver->nvars != _MODEL_NVARS_) {
    fprintf(stderr,"Error in ShallowWater2DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
    return(1);
  }
  if (solver->ndims != _MODEL_NDIMS_) {
    fprintf(stderr,"Error in ShallowWater2DInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
    return(1);
  }

  /* default values */
  physics->g       = 1.0; 
  physics->bt_type = 0;
  physics->fhat    = 0.0;
  physics->beta    = 0.0;
  physics->D       = 0.0;
  strcpy(physics->upw_choice,"roe");

  /* reading physical model specific inputs */
  if (!mpi->rank) {
    FILE *in;
    if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
    in = fopen("physics.inp","r");
    if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
    else {
      char word[_MAX_STRING_SIZE_];
      ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
      if (!strcmp(word, "begin")){
          while (strcmp(word, "end")){
              ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
          if (!strcmp(word, "gravity")) { 
            ferr = fscanf(in,"%lf",&physics->g); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "topography_type")) { 
            ferr = fscanf(in,"%d",&physics->bt_type); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "coriolis_fhat")) { 
            ferr = fscanf(in,"%lf",&physics->fhat); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "coriolis_beta")) { 
            ferr = fscanf(in,"%lf",&physics->beta); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word, "coriolis_D")) { 
            ferr = fscanf(in,"%lf",&physics->D); 
            if (ferr != 1) return(1);
          } else if (!strcmp(word,"upwinding")) {
            ferr = fscanf(in,"%s",physics->upw_choice); 
            if (ferr != 1) return(1);
          } else if (strcmp(word,"end")) {
            char useless[_MAX_STRING_SIZE_];
            ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
            printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
            printf("recognized or extraneous. Ignoring.\n");
          }
        }
        } else {
          fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
        return(1);
        }
    }
    fclose(in);
  }

#ifndef serial
  IERR MPIBroadcast_double    (&physics->g        ,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_integer   (&physics->bt_type  ,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->fhat     ,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->beta     ,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_double    (&physics->D        ,1,0,&mpi->world);                  CHECKERR(ierr);
  IERR MPIBroadcast_character (physics->upw_choice,_MAX_STRING_SIZE_,0,&mpi->world);  CHECKERR(ierr);
#endif

  /* initializing physical model-specific functions */
  solver->ComputeCFL = ShallowWater2DComputeCFL;
  solver->FFunction  = ShallowWater2DFlux;
  solver->SFunction  = ShallowWater2DSource;
  solver->UFunction  = ShallowWater2DModifiedSolution;
  solver->JFunction  = ShallowWater2DJacobian;
  if      (!strcmp(physics->upw_choice,_ROE_ )) solver->Upwind = ShallowWater2DUpwindRoe;
  else if (!strcmp(physics->upw_choice,_LLF_ )) solver->Upwind = ShallowWater2DUpwindLLF;
  else {
    if (!mpi->rank) fprintf(stderr,"Error in ShallowWater2DInitialize(): %s is not a valid upwinding scheme.\n",
                            physics->upw_choice);
    return(1);
  }
  solver->AveragingFunction     = ShallowWater2DRoeAverage;
  solver->GetLeftEigenvectors   = ShallowWater2DLeftEigenvectors;
  solver->GetRightEigenvectors  = ShallowWater2DRightEigenvectors;
  solver->PhysicsOutput         = ShallowWater2DWriteTopography;
   
  if      (!strcmp(physics->upw_choice,_LLF_ )) physics->SourceUpwind = ShallowWater2DSourceUpwindLLF;
  else if (!strcmp(physics->upw_choice,_ROE_ )) physics->SourceUpwind = ShallowWater2DSourceUpwindRoe;

  /* allocate array to hold the bottom topography field */
  physics->b = (double*) calloc (solver->npoints_local_wghosts, sizeof(double));
  solver->PhysicsInput = ShallowWater2DTopography;

  count++;
  return(0);
}

ShallowWater2DCleanup()

int ShallowWater2DCleanup

(

void *

s

)

Function to clean up the 2D ShallowWater module

Function to clean up all physics-related allocations for the 2D shallow water equations

Parameters

s	Solver object of type HyPar

Definition at line 10 of file ShallowWater2DCleanup.c.

{
  ShallowWater2D *param  = (ShallowWater2D*) s;
  free(param->b);
  return(0);
}