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
ShallowWater2DUpwind.c File Reference

Contains functions to compute the upwind flux at grid interfaces for the 2D shallow water equations. More...

#include <stdlib.h>
#include <math.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <matmult_native.h>
#include <physicalmodels/shallowwater2d.h>
#include <hypar.h>

Go to the source code of this file.

Functions

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

Detailed Description

Contains functions to compute the upwind flux at grid interfaces for the 2D shallow water equations.

Author
Debojyoti Ghosh

Definition in file ShallowWater2DUpwind.c.

Function Documentation

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}

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
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension
sSolver object of type HyPar
tCurrent solution time

Definition at line 31 of file ShallowWater2DUpwind.c.

42 {
43  HyPar *solver = (HyPar*) s;
44  ShallowWater2D *param = (ShallowWater2D*) solver->physics;
45  int done,k;
46 
47  int ndims = solver->ndims;
48  int ghosts= solver->ghosts;
49  int *dim = solver->dim_local;
50 
51  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
52  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
53  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
54  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_],
55  DL[_MODEL_NVARS_*_MODEL_NVARS_], modA[_MODEL_NVARS_*_MODEL_NVARS_];
56 
57  done = 0; _ArraySetValue_(index_outer,ndims,0);
58  while (!done) {
59  _ArrayCopy1D_(index_outer,index_inter,ndims);
60  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
61  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
62  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
63  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
64  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
65  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
66  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
67 
68  /* Roe's upwinding scheme */
69 
70  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
71  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
72  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
73 
74  _ShallowWater2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
75  _ShallowWater2DEigenvalues_ (uavg,D,param,dir);
76  _ShallowWater2DLeftEigenvectors_ (uavg,L,param,dir);
77  _ShallowWater2DRightEigenvectors_ (uavg,R,param,dir);
78 
79  k = 0; D[k] = absolute(D[k]);
80  k = 4; D[k] = absolute(D[k]);
81  k = 8; D[k] = absolute(D[k]);
82 
83  MatMult3(_MODEL_NVARS_,DL,D,L);
84  MatMult3(_MODEL_NVARS_,modA,R,DL);
85 
86  udiss[0] = modA[0*_MODEL_NVARS_+0]*udiff[0] + modA[0*_MODEL_NVARS_+1]*udiff[1] + modA[0*_MODEL_NVARS_+2]*udiff[2];
87  udiss[1] = modA[1*_MODEL_NVARS_+0]*udiff[0] + modA[1*_MODEL_NVARS_+1]*udiff[1] + modA[1*_MODEL_NVARS_+2]*udiff[2];
88  udiss[2] = modA[2*_MODEL_NVARS_+0]*udiff[0] + modA[2*_MODEL_NVARS_+1]*udiff[1] + modA[2*_MODEL_NVARS_+2]*udiff[2];
89 
90  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
91  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
92  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
93  }
94  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
95  }
96 
97  return(0);
98 }
_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
MatMult3
#define MatMult3(N, A, X, Y)
Definition: matmult_native.h:43
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ShallowWater2DLeftEigenvectors_
#define _ShallowWater2DLeftEigenvectors_(u, L, p, dir)
Definition: shallowwater2d.h:154
_ShallowWater2DEigenvalues_
#define _ShallowWater2DEigenvalues_(u, D, p, dir)
Definition: shallowwater2d.h:129
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ShallowWater2DRightEigenvectors_
#define _ShallowWater2DRightEigenvectors_(u, R, p, dir)
Definition: shallowwater2d.h:197
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
_ShallowWater2DRoeAverage_
#define _ShallowWater2DRoeAverage_(uavg, uL, uR, p)
Definition: shallowwater2d.h:108
HyPar::ndims
int ndims
Definition: hypar.h:26
absolute
#define absolute(a)
Definition: math_ops.h:32
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
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.

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
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension
sSolver object of type HyPar
tCurrent solution time

Definition at line 118 of file ShallowWater2DUpwind.c.

129 {
130  HyPar *solver = (HyPar*) s;
131  ShallowWater2D *param = (ShallowWater2D*) solver->physics;
132  int done,k;
133 
134  int ndims = solver->ndims;
135  int *dim = solver->dim_local;
136  int ghosts = solver->ghosts;
137 
138  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
139  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
140  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
141  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
142 
143  done = 0; _ArraySetValue_(index_outer,ndims,0);
144  while (!done) {
145  _ArrayCopy1D_(index_outer,index_inter,ndims);
146  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
147  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
148  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
149  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
150  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
151  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
152  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
153  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
154 
155  /* Local Lax-Friedrich upwinding scheme */
156 
157  _ShallowWater2DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
158  _ShallowWater2DEigenvalues_ (uavg,D,param,dir);
159  _ShallowWater2DLeftEigenvectors_ (uavg,L,param,dir);
160  _ShallowWater2DRightEigenvectors_(uavg,R,param,dir);
161 
162  /* calculate characteristic fluxes and variables */
163  MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
164  MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
165  MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
166  MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
167 
168  for (k = 0; k < _MODEL_NVARS_; k++) {
169  double eigL,eigC,eigR;
170  _ShallowWater2DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,dir);
171  eigL = D[k*_MODEL_NVARS_+k];
172  _ShallowWater2DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,dir);
173  eigR = D[k*_MODEL_NVARS_+k];
174  _ShallowWater2DEigenvalues_(uavg,D,param,dir);
175  eigC = D[k*_MODEL_NVARS_+k];
176 
177  double alpha = max3(absolute(eigL),absolute(eigC),absolute(eigR));
178  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
179  }
180 
181  /* calculate the interface flux from the characteristic flux */
182  MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
183  }
184  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
185  }
186 
187  return(0);
188 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
max3
#define max3(a, b, c)
Definition: math_ops.h:27
_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
_ShallowWater2DLeftEigenvectors_
#define _ShallowWater2DLeftEigenvectors_(u, L, p, dir)
Definition: shallowwater2d.h:154
_ShallowWater2DEigenvalues_
#define _ShallowWater2DEigenvalues_(u, D, p, dir)
Definition: shallowwater2d.h:129
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ShallowWater2DRightEigenvectors_
#define _ShallowWater2DRightEigenvectors_(u, R, p, dir)
Definition: shallowwater2d.h:197
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
_ShallowWater2DRoeAverage_
#define _ShallowWater2DRoeAverage_(uavg, uL, uR, p)
Definition: shallowwater2d.h:108
HyPar::ndims
int ndims
Definition: hypar.h:26
absolute
#define absolute(a)
Definition: math_ops.h:32
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
MatVecMult3
#define MatVecMult3(N, y, A, x)
Definition: matmult_native.h:137