ShallowWater2DInitialize.c File Reference

Initialize the 2D shallow water equations module. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <physicalmodels/shallowwater2d.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions
double	ShallowWater2DComputeCFL (void *, void *, double, double)
int	ShallowWater2DFlux (double *, double *, int, void *, double)
int	ShallowWater2DSource (double *, double *, void *, void *, double)
int	ShallowWater2DUpwindRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
int	ShallowWater2DUpwindLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	ShallowWater2DJacobian (double *, double *, void *, int, int, int)
int	ShallowWater2DRoeAverage (double *, double *, double *, void *)
int	ShallowWater2DLeftEigenvectors (double *, double *, void *, int)
int	ShallowWater2DRightEigenvectors (double *, double *, void *, int)
int	ShallowWater2DTopography (void *, void *, int, int, int *)
int	ShallowWater2DSourceUpwindLLF (double *, double *, double *, double *, int, void *, double)
int	ShallowWater2DSourceUpwindRoe (double *, double *, double *, double *, int, void *, double)
int	ShallowWater2DModifiedSolution (double *, double *, int, void *, void *, double)
int	ShallowWater2DWriteTopography (void *, void *, double)
int	ShallowWater2DInitialize (void *s, void *m)

Detailed Description

Initialize the 2D shallow water equations module.

Author

Debojyoti Ghosh

Definition in file ShallowWater2DInitialize.c.

Function Documentation

double ShallowWater2DComputeCFL	(	void *	s,
		void *	m,
		double	dt,
		double	t
	)

Computes the maximum CFL number over the domain. Note that the CFL is computed over the local domain on this processor only.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
dt	Time step size for which to compute the CFL
t	Time

Definition at line 17 of file ShallowWater2DComputeCFL.c.

{
  HyPar           *solver = (HyPar*)   s;
  ShallowWater2D  *param  = (ShallowWater2D*) solver->physics;

  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;
  int ndims   = solver->ndims;
  int index[ndims];
  double *u   = solver->u;

  double max_cfl = 0;
  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
    double h, uvel, vvel, c, dxinv, dyinv, local_cfl[_MODEL_NDIMS_];
    _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);

    _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
    _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv); /* 1/dy */
    c = sqrt(param->g*h); /* speed of gravity waves */
    local_cfl[_XDIR_]=(absolute(uvel)+c)*dt*dxinv;/* local cfl for this grid point (x) */
    local_cfl[_YDIR_]=(absolute(vvel)+c)*dt*dyinv;/* local cfl for this grid point (y) */
    if (local_cfl[_XDIR_] > max_cfl) max_cfl = local_cfl[_XDIR_];
    if (local_cfl[_YDIR_] > max_cfl) max_cfl = local_cfl[_YDIR_];

    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(max_cfl);
}

int ShallowWater2DFlux	(	double *	f,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the hyperbolic flux over the local domain.

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

Parameters

f	Array to hold the computed flux (same size and layout as u)
u	Array containing the conserved solution
dir	Spatial dimension
s	Solver object of type HyPar
t	Current time

Definition at line 20 of file ShallowWater2DFlux.c.

{
  HyPar             *solver = (HyPar*)   s;
  ShallowWater2D    *param  = (ShallowWater2D*) solver->physics;
  int               *dim    = solver->dim_local;
  int               ghosts  = solver->ghosts;
  static const int  ndims   = _MODEL_NDIMS_;
  static const int  nvars   = _MODEL_NVARS_;
  static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];

  /* set bounds for array index to include ghost points */
  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    double h, uvel, vvel;
    _ShallowWater2DGetFlowVar_((u+nvars*p),h,uvel,vvel);
    _ShallowWater2DSetFlux_((f+nvars*p),h,uvel,vvel,param->g,dir);
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

int ShallowWater2DSource	(	double *	source,
		double *	u,
		void *	s,
		void *	m,
		double	t
	)

Compute the source terms for the 2D shallow water equations.

• The topography gradient source term is computed according to the balanced formulation introduced in the reference below. The source term is reformulated and "discretized" in a similar fashion as the hyperbolic flux to ensure that the hydrostatic balance is maintained to machine precision.
  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
• The Coriolis source term is added in a naive way – a balanced formulation may be needed. The Coriolis source term is implemented as described in:
  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)

Parameters

source	Computed source terms (array size & layout same as u)
u	Solution (conserved variables)
s	Solver object of type HyPar
m	MPI object of type MPIVariables
t	Current solution time

Definition at line 33 of file ShallowWater2DSource.c.

{
  HyPar          *solver = (HyPar* ) s;
  MPIVariables   *mpi = (MPIVariables*) m;
  ShallowWater2D *param  = (ShallowWater2D*) solver->physics;

  int     v, done, p, p1, p2;
  double  *SourceI = solver->fluxI; /* interace source term       */
  double  *SourceC = solver->fluxC; /* cell-centered source term  */
  double  *SourceL = solver->fL;
  double  *SourceR = solver->fR;

  int     ndims   = solver->ndims;
  int     ghosts  = solver->ghosts;
  int     *dim    = solver->dim_local;
  double  *x      = solver->x;
  double  *dxinv  = solver->dxinv;
  int     index[ndims],index1[ndims],index2[ndims],dim_interface[ndims];

  /* Along X */
  
  /* set interface dimensions */
  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_XDIR_]++;

  /* calculate the first source function */
  IERR ShallowWater2DSourceFunction1(SourceC,u,x,solver,mpi,t,_XDIR_); CHECKERR(ierr);
  /* calculate the left and right interface source terms */
  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
  /* calculate the final interface source term */
  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_XDIR_,solver,t);
  /* calculate the final cell-centered source term */
  done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index,index1,ndims);
    _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
    _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p );
    _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
    _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);
    double dx_inverse;   _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
    for (v=0; v<_MODEL_NVARS_; v++)
      source[_MODEL_NVARS_*p+v] += (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse;
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }
  /* calculate the second source function */
  IERR ShallowWater2DSourceFunction2(SourceC,u,x,solver,mpi,t,_XDIR_); CHECKERR(ierr);
  /* calculate the left and right interface source terms */
  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
  /* calculate the final interface source term */
  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_XDIR_,solver,t);
  /* calculate the final cell-centered source term */
  done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index,index1,ndims);
    _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
    _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p );
    _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
    _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);
    double dx_inverse;  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
    double h, uvel, vvel; _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);
    double term[_MODEL_NVARS_] = { 0.0, -param->g * (h + param->b[p]), 0.0 };
    for (v=0; v<_MODEL_NVARS_; v++)
      source[_MODEL_NVARS_*p+v] += term[v]*(SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse;
    uvel = vvel = h; /* useless statement to avoid compiler warning */
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  /* Along Y */

  /* set interface dimensions */
  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_YDIR_]++;

  /* calculate the first source function */
  IERR ShallowWater2DSourceFunction1(SourceC,u,x,solver,mpi,t,_YDIR_); CHECKERR(ierr);
  /* calculate the left and right interface source terms */
  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
  /* calculate the final interface source term */
  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_YDIR_,solver,t);
  /* calculate the final cell-centered source term */
  done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index,index1,ndims);
    _ArrayCopy1D_(index,index2,ndims); index2[_YDIR_]++;
    _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p );
    _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
    _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);
    double dy_inverse;   _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,dxinv,dy_inverse);
    for (v=0; v<_MODEL_NVARS_; v++)
      source[_MODEL_NVARS_*p+v] += (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dy_inverse;
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }
  /* calculate the second source function */
  IERR ShallowWater2DSourceFunction2(SourceC,u,x,solver,mpi,t,_YDIR_); CHECKERR(ierr);
  /* calculate the left and right interface source terms */
  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
  /* calculate the final interface source term */
  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_YDIR_,solver,t);
  /* calculate the final cell-centered source term */
  done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index,index1,ndims);
    _ArrayCopy1D_(index,index2,ndims); index2[_YDIR_]++;
    _ArrayIndex1D_(ndims,dim          ,index ,ghosts,p );
    _ArrayIndex1D_(ndims,dim_interface,index1,0     ,p1);
    _ArrayIndex1D_(ndims,dim_interface,index2,0     ,p2);
    double dy_inverse;  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,dxinv,dy_inverse);
    double h, uvel, vvel; _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);
    double term[_MODEL_NVARS_] = { 0.0, 0.0, -param->g * (h + param->b[p]) };
    for (v=0; v<_MODEL_NVARS_; v++)
      source[_MODEL_NVARS_*p+v] += term[v]*(SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dy_inverse;
    uvel = vvel = h; /* useless statement to avoid compiler warning */
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  /* Naive addition of Coriolis force terms */
  done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    _ArrayIndex1D_(ndims,dim,index ,ghosts,p);
    double h, uvel, vvel;
    _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);
    double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->x,ycoord);
    double coeff = param->fhat + param->beta * (ycoord - 0.5*param->D);
    double coriolis_x =  coeff * vvel,
           coriolis_y = -coeff * uvel;
    source[_MODEL_NVARS_*p+1] += coriolis_x;
    source[_MODEL_NVARS_*p+2] += coriolis_y;
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(0);
}

int ShallowWater2DUpwindRoe	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Roe's upwinding scheme.

\begin{equation} {\bf f}_{j+1/2} = \frac{1}{2}\left[ {\bf f}_{j+1/2}^L + {\bf f}_{j+1/2}^R - \left| A\left({\bf u}_{j+1/2}^L,{\bf u}_{j+1/2}^R\right) \right| \left( {\bf u}_{j+1/2}^R - {\bf u}_{j+1/2}^L \right)\right] \end{equation}

• Roe, P. L., “Approximate Riemann solvers, parameter vectors, and difference schemes,” Journal of Computational Physics, Vol. 43, No. 2, 1981, pp. 357–372, http://dx.doi.org/10.1016/0021-9991(81)90128-5.

This upwinding scheme is modified for the balanced discretization of the 2D shallow water equations. See the following reference:

• 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

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension
s	Solver object of type HyPar
t	Current solution time

Definition at line 31 of file ShallowWater2DUpwind.c.

{
  HyPar           *solver = (HyPar*) s;
  ShallowWater2D  *param  = (ShallowWater2D*) solver->physics;
  int             done,k;

  int ndims = solver->ndims;
  int ghosts= solver->ghosts;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_], 
                DL[_MODEL_NVARS_*_MODEL_NVARS_], modA[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];

      /* Roe's upwinding scheme */

      udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
      udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
      udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);

      _ShallowWater2DRoeAverage_         (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param); 
      _ShallowWater2DEigenvalues_        (uavg,D,param,dir);
      _ShallowWater2DLeftEigenvectors_   (uavg,L,param,dir);
      _ShallowWater2DRightEigenvectors_  (uavg,R,param,dir);

      k = 0; D[k] = absolute(D[k]);
      k = 4; D[k] = absolute(D[k]);
      k = 8; D[k] = absolute(D[k]);

      MatMult3(_MODEL_NVARS_,DL,D,L);
      MatMult3(_MODEL_NVARS_,modA,R,DL);

      udiss[0] = modA[0*_MODEL_NVARS_+0]*udiff[0] + modA[0*_MODEL_NVARS_+1]*udiff[1] + modA[0*_MODEL_NVARS_+2]*udiff[2];
      udiss[1] = modA[1*_MODEL_NVARS_+0]*udiff[0] + modA[1*_MODEL_NVARS_+1]*udiff[1] + modA[1*_MODEL_NVARS_+2]*udiff[2];
      udiss[2] = modA[2*_MODEL_NVARS_+0]*udiff[0] + modA[2*_MODEL_NVARS_+1]*udiff[1] + modA[2*_MODEL_NVARS_+2]*udiff[2];

      fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
      fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
      fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

int ShallowWater2DUpwindLLF	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	uL,
		double *	uR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Characteristic-based local Lax-Friedrich upwinding scheme.

\begin{align} \alpha_{j+1/2}^{k,L} &= \sum_{k=1}^2 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,L}, \\ \alpha_{j+1/2}^{k,R} &= \sum_{k=1}^2 {\bf l}_{j+1/2}^k \cdot {\bf f}_{j+1/2}^{k,R}, \\ v_{j+1/2}^{k,L} &= \sum_{k=1}^2 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,L}, \\ v_{j+1/2}^{k,R} &= \sum_{k=1}^2 {\bf l}_{j+1/2}^k \cdot {\bf u}_{j+1/2}^{k,R}, \\ \alpha_{j+1/2}^k &= \frac{1}{2}\left[ \alpha_{j+1/2}^{k,L} + \alpha_{j+1/2}^{k,R} - \left(\max_{\left[j,j+1\right]} \lambda\right) \left( v_{j+1/2}^{k,R} - v_{j+1/2}^{k,L} \right) \right], \\ {\bf f}_{j+1/2} &= \sum_{k=1}^2 \alpha_{j+1/2}^k {\bf r}_{j+1/2}^k \end{align}

where \({\bf l}\), \({\bf r}\), and \(\lambda\) are the left-eigenvectors, right-eigenvectors and eigenvalues. The subscripts denote the grid locations.

• C.-W. Shu, and S. Osher, "Efficient implementation of essentially non-oscillatory schemes, II", J. Comput. Phys., 83 (1989), pp. 32–78, http://dx.doi.org/10.1016/0021-9991(89)90222-2.

This upwinding scheme is modified for the balanced discretization of the 1D shallow water equations. See the following reference:

• 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

Parameters

fI	Computed upwind interface flux
fL	Left-biased reconstructed interface flux
fR	Right-biased reconstructed interface flux
uL	Left-biased reconstructed interface solution
uR	Right-biased reconstructed interface solution
u	Cell-centered solution
dir	Spatial dimension
s	Solver object of type HyPar
t	Current solution time

Definition at line 118 of file ShallowWater2DUpwind.c.

{
  HyPar           *solver = (HyPar*) s;
  ShallowWater2D  *param  = (ShallowWater2D*) solver->physics;
  int             done,k;

  int ndims   = solver->ndims;
  int *dim    = solver->dim_local;
  int ghosts  = solver->ghosts;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
      int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
      int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
      int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
      double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_], 
             ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];

      /* Local Lax-Friedrich upwinding scheme */

      _ShallowWater2DRoeAverage_       (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
      _ShallowWater2DEigenvalues_      (uavg,D,param,dir);
      _ShallowWater2DLeftEigenvectors_ (uavg,L,param,dir);
      _ShallowWater2DRightEigenvectors_(uavg,R,param,dir);

      /* calculate characteristic fluxes and variables */
      MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
      MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));

      for (k = 0; k < _MODEL_NVARS_; k++) {
        double eigL,eigC,eigR;
        _ShallowWater2DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,dir);
        eigL = D[k*_MODEL_NVARS_+k];
        _ShallowWater2DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,dir);
        eigR = D[k*_MODEL_NVARS_+k];
        _ShallowWater2DEigenvalues_(uavg,D,param,dir);
        eigC = D[k*_MODEL_NVARS_+k];

        double alpha = max3(absolute(eigL),absolute(eigC),absolute(eigR));
        fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
      }

      /* calculate the interface flux from the characteristic flux */
      MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

int ShallowWater2DJacobian	(	double *	Jac,
		double *	u,
		void *	p,
		int	dir,
		int	nvars,
		int	upw
	)

Function to compute the flux Jacobian of the 2D shallow water equations, given the solution at a grid point. The Jacobian is square matrix of size nvar=2, and is returned as a 1D array (double) of 4 elements in row-major format.

Parameters

Jac	Jacobian matrix: 1D array of size nvar^2 = 9
u	solution at a grid point (array of size nvar = 3)
p	object containing the physics-related parameters
dir	spatial dimension (x/y)
nvars	number of vector components
upw	0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 14 of file ShallowWater2DJacobian.c.

{
  ShallowWater2D  *param = (ShallowWater2D*) p;
  static double   R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], 
                  L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_];

  /* get the eigenvalues and left,right eigenvectors */
  _ShallowWater2DEigenvalues_      (u,D,param,dir);
  _ShallowWater2DLeftEigenvectors_ (u,L,param,dir);
  _ShallowWater2DRightEigenvectors_(u,R,param,dir);

  int aupw = absolute(upw), k;
  k = 0; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 4; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 8; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );

  MatMult3(_MODEL_NVARS_,DL,D,L);
  MatMult3(_MODEL_NVARS_,Jac,R,DL);

  return(0);
}

int ShallowWater2DRoeAverage	(	double *	uavg,
		double *	uL,
		double *	uR,
		void *	p
	)

Compute the Roe-averaged state for the 2D shallow water equations. This function just calls the macro _ShallowWater2DRoeAverage_ and is not used by any functions within the 2D shallow water module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

uavg	The computed Roe-averaged state
uL	Left state (conserved variables)
uR	Right state (conserved variables)
p	Object of type ShallowWater2D with physics-related variables

Definition at line 16 of file ShallowWater2DFunctions.c.

{
  ShallowWater2D *param  = (ShallowWater2D*) p;
  _ShallowWater2DRoeAverage_(uavg,uL,uR,param);
  return(0);
}

int ShallowWater2DLeftEigenvectors	(	double *	u,
		double *	L,
		void *	p,
		int	dir
	)

Compute the left eigenvections for the 2D shallow water equations. This function just calls the macro _ShallowWater2DLeftEigenvectors_ and is not used by any functions within the 2D shallow water module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

u	Conserved solution at a grid point
L	Array of size nvar^2 = 3^2 to save the matrix of left eigenvectors in (row-major format).
p	Object of type ShallowWater2D with physics-related variables
dir	Spatial dimension

Definition at line 19 of file ShallowWater2DEigen.c.

{
  ShallowWater2D *param  = (ShallowWater2D*)  p;
  _ShallowWater2DLeftEigenvectors_(u,L,param,dir);
  return(0);
}

int ShallowWater2DRightEigenvectors	(	double *	u,
		double *	R,
		void *	p,
		int	dir
	)

Compute the right eigenvections for the 2D shallow water equations. This function just calls the macro _ShallowWater2DRightEigenvectors_ and is not used by any functions within the 2D shallow water module. However, it's necessary to define it and provide it to the the solver object (HyPar) so that it can then send it to interpolation functions for a characteristic-based reconstruction.

Parameters

u	Conserved solution at a grid point
R	Array of size nvar^2 = 3^2 to save the matrix of right eigenvectors in (row-major format).
p	Object of type ShallowWater2D with physics-related variables
dir	Spatial dimension

Definition at line 38 of file ShallowWater2DEigen.c.

{
  ShallowWater2D *param  = (ShallowWater2D*)  p;
  _ShallowWater2DRightEigenvectors_(u,R,param,dir);
  return(0);
}

int ShallowWater2DTopography	(	void *	s,
		void *	m,
		int	idx,
		int	nsims,
		int *	dim_topo
	)

Set the bottom topography over the domain - reads the topography data from a file, if available, else sets it to a constant

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
idx	Index of this simulation
nsims	Total number of simulations
dim_topo	Grid dimensions of the advection field stored in file

Definition at line 20 of file ShallowWater2DTopography.c.

{
  HyPar          *solver = (HyPar*) s;
  MPIVariables   *mpi    = (MPIVariables*) m;
  ShallowWater2D *param  = (ShallowWater2D*) solver->physics;
  double         *S      = param->b;
  int            d, done, *dim = solver->dim_local, 
                 ghosts = solver->ghosts;
  _DECLARE_IERR_;

  char fname_root[_MAX_STRING_SIZE_] = "topography";
  if (idx >= 0) {
    if (nsims > 1) {
      char index[_MAX_STRING_SIZE_];
      GetStringFromInteger(idx, index, (int)log10(nsims)+1);
      strcat(fname_root, "_");
      strcat(fname_root, index);
    }
  }

  /* read topography from provided file, if available */
  if (dim_topo == NULL) {
    int ierr = ReadArray( solver->ndims,
                          1,
                          solver->dim_global,
                          solver->dim_local,
                          solver->ghosts,
                          solver,
                          mpi,
                          NULL,
                          S,
                          fname_root,
                          &param->topo_flag);
    if (ierr) {
      fprintf(stderr,"Error in ShallowWater2DTopography()\n");
      fprintf(stderr,"  ReadArray() returned error!\n");
      return ierr;
    }
  } else {
    int ierr = ReadArraywInterp(  solver->ndims,
                                  1,
                                  solver->dim_global,
                                  solver->dim_local,
                                  dim_topo,
                                  solver->ghosts,
                                  solver,
                                  mpi,
                                  NULL,
                                  S,
                                  fname_root,
                                  &param->topo_flag);
    if (ierr) {
      fprintf(stderr,"Error in ShallowWater2DTopography()\n");
      fprintf(stderr,"  ReadArraywInterp() returned error!\n");
      return ierr;
    }
  }

  if (!param->topo_flag) {
    /* if topography file not available, set it to zero */
    _ArraySetValue_(S,solver->npoints_local_wghosts,0.0);
  }

  /* if parallel, exchange MPI-boundary ghost point data */
  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,1,solver->dim_local,
                                solver->ghosts,mpi,S); CHECKERR(ierr);


  if (param->bt_type) {
    /* if topography is periodic, then the overall problem must also be periodic
       (i.e. boundary conditions will be specified as periodic). Hence, 
       MPIExchangeBoundariesnD() will take care of setting the ghosts points 
       for multi-processor simulations. For single processor, set the ghost
       points accordingly. */
    int indexb[_MODEL_NDIMS_], indexi[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], 
        offset[_MODEL_NDIMS_];
    for (d = 0; d < _MODEL_NDIMS_; d++) {
      if (mpi->iproc[d] == 1) {
        _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
        /* left boundary */
        done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
        _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = -ghosts;
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = indexb[d] + dim[d] - ghosts;
          int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
          S[p1] = S[p2];
          _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
        }
        /* right boundary */
        done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
        _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = dim[d];
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_);
          int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
          S[p1] = S[p2];
          _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
        }
      }
    }
  } else {
    /* if topography is not periodic, extrapolate it at the boundaries */
    int indexb[_MODEL_NDIMS_], indexi[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], 
        offset[_MODEL_NDIMS_];
    for (d = 0; d < _MODEL_NDIMS_; d++) {
      /* left boundary */
      if (!mpi->ip[d]) {
        _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
        _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = -ghosts;
        done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = ghosts-1-indexb[d];
          int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
          S[p1] = S[p2];
          _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
        }
      }
      /* right boundary */
      if (mpi->ip[d] == mpi->iproc[d]-1) {
        _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); bounds[d] = ghosts;
        _ArraySetValue_(offset,_MODEL_NDIMS_,0); offset[d] = dim[d];
        done = 0; _ArraySetValue_(indexb,_MODEL_NDIMS_,0);
        while (!done) {
          _ArrayCopy1D_(indexb,indexi,_MODEL_NDIMS_); indexi[d] = dim[d]-1-indexb[d];
          int p1; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,indexb,offset,ghosts,p1);
          int p2; _ArrayIndex1D_  (_MODEL_NDIMS_,dim,indexi,ghosts,p2);
          S[p1] = S[p2];
          _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,indexb,done);
        }
      }
    }
  }

  return(0);
}

int ShallowWater2DSourceUpwindLLF	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the "upwind" source term in the balanced formulation introduced in the reference below. The "upwind" state is just the arithmetic average of the left and right states.

• 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

Parameters

fI	Computed interface source term ("upwinded")
fL	Left-biased interface source term
fR	Right-biased interface source term
u	Solution (conserved variables)
dir	Spatial dimension
s	Solver object of type HyPar
t	Current solution time

Definition at line 22 of file ShallowWater2DSourceUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  int       done,k;

  int ndims = solver->ndims;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p;  _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      for (k = 0; k < _MODEL_NVARS_; k++) 
        (fI+_MODEL_NVARS_*p)[k] = 0.5 * ((fL+_MODEL_NVARS_*p)[k] + (fR+_MODEL_NVARS_*p)[k]);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

int ShallowWater2DSourceUpwindRoe	(	double *	fI,
		double *	fL,
		double *	fR,
		double *	u,
		int	dir,
		void *	s,
		double	t
	)

Compute the "upwind" source term in the balanced formulation introduced in the reference below. The "upwind" state is just the arithmetic average of the left and right states.

• 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

Parameters

fI	Computed interface source term ("upwinded")
fL	Left-biased interface source term
fR	Right-biased interface source term
u	Solution (conserved variables)
dir	Spatial dimension
s	Solver object of type HyPar
t	Current solution time

Definition at line 64 of file ShallowWater2DSourceUpwind.c.

{
  HyPar     *solver = (HyPar*)    s;
  int       done,k;

  int ndims = solver->ndims;
  int *dim  = solver->dim_local;

  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;

  done = 0; _ArraySetValue_(index_outer,ndims,0);
  while (!done) {
    _ArrayCopy1D_(index_outer,index_inter,ndims);
    for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
      int p;  _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
      for (k = 0; k < _MODEL_NVARS_; k++) 
        (fI+_MODEL_NVARS_*p)[k] = 0.5 * ((fL+_MODEL_NVARS_*p)[k] + (fR+_MODEL_NVARS_*p)[k]);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

int ShallowWater2DModifiedSolution	(	double *	uC,
		double *	u,
		int	d,
		void *	s,
		void *	m,
		double	waqt
	)

Compute the modified solution for the upwinding step in a balanced conservative finite-difference algorithm for the 2D shallow water equations.

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

Parameters

uC	The modified solution (same array size and layout as u)
u	The solution (conserved variables)
d	Spatial dimension
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current solution time

Definition at line 21 of file ShallowWater2DModifiedSolution.c.

{
  HyPar           *solver = (HyPar*) s;
  ShallowWater2D  *param  = (ShallowWater2D*) solver->physics;

  int     ghosts  = solver->ghosts;
  int     *dim    = solver->dim_local;
  int     ndims   = solver->ndims;
  int     index[ndims], bounds[ndims], offset[ndims];

  /* set bounds for array index to include ghost points */
  _ArrayCopy1D_(dim,bounds,ndims);
  int i; for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;

  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    double h,uvel,vvel;
    _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);
    uC[_MODEL_NVARS_*p+0] = h + param->b[p];
    uC[_MODEL_NVARS_*p+1] = h * uvel;
    uC[_MODEL_NVARS_*p+2] = h * vvel;
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }
  return(0);
}

int ShallowWater2DWriteTopography	(	void *	s,
		void *	m,
		double	a_t
	)

Write out the topography data to file

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
a_t	Current simulation time

Definition at line 17 of file ShallowWater2DWriteTopography.c.

{
  HyPar           *solver = (HyPar*)          s;
  MPIVariables    *mpi    = (MPIVariables*)   m;
  ShallowWater2D  *params = (ShallowWater2D*) solver->physics;
  _DECLARE_IERR_;

  if (params->topo_flag) {

    char fname_root[_MAX_STRING_SIZE_] = "topography";
    if (solver->nsims > 1) {
      char index[_MAX_STRING_SIZE_];
      GetStringFromInteger(solver->my_idx, index, (int)log10(solver->nsims)+1);
      strcat(fname_root, "_");
      strcat(fname_root, index);
      strcat(fname_root, "_");
    }

    WriteArray(  solver->ndims,
                 1,
                 solver->dim_global,
                 solver->dim_local,
                 solver->ghosts,
                 solver->x,
                 params->b,
                 solver,
                 mpi,
                 fname_root );

    if (!strcmp(solver->plot_solution, "yes")) {
      PlotArray( solver->ndims,
                 1,
                 solver->dim_global,
                 solver->dim_local,
                 solver->ghosts,
                 solver->x,
                 params->b,
                 a_t,
                 solver,
                 mpi,
                 fname_root );
    }

  }

  return(0);
}

int ShallowWater2DInitialize	(	void *	s,
		void *	m
	)

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