HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
LinearADRUpwind.c
Go to the documentation of this file.
1 
7 #include <stdlib.h>
8 #include <math.h>
9 #include <basic.h>
10 #include <arrayfunctions.h>
11 #include <mathfunctions.h>
13 #include <hypar.h>
14 
16 int LinearADRUpwind( double *fI,
17  double *fL,
18  double *fR,
19  double *uL,
20  double *uR,
21  double *u,
22  int dir,
23  void *s,
24  double t
25  )
26 {
27  HyPar *solver = (HyPar*) s;
28  LinearADR *param = (LinearADR*) solver->physics;
29  int done,v;
30 
31  int ndims = solver->ndims;
32  int nvars = solver->nvars;
33  int ghosts= solver->ghosts;
34  int *dim = solver->dim_local;
35 
36  double *a = param->a;
37 
38  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
39  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
40  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
41 
42  if (param->constant_advection == 1) {
43 
44  done = 0; _ArraySetValue_(index_outer,ndims,0);
45  while (!done) {
46  _ArrayCopy1D_(index_outer,index_inter,ndims);
47  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
48  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
49  for (v = 0; v < nvars; v++) {
50  fI[nvars*p+v] = (a[nvars*dir+v] > 0 ? fL[nvars*p+v] : fR[nvars*p+v] );
51  }
52  }
53  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
54  }
55 
56  } else if (param->constant_advection == 0) {
57 
58  done = 0; _ArraySetValue_(index_outer,ndims,0);
59  while (!done) {
60  _ArrayCopy1D_(index_outer,index_inter,ndims);
61  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
62  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
63  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
64  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
65  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
66  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
67  for (v = 0; v < nvars; v++) {
68  double eigL = a[nvars*ndims*pL+nvars*dir+v],
69  eigR = a[nvars*ndims*pR+nvars*dir+v];
70  if ((eigL > 0) && (eigR > 0)) {
71  fI[nvars*p+v] = fL[nvars*p+v];
72  } else if ((eigL < 0) && (eigR < 0)) {
73  fI[nvars*p+v] = fR[nvars*p+v];
74  } else {
75  double alpha = max(absolute(eigL), absolute(eigR));
76  fI[nvars*p+v] = 0.5 * (fL[nvars*p+v] + fR[nvars*p+v] - alpha * (uR[nvars*p+v] - uL[nvars*p+v]));
77  }
78  }
79  }
80  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
81  }
82 
83  } else {
84 
85  done = 0; _ArraySetValue_(index_outer,ndims,0);
86  while (!done) {
87  _ArrayCopy1D_(index_outer,index_inter,ndims);
88  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
89  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
90  for (v = 0; v < nvars; v++) {
91  fI[nvars*p+v] = 0.0;
92  }
93  }
94  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
95  }
96 
97  }
98 
99  return(0);
100 }
101 
103 int LinearADRCenteredFlux( double *fI,
104  double *fL,
105  double *fR,
106  double *uL,
107  double *uR,
108  double *u,
109  int dir,
110  void *s,
111  double t
112  )
113 {
114  HyPar *solver = (HyPar*) s;
115  LinearADR *param = (LinearADR*) solver->physics;
116 
117  int ndims = solver->ndims;
118  int nvars = solver->nvars;
119  int ghosts= solver->ghosts;
120  int *dim = solver->dim_local;
121 
122  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
123  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
124  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
125 
126  int done = 0; _ArraySetValue_(index_outer,ndims,0);
127  while (!done) {
128  _ArrayCopy1D_(index_outer,index_inter,ndims);
129  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
130  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
131  for (int v = 0; v < nvars; v++) {
132  fI[nvars*p+v] = 0.5 * (fL[nvars*p+v] + fR[nvars*p+v]);
133  }
134  }
135  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
136  }
137 
138  return(0);
139 }
int constant_advection
Definition: linearadr.h:40
#define absolute(a)
Definition: math_ops.h:32
int nvars
Definition: hypar.h:29
Contains function definitions for common mathematical functions.
Some basic definitions and macros.
int ndims
Definition: hypar.h:26
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
double * a
Definition: linearadr.h:50
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Contains structure definition for hypar.
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Linear Advection-Diffusion-Reaction model.
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
int LinearADRCenteredFlux(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
#define _ArrayCopy1D_(x, y, size)
#define max(a, b)
Definition: math_ops.h:18
Contains macros and function definitions for common array operations.
int LinearADRUpwind(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)