HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
ShallowWater1DUpwind.c
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <math.h>
8 #include <basic.h>
9 #include <arrayfunctions.h>
10 #include <mathfunctions.h>
11 #include <matmult_native.h>
12 #include <physicalmodels/shallowwater1d.h>
13 #include <hypar.h>
14 
31 int ShallowWater1DUpwindRoe(
32  double *fI,
33  double *fL,
34  double *fR,
35  double *uL,
36  double *uR,
37  double *u,
38  int dir,
39  void *s,
40  double t
41  )
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 }
95 
114 int ShallowWater1DUpwindLLF(
115  double *fI,
116  double *fL,
117  double *fR,
118  double *uL,
119  double *uR,
120  double *u,
121  int dir,
122  void *s,
123  double t
124  )
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 }
185 
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
max3
#define max3(a, b, c)
Definition: math_ops.h:27
absolute
#define absolute(a)
Definition: math_ops.h:32
shallowwater1d.h
1D Shallow Water Equations
mathfunctions.h
Contains function definitions for common mathematical functions.
_ShallowWater1DRightEigenvectors_
#define _ShallowWater1DRightEigenvectors_(u, R, p, dir)
Definition: shallowwater1d.h:117
basic.h
Some basic definitions and macros.
MatVecMult2
#define MatVecMult2(N, y, A, x)
Definition: matmult_native.h:128
HyPar::ndims
int ndims
Definition: hypar.h:26
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
hypar.h
Contains structure definition for hypar.
_ShallowWater1DRoeAverage_
#define _ShallowWater1DRoeAverage_(uavg, uL, uR, p)
Definition: shallowwater1d.h:74
_ShallowWater1DEigenvalues_
#define _ShallowWater1DEigenvalues_(u, D, p, dir)
Definition: shallowwater1d.h:93
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_ArraySetValue_
#define _ArraySetValue_(x, size, value)
Definition: arrayfunctions.h:169
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
MatMult2
#define MatMult2(N, A, X, Y)
Definition: matmult_native.h:32
_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
ShallowWater1DUpwindLLF
int ShallowWater1DUpwindLLF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: ShallowWater1DUpwind.c:114
matmult_native.h
Contains macros and function definitions for common matrix multiplication.
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
arrayfunctions.h
Contains macros and function definitions for common array operations.
_ShallowWater1DLeftEigenvectors_
#define _ShallowWater1DLeftEigenvectors_(u, L, p, dir)
Definition: shallowwater1d.h:106
ShallowWater1DUpwindRoe
int ShallowWater1DUpwindRoe(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: ShallowWater1DUpwind.c:31