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

Contains functions to compute the upwind flux at grid interfaces for the 1D 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/shallowwater1d.h>
#include <hypar.h>

Go to the source code of this file.

Functions

int ShallowWater1DUpwindRoe (double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
 
int ShallowWater1DUpwindLLF (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 1D shallow water equations.

Author
Debojyoti Ghosh

Definition in file ShallowWater1DUpwind.c.

Function Documentation

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}

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 (unused since this is a 1D system)
sSolver object of type HyPar
tCurrent solution time

Definition at line 31 of file ShallowWater1DUpwind.c.

42 {
43  HyPar *solver = (HyPar*) s;
44  ShallowWater1D *param = (ShallowWater1D*) 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 
73  _ShallowWater1DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
74  _ShallowWater1DEigenvalues_ (uavg,D,param,0);
75  _ShallowWater1DLeftEigenvectors_ (uavg,L,param,0);
76  _ShallowWater1DRightEigenvectors_ (uavg,R,param,0);
77 
78  k = 0; D[k] = absolute(D[k]);
79  k = 3; D[k] = absolute(D[k]);
80 
81  MatMult2(2,DL,D,L);
82  MatMult2(2,modA,R,DL);
83 
84  udiss[0] = modA[0*_MODEL_NVARS_+0]*udiff[0] + modA[0*_MODEL_NVARS_+1]*udiff[1];
85  udiss[1] = modA[1*_MODEL_NVARS_+0]*udiff[0] + modA[1*_MODEL_NVARS_+1]*udiff[1];
86 
87  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
88  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
89  }
90  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
91  }
92 
93  return(0);
94 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ShallowWater1DRoeAverage_
#define _ShallowWater1DRoeAverage_(uavg, uL, uR, p)
Definition: shallowwater1d.h:74
ShallowWater1D
Structure containing variables and parameters specific to the 1D Shallow Water equations. This structure contains the physical parameters, variables, and function pointers specific to the 1D ShallowWater equations.
Definition: shallowwater1d.h:137
_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
MatMult2
#define MatMult2(N, A, X, Y)
Definition: matmult_native.h:32
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ShallowWater1DLeftEigenvectors_
#define _ShallowWater1DLeftEigenvectors_(u, L, p, dir)
Definition: shallowwater1d.h:106
_ShallowWater1DEigenvalues_
#define _ShallowWater1DEigenvalues_(u, D, p, dir)
Definition: shallowwater1d.h:93
_ShallowWater1DRightEigenvectors_
#define _ShallowWater1DRightEigenvectors_(u, R, p, dir)
Definition: shallowwater1d.h:117
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::ndims
int ndims
Definition: hypar.h:26
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
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.

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 (unused since this is a 1D system)
sSolver object of type HyPar
tCurrent solution time

Definition at line 114 of file ShallowWater1DUpwind.c.

125 {
126  HyPar *solver = (HyPar*) s;
127  ShallowWater1D *param = (ShallowWater1D*) solver->physics;
128  int done,k;
129 
130  int ndims = solver->ndims;
131  int *dim = solver->dim_local;
132  int ghosts = solver->ghosts;
133 
134  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
135  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
136  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
137  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
138 
139  done = 0; _ArraySetValue_(index_outer,ndims,0);
140  while (!done) {
141  _ArrayCopy1D_(index_outer,index_inter,ndims);
142  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
143  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
144  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
145  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
146  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
147  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
148  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
149  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
150 
151  /* Local Lax-Friedrich upwinding scheme */
152 
153  _ShallowWater1DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
154  _ShallowWater1DEigenvalues_ (uavg,D,param,0);
155  _ShallowWater1DLeftEigenvectors_ (uavg,L,param,0);
156  _ShallowWater1DRightEigenvectors_(uavg,R,param,0);
157 
158  /* calculate characteristic fluxes and variables */
159  MatVecMult2(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
160  MatVecMult2(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
161  MatVecMult2(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
162  MatVecMult2(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
163 
164  for (k = 0; k < _MODEL_NVARS_; k++) {
165  double eigL,eigC,eigR;
166  _ShallowWater1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
167  eigL = D[k*_MODEL_NVARS_+k];
168  _ShallowWater1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
169  eigR = D[k*_MODEL_NVARS_+k];
170  _ShallowWater1DEigenvalues_(uavg,D,param,0);
171  eigC = D[k*_MODEL_NVARS_+k];
172 
173  double alpha = max3(absolute(eigL),absolute(eigC),absolute(eigR));
174  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
175  }
176 
177  /* calculate the interface flux from the characteristic flux */
178  MatVecMult2(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
179  }
180  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
181  }
182 
183  return(0);
184 }
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
_ShallowWater1DRoeAverage_
#define _ShallowWater1DRoeAverage_(uavg, uL, uR, p)
Definition: shallowwater1d.h:74
max3
#define max3(a, b, c)
Definition: math_ops.h:27
ShallowWater1D
Structure containing variables and parameters specific to the 1D Shallow Water equations. This structure contains the physical parameters, variables, and function pointers specific to the 1D ShallowWater equations.
Definition: shallowwater1d.h:137
_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
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ShallowWater1DLeftEigenvectors_
#define _ShallowWater1DLeftEigenvectors_(u, L, p, dir)
Definition: shallowwater1d.h:106
_ShallowWater1DEigenvalues_
#define _ShallowWater1DEigenvalues_(u, D, p, dir)
Definition: shallowwater1d.h:93
_ShallowWater1DRightEigenvectors_
#define _ShallowWater1DRightEigenvectors_(u, R, p, dir)
Definition: shallowwater1d.h:117
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
MatVecMult2
#define MatVecMult2(N, y, A, x)
Definition: matmult_native.h:128
HyPar::ndims
int ndims
Definition: hypar.h:26
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23