HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
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

◆ ShallowWater2DUpwindRoe()

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 }
absolute
#define absolute(a)
Definition: math_ops.h:32
_ShallowWater2DRoeAverage_
#define _ShallowWater2DRoeAverage_(uavg, uL, uR, p)
Definition: shallowwater2d.h:108
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
_ShallowWater2DLeftEigenvectors_
#define _ShallowWater2DLeftEigenvectors_(u, L, p, dir)
Definition: shallowwater2d.h:154
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ShallowWater2DEigenvalues_
#define _ShallowWater2DEigenvalues_(u, D, p, dir)
Definition: shallowwater2d.h:129
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::ghosts
int ghosts
Definition: hypar.h:52
MatMult3
#define MatMult3(N, A, X, Y)
Definition: matmult_native.h:43
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ShallowWater2DRightEigenvectors_
#define _ShallowWater2DRightEigenvectors_(u, R, p, dir)
Definition: shallowwater2d.h:197
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319

◆ ShallowWater2DUpwindLLF()

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 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
_ShallowWater2DRoeAverage_
#define _ShallowWater2DRoeAverage_(uavg, uL, uR, p)
Definition: shallowwater2d.h:108
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
_ShallowWater2DLeftEigenvectors_
#define _ShallowWater2DLeftEigenvectors_(u, L, p, dir)
Definition: shallowwater2d.h:154
MatVecMult3
#define MatVecMult3(N, y, A, x)
Definition: matmult_native.h:137
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ShallowWater2DEigenvalues_
#define _ShallowWater2DEigenvalues_(u, D, p, dir)
Definition: shallowwater2d.h:129
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
_ArrayIncrementIndex_
#define _ArrayIncrementIndex_(N, imax, i, done)
Definition: arrayfunctions.h:126
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ShallowWater2DRightEigenvectors_
#define _ShallowWater2DRightEigenvectors_(u, R, p, dir)
Definition: shallowwater2d.h:197
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319