HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ShallowWater2DSource.c File Reference

Contains the functions to compute the source terms for the 2D shallow water equations. More...

#include <stdlib.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

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

Detailed Description

Contains the functions to compute the source terms for the 2D shallow water equations.

Author
Debojyoti Ghosh

Definition in file ShallowWater2DSource.c.

Function Documentation

int ShallowWater2DSourceFunction1 ( double *  f,
double *  u,
double *  x,
void *  s,
void *  m,
double  t,
int  dir 
)
static

Compute the first source function that is then "discretized" in a way similar to the hyperbolic flux function for 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
    This function computes:

    \begin{equation} {\bf S}_2 = \left\{ \begin{array}{cc} \left[ \begin{array}{c} 0 \\ \frac{1}{2}gb^2 \\ 0 \end{array}\right] & {\rm dir} = x \\ \left[ \begin{array}{c} 0 \\ 0 \\ \frac{1}{2}gb^2 \end{array}\right] & {\rm dir} = y \end{array} \right. \end{equation}

Parameters
fComputed source function (array size and layout same as u)
uSolution (conserved variables)
xSpatial coordinates
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent solution time
dirSpatial dimension

Definition at line 194 of file ShallowWater2DSource.c.

203 {
204  HyPar *solver = (HyPar* ) s;
205  ShallowWater2D *param = (ShallowWater2D*) solver->physics;
206 
207  int ghosts = solver->ghosts;
208  int *dim = solver->dim_local;
209  int ndims = solver->ndims;
210  int index[ndims], bounds[ndims], offset[ndims];
211 
212  /* set bounds for array index to include ghost points */
213  _ArrayCopy1D_(dim,bounds,ndims);
214  int i; for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
215 
216  /* set offset such that index is compatible with ghost point arrangement */
217  _ArraySetValue_(offset,ndims,-ghosts);
218 
219  int done = 0; _ArraySetValue_(index,ndims,0);
220  while (!done) {
221  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
222  (f+_MODEL_NVARS_*p)[0] = 0.0;
223  (f+_MODEL_NVARS_*p)[1] = (dir == _XDIR_ ? 0.5 * param->g * param->b[p] * param->b[p] : 0.0 );
224  (f+_MODEL_NVARS_*p)[2] = (dir == _YDIR_ ? 0.5 * param->g * param->b[p] * param->b[p] : 0.0 );
225  _ArrayIncrementIndex_(ndims,bounds,index,done);
226  }
227 
228  return(0);
229 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::ndims
int ndims
Definition: hypar.h:26
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.
Definition: shallowwater2d.h:240
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
int ShallowWater2DSourceFunction2 ( double *  f,
double *  u,
double *  x,
void *  s,
void *  m,
double  t,
int  dir 
)
static

Compute the second source function that is then "discretized" in a way similar to the hyperbolic flux function for 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
    This function computes:

    \begin{equation} {\bf S}_2 = \left\{ \begin{array}{cc} \left[ \begin{array}{c} 0 \\ b \\ 0 \end{array}\right] & {\rm dir} = x \\ \left[ \begin{array}{c} 0 \\ 0 \\ b \end{array}\right] & {\rm dir} = y \end{array} \right. \end{equation}

Parameters
fComputed source function (array size and layout same as u)
uSolution (conserved variables)
xSpatial coordinates
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent solution time
dirSpatial dimension

Definition at line 251 of file ShallowWater2DSource.c.

260 {
261  HyPar *solver = (HyPar* ) s;
262  ShallowWater2D *param = (ShallowWater2D*) solver->physics;
263 
264  int ghosts = solver->ghosts;
265  int *dim = solver->dim_local;
266  int ndims = solver->ndims;
267  int index[ndims], bounds[ndims], offset[ndims];
268 
269  /* set bounds for array index to include ghost points */
270  _ArrayCopy1D_(dim,bounds,ndims);
271  int i; for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
272 
273  /* set offset such that index is compatible with ghost point arrangement */
274  _ArraySetValue_(offset,ndims,-ghosts);
275 
276  int done = 0; _ArraySetValue_(index,ndims,0);
277  while (!done) {
278  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
279  (f+_MODEL_NVARS_*p)[0] = 0.0;
280  (f+_MODEL_NVARS_*p)[1] = (dir == _XDIR_ ? param->b[p] : 0.0 );
281  (f+_MODEL_NVARS_*p)[2] = (dir == _YDIR_ ? param->b[p] : 0.0 );
282  _ArrayIncrementIndex_(ndims,bounds,index,done);
283  }
284 
285  return(0);
286 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1DWO_
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
Definition: arrayfunctions.h:83
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::ndims
int ndims
Definition: hypar.h:26
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.
Definition: shallowwater2d.h:240
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
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
sourceComputed source terms (array size & layout same as u)
uSolution (conserved variables)
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent solution time

Definition at line 33 of file ShallowWater2DSource.c.

40 {
41  HyPar *solver = (HyPar* ) s;
42  MPIVariables *mpi = (MPIVariables*) m;
43  ShallowWater2D *param = (ShallowWater2D*) solver->physics;
44 
45  int v, done, p, p1, p2;
46  double *SourceI = solver->fluxI; /* interace source term */
47  double *SourceC = solver->fluxC; /* cell-centered source term */
48  double *SourceL = solver->fL;
49  double *SourceR = solver->fR;
50 
51  int ndims = solver->ndims;
52  int ghosts = solver->ghosts;
53  int *dim = solver->dim_local;
54  double *x = solver->x;
55  double *dxinv = solver->dxinv;
56  int index[ndims],index1[ndims],index2[ndims],dim_interface[ndims];
57 
58  /* Along X */
59 
60  /* set interface dimensions */
61  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_XDIR_]++;
62 
63  /* calculate the first source function */
64  IERR ShallowWater2DSourceFunction1(SourceC,u,x,solver,mpi,t,_XDIR_); CHECKERR(ierr);
65  /* calculate the left and right interface source terms */
66  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
67  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
68  /* calculate the final interface source term */
69  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_XDIR_,solver,t);
70  /* calculate the final cell-centered source term */
71  done = 0; _ArraySetValue_(index,ndims,0);
72  while (!done) {
73  _ArrayCopy1D_(index,index1,ndims);
74  _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
75  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
76  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
77  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
78  double dx_inverse; _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
79  for (v=0; v<_MODEL_NVARS_; v++)
80  source[_MODEL_NVARS_*p+v] += (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse;
81  _ArrayIncrementIndex_(ndims,dim,index,done);
82  }
83  /* calculate the second source function */
84  IERR ShallowWater2DSourceFunction2(SourceC,u,x,solver,mpi,t,_XDIR_); CHECKERR(ierr);
85  /* calculate the left and right interface source terms */
86  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
87  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_XDIR_,solver,mpi,0); CHECKERR(ierr);
88  /* calculate the final interface source term */
89  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_XDIR_,solver,t);
90  /* calculate the final cell-centered source term */
91  done = 0; _ArraySetValue_(index,ndims,0);
92  while (!done) {
93  _ArrayCopy1D_(index,index1,ndims);
94  _ArrayCopy1D_(index,index2,ndims); index2[_XDIR_]++;
95  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
96  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
97  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
98  double dx_inverse; _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,dxinv,dx_inverse);
99  double h, uvel, vvel; _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);
100  double term[_MODEL_NVARS_] = { 0.0, -param->g * (h + param->b[p]), 0.0 };
101  for (v=0; v<_MODEL_NVARS_; v++)
102  source[_MODEL_NVARS_*p+v] += term[v]*(SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dx_inverse;
103  uvel = vvel = h; /* useless statement to avoid compiler warning */
104  _ArrayIncrementIndex_(ndims,dim,index,done);
105  }
106 
107  /* Along Y */
108 
109  /* set interface dimensions */
110  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[_YDIR_]++;
111 
112  /* calculate the first source function */
113  IERR ShallowWater2DSourceFunction1(SourceC,u,x,solver,mpi,t,_YDIR_); CHECKERR(ierr);
114  /* calculate the left and right interface source terms */
115  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
116  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
117  /* calculate the final interface source term */
118  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_YDIR_,solver,t);
119  /* calculate the final cell-centered source term */
120  done = 0; _ArraySetValue_(index,ndims,0);
121  while (!done) {
122  _ArrayCopy1D_(index,index1,ndims);
123  _ArrayCopy1D_(index,index2,ndims); index2[_YDIR_]++;
124  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
125  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
126  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
127  double dy_inverse; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,dxinv,dy_inverse);
128  for (v=0; v<_MODEL_NVARS_; v++)
129  source[_MODEL_NVARS_*p+v] += (SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dy_inverse;
130  _ArrayIncrementIndex_(ndims,dim,index,done);
131  }
132  /* calculate the second source function */
133  IERR ShallowWater2DSourceFunction2(SourceC,u,x,solver,mpi,t,_YDIR_); CHECKERR(ierr);
134  /* calculate the left and right interface source terms */
135  IERR solver->InterpolateInterfacesHyp(SourceL,SourceC,u,x, 1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
136  IERR solver->InterpolateInterfacesHyp(SourceR,SourceC,u,x,-1,_YDIR_,solver,mpi,0); CHECKERR(ierr);
137  /* calculate the final interface source term */
138  IERR param->SourceUpwind(SourceI,SourceL,SourceR,u,_YDIR_,solver,t);
139  /* calculate the final cell-centered source term */
140  done = 0; _ArraySetValue_(index,ndims,0);
141  while (!done) {
142  _ArrayCopy1D_(index,index1,ndims);
143  _ArrayCopy1D_(index,index2,ndims); index2[_YDIR_]++;
144  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p );
145  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
146  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
147  double dy_inverse; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,dxinv,dy_inverse);
148  double h, uvel, vvel; _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);
149  double term[_MODEL_NVARS_] = { 0.0, 0.0, -param->g * (h + param->b[p]) };
150  for (v=0; v<_MODEL_NVARS_; v++)
151  source[_MODEL_NVARS_*p+v] += term[v]*(SourceI[_MODEL_NVARS_*p2+v]-SourceI[_MODEL_NVARS_*p1+v])*dy_inverse;
152  uvel = vvel = h; /* useless statement to avoid compiler warning */
153  _ArrayIncrementIndex_(ndims,dim,index,done);
154  }
155 
156  /* Naive addition of Coriolis force terms */
157  done = 0; _ArraySetValue_(index,ndims,0);
158  while (!done) {
159  _ArrayIndex1D_(ndims,dim,index ,ghosts,p);
160  double h, uvel, vvel;
161  _ShallowWater2DGetFlowVar_((u+_MODEL_NVARS_*p),h,uvel,vvel);
162  double ycoord; _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->x,ycoord);
163  double coeff = param->fhat + param->beta * (ycoord - 0.5*param->D);
164  double coriolis_x = coeff * vvel,
165  coriolis_y = -coeff * uvel;
166  source[_MODEL_NVARS_*p+1] += coriolis_x;
167  source[_MODEL_NVARS_*p+2] += coriolis_y;
168  _ArrayIncrementIndex_(ndims,dim,index,done);
169  }
170 
171  return(0);
172 }
_YDIR_
#define _YDIR_
Definition: euler2d.h:41
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
HyPar::fL
double * fL
Definition: hypar.h:139
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_GetCoordinate_
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
ShallowWater2DSourceFunction2
static int ShallowWater2DSourceFunction2(double *, double *, double *, void *, void *, double, int)
Definition: ShallowWater2DSource.c:251
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
HyPar::fluxI
double * fluxI
Definition: hypar.h:136
ShallowWater2D::SourceUpwind
int(* SourceUpwind)(double *, double *, double *, double *, int, void *, double)
Definition: shallowwater2d.h:256
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::fluxC
double * fluxC
Definition: hypar.h:128
HyPar::InterpolateInterfacesHyp
int(* InterpolateInterfacesHyp)(double *, double *, double *, double *, int, int, void *, void *, int)
Definition: hypar.h:224
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
_ShallowWater2DGetFlowVar_
#define _ShallowWater2DGetFlowVar_(u, h, uvel, vvel)
Definition: shallowwater2d.h:73
HyPar::dxinv
double * dxinv
Definition: hypar.h:110
ShallowWater2DSourceFunction1
static int ShallowWater2DSourceFunction1(double *, double *, double *, void *, void *, double, int)
Definition: ShallowWater2DSource.c:194
HyPar::ndims
int ndims
Definition: hypar.h:26
IERR
#define IERR
Definition: basic.h:16
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar::x
double * x
Definition: hypar.h:107
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.
Definition: shallowwater2d.h:240
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_XDIR_
#define _XDIR_
Definition: euler1d.h:75
HyPar::fR
double * fR
Definition: hypar.h:139