HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
ShallowWater1DUpwindRoe
int ShallowWater1DUpwindRoe(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: ShallowWater1DUpwind.c:31
_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
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
ShallowWater1DUpwindLLF
int ShallowWater1DUpwindLLF(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: ShallowWater1DUpwind.c:114
mathfunctions.h
Contains function definitions for common mathematical functions.
_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.h
Contains structure definition for hypar.
shallowwater1d.h
1D Shallow Water Equations
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
arrayfunctions.h
Contains macros and function definitions for common array operations.
absolute
#define absolute(a)
Definition: math_ops.h:32
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
matmult_native.h
Contains macros and function definitions for common matrix multiplication.