ShallowWater1DInitialize.c File Reference

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

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

Go to the source code of this file.

Functions
double	ShallowWater1DComputeCFL (void *, void *, double, double)
int	ShallowWater1DFlux (double *, double *, int, void *, double)
int	ShallowWater1DSource (double *, double *, void *, void *, double)
int	ShallowWater1DUpwindRoe (double *, double *, double *, double *, double *, double *, int, void *, double)
int	ShallowWater1DUpwindLLF (double *, double *, double *, double *, double *, double *, int, void *, double)
int	ShallowWater1DJacobian (double *, double *, void *, int, int, int)
int	ShallowWater1DRoeAverage (double *, double *, double *, void *)
int	ShallowWater1DLeftEigenvectors (double *, double *, void *, int)
int	ShallowWater1DRightEigenvectors (double *, double *, void *, int)
int	ShallowWater1DTopography (void *, void *, int, int, int *)
int	ShallowWater1DSourceUpwindLLF (double *, double *, double *, double *, int, void *, double)
int	ShallowWater1DSourceUpwindRoe (double *, double *, double *, double *, int, void *, double)
int	ShallowWater1DModifiedSolution (double *, double *, int, void *, void *, double)
int	ShallowWater1DWriteTopography (void *, void *, double)
int	ShallowWater1DInitialize (void *s, void *m)

Detailed Description

Initialize the 1D shallow water equations module.

Author

Debojyoti Ghosh

Definition in file ShallowWater1DInitialize.c.

Function Documentation

ShallowWater1DComputeCFL()

double ShallowWater1DComputeCFL	(	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 ShallowWater1DComputeCFL.c.

{
  HyPar           *solver = (HyPar*)   s;
  ShallowWater1D  *param  = (ShallowWater1D*) 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, v, c, dxinv, local_cfl;
    _ShallowWater1DGetFlowVar_((u+_MODEL_NVARS_*p),h,v);
    _GetCoordinate_(0,index[0],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
    c = sqrt(param->g*h); /* speed of gravity waves */
    local_cfl = (absolute(v)+c)*dt*dxinv; /* local cfl for this grid point */
    if (local_cfl > max_cfl) max_cfl = local_cfl;
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(max_cfl);
}

ShallowWater1DFlux()

int ShallowWater1DFlux	(	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}{c} hu \\ hu^2 + \frac{1}{2} gh^2 \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 (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current time

Definition at line 16 of file ShallowWater1DFlux.c.

{
  HyPar             *solver = (HyPar*)   s;
  ShallowWater1D    *param  = (ShallowWater1D*) 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, v;
    _ShallowWater1DGetFlowVar_((u+nvars*p),h,v);
    _ShallowWater1DSetFlux_((f+nvars*p),h,v,param->g);
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }

  return(0);
}

ShallowWater1DSource()

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

Compute the source terms for the 1D shallow water equations. The 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

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 25 of file ShallowWater1DSource.c.

{
  HyPar          *solver = (HyPar* ) s;
  MPIVariables   *mpi = (MPIVariables*) m;
  ShallowWater1D *param  = (ShallowWater1D*) 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];

  /* set interface dimensions */
  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_XDIR_]++;
  
  /* calculate the first source function */
  IERR ShallowWater1DSourceFunction1(SourceC,u,x,solver,mpi,t); 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 ShallowWater1DSourceFunction2(SourceC,u,x,solver,mpi,t); 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, vel;      _ShallowWater1DGetFlowVar_((u+_MODEL_NVARS_*p),h,vel);
    double term[_MODEL_NVARS_] = { 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])*dx_inverse;
    }
    vel = h; /* useless statement to avoid compiler warning */
    _ArrayIncrementIndex_(ndims,dim,index,done);
  }

  return(0);
}

ShallowWater1DUpwindRoe()

int ShallowWater1DUpwindRoe	(	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 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 (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 31 of file ShallowWater1DUpwind.c.

{
  HyPar           *solver = (HyPar*) s;
  ShallowWater1D  *param  = (ShallowWater1D*) 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]);

      _ShallowWater1DRoeAverage_         (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param); 
      _ShallowWater1DEigenvalues_        (uavg,D,param,0);
      _ShallowWater1DLeftEigenvectors_   (uavg,L,param,0);
      _ShallowWater1DRightEigenvectors_  (uavg,R,param,0);

      k = 0; D[k] = absolute(D[k]);
      k = 3; D[k] = absolute(D[k]);

      MatMult2(2,DL,D,L);
      MatMult2(2,modA,R,DL);

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

      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];
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

ShallowWater1DUpwindLLF()

int ShallowWater1DUpwindLLF	(	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 (unused since this is a 1D system)
s	Solver object of type HyPar
t	Current solution time

Definition at line 114 of file ShallowWater1DUpwind.c.

{
  HyPar           *solver = (HyPar*) s;
  ShallowWater1D  *param  = (ShallowWater1D*) 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 */

      _ShallowWater1DRoeAverage_       (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
      _ShallowWater1DEigenvalues_      (uavg,D,param,0);
      _ShallowWater1DLeftEigenvectors_ (uavg,L,param,0);
      _ShallowWater1DRightEigenvectors_(uavg,R,param,0);

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

      for (k = 0; k < _MODEL_NVARS_; k++) {
        double eigL,eigC,eigR;
        _ShallowWater1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
        eigL = D[k*_MODEL_NVARS_+k];
        _ShallowWater1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
        eigR = D[k*_MODEL_NVARS_+k];
        _ShallowWater1DEigenvalues_(uavg,D,param,0);
        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 */
      MatVecMult2(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
    }
    _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
  }

  return(0);
}

ShallowWater1DJacobian()

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

Function to compute the flux Jacobian of the 1D 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 = 4
u	solution at a grid point (array of size nvar = 2)
p	object containing the physics-related parameters
dir	dimension (x/y/z) (not used, since this is 1D system)
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 ShallowWater1DJacobian.c.

{
  ShallowWater1D  *param = (ShallowWater1D*) 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 */
  _ShallowWater1DEigenvalues_      (u,D,param,0);
  _ShallowWater1DLeftEigenvectors_ (u,L,param,0);
  _ShallowWater1DRightEigenvectors_(u,R,param,0);

  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 = 3; 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]) );

  MatMult2(2,DL,D,L);
  MatMult2(2,Jac,R,DL);

  return(0);
}

ShallowWater1DRoeAverage()

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

Compute the Roe-averaged state for the 1D shallow water equations. This function just calls the macro _ShallowWater1DRoeAverage_ and is not used by any functions within the 1D 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 ShallowWater1D with physics-related variables

Definition at line 16 of file ShallowWater1DFunctions.c.

{
  ShallowWater1D *param  = (ShallowWater1D*) p;
  _ShallowWater1DRoeAverage_(uavg,uL,uR,param);
  return(0);
}

ShallowWater1DLeftEigenvectors()

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

Compute the left eigenvections for the 1D shallow water equations. This function just calls the macro _ShallowWater1DLeftEigenvectors_ and is not used by any functions within the 1D 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 = 2^2 to save the matrix of left eigenvectors in (row-major format).
p	Object of type ShallowWater1D with physics-related variables
dir	Spatial dimension (not used, since this is a 1D system)

Definition at line 19 of file ShallowWater1DEigen.c.

{
  ShallowWater1D *param  = (ShallowWater1D*)  p;
  _ShallowWater1DLeftEigenvectors_(u,L,param,dir);
  return(0);
}

ShallowWater1DRightEigenvectors()

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

Compute the right eigenvections for the 1D shallow water equations. This function just calls the macro _ShallowWater1DRightEigenvectors_ and is not used by any functions within the 1D 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 = 2^2 to save the matrix of right eigenvectors in (row-major format).
p	Object of type ShallowWater1D with physics-related variables
dir	Spatial dimension (not used, since this is a 1D system)

Definition at line 38 of file ShallowWater1DEigen.c.

{
  ShallowWater1D *param  = (ShallowWater1D*)  p;
  _ShallowWater1DRightEigenvectors_(u,R,param,dir);
  return(0);
}

ShallowWater1DTopography()

int ShallowWater1DTopography	(	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 ShallowWater1DTopography.c.

{
  HyPar          *solver = (HyPar*) s;
  MPIVariables   *mpi    = (MPIVariables*) m;
  ShallowWater1D *param  = (ShallowWater1D*) 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 ShallowWater1DTopography()\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 ShallowWater1DTopography()\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);
}

ShallowWater1DSourceUpwindLLF()

int ShallowWater1DSourceUpwindLLF	(	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 (unused since this is a 1D case)
s	Solver object of type HyPar
t	Current solution time

Definition at line 22 of file ShallowWater1DSourceUpwind.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);
}

ShallowWater1DSourceUpwindRoe()

int ShallowWater1DSourceUpwindRoe	(	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 (unused since this is a 1D case)
s	Solver object of type HyPar
t	Current solution time

Definition at line 64 of file ShallowWater1DSourceUpwind.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);
}

ShallowWater1DModifiedSolution()

int ShallowWater1DModifiedSolution	(	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 1D 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 (unused since this is a 1D system)
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current solution time

Definition at line 21 of file ShallowWater1DModifiedSolution.c.

{
  HyPar           *solver = (HyPar*) s;
  ShallowWater1D  *param  = (ShallowWater1D*) 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,v;
    _ShallowWater1DGetFlowVar_((u+_MODEL_NVARS_*p),h,v);
    uC[_MODEL_NVARS_*p+0] = h + param->b[p];
    uC[_MODEL_NVARS_*p+1] = h * v;
    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }
  return(0);
}

ShallowWater1DWriteTopography()

int ShallowWater1DWriteTopography	(	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 ShallowWater1DWriteTopography.c.

{
  HyPar           *solver = (HyPar*)          s;
  MPIVariables    *mpi    = (MPIVariables*)   m;
  ShallowWater1D  *params = (ShallowWater1D*) 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);
}

ShallowWater1DInitialize()

int ShallowWater1DInitialize	(	void *	s,
		void *	m
	)

Function to initialize the 1D shallow water equations (ShallowWater1D) 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 ShallowWater1DInitialize.c.

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

  static int count = 0;

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

  /* default values */
  physics->g       = 1.0; 
  physics->bt_type = 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,"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_character (physics->upw_choice,_MAX_STRING_SIZE_,0,&mpi->world);  CHECKERR(ierr);
#endif

  /* initializing physical model-specific functions */
  solver->ComputeCFL = ShallowWater1DComputeCFL;
  solver->FFunction  = ShallowWater1DFlux;
  solver->SFunction  = ShallowWater1DSource;
  solver->UFunction  = ShallowWater1DModifiedSolution;
  solver->JFunction  = ShallowWater1DJacobian;
  if      (!strcmp(physics->upw_choice,_ROE_ )) solver->Upwind = ShallowWater1DUpwindRoe;
  else if (!strcmp(physics->upw_choice,_LLF_ )) solver->Upwind = ShallowWater1DUpwindLLF;
  else {
    if (!mpi->rank) fprintf(stderr,"Error in ShallowWater1DInitialize(): %s is not a valid upwinding scheme.\n",
                            physics->upw_choice);
    return(1);
  }
  solver->AveragingFunction     = ShallowWater1DRoeAverage;
  solver->GetLeftEigenvectors   = ShallowWater1DLeftEigenvectors;
  solver->GetRightEigenvectors  = ShallowWater1DRightEigenvectors;
  solver->PhysicsOutput         = ShallowWater1DWriteTopography;
   
  if      (!strcmp(physics->upw_choice,_LLF_ )) physics->SourceUpwind = ShallowWater1DSourceUpwindLLF;
  else if (!strcmp(physics->upw_choice,_ROE_ )) physics->SourceUpwind = ShallowWater1DSourceUpwindRoe;

  /* allocate array to hold the bottom topography field */
  physics->b = (double*) calloc (solver->npoints_local_wghosts, sizeof(double));
  /* set function pointer to read this topography */
  solver->PhysicsInput = ShallowWater1DTopography;

  count++;
  return(0);
}