HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Euler1DUpwind.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/euler1d.h>
13 #include <hypar.h>
14 
30 int Euler1DUpwindRoe(
31  double *fI,
32  double *fL,
33  double *fR,
34  double *uL,
35  double *uR,
36  double *u,
37  int dir,
38  void *s,
39  double t
40  )
41 {
42  HyPar *solver = (HyPar*) s;
43  Euler1D *param = (Euler1D*) solver->physics;
44  int done,k;
45  _DECLARE_IERR_;
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  _Euler1DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
75  _Euler1DEigenvalues_ (uavg,D,param,0);
76  _Euler1DLeftEigenvectors_ (uavg,L,param,0);
77  _Euler1DRightEigenvectors_ (uavg,R,param,0);
78 
79  double kappa = max(param->grav_field[pL],param->grav_field[pR]);
80  k = 0; D[k] = kappa*absolute(D[k]);
81  k = 4; D[k] = kappa*absolute(D[k]);
82  k = 8; D[k] = kappa*absolute(D[k]);
83 
84  MatMult3(3,DL,D,L);
85  MatMult3(3,modA,R,DL);
86 
87  udiss[0] = modA[0*_MODEL_NVARS_+0]*udiff[0] + modA[0*_MODEL_NVARS_+1]*udiff[1] + modA[0*_MODEL_NVARS_+2]*udiff[2];
88  udiss[1] = modA[1*_MODEL_NVARS_+0]*udiff[0] + modA[1*_MODEL_NVARS_+1]*udiff[1] + modA[1*_MODEL_NVARS_+2]*udiff[2];
89  udiss[2] = modA[2*_MODEL_NVARS_+0]*udiff[0] + modA[2*_MODEL_NVARS_+1]*udiff[1] + modA[2*_MODEL_NVARS_+2]*udiff[2];
90 
91  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
92  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
93  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
94  }
95  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
96  }
97 
98  return(0);
99 }
100 
115 int Euler1DUpwindRF(
116  double *fI,
117  double *fL,
118  double *fR,
119  double *uL,
120  double *uR,
121  double *u,
122  int dir,
123  void *s,
124  double t
125  )
126 {
127  HyPar *solver = (HyPar*) s;
128  Euler1D *param = (Euler1D*) solver->physics;
129  int done,k;
130  _DECLARE_IERR_;
131 
132  int ndims = solver->ndims;
133  int *dim = solver->dim_local;
134  int ghosts = solver->ghosts;
135 
136  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
137  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
138  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
139  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
140 
141  done = 0; _ArraySetValue_(index_outer,ndims,0);
142  while (!done) {
143  _ArrayCopy1D_(index_outer,index_inter,ndims);
144  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
145  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
146  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
147  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
148  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
149  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
150  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
151  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
152  double kappa = max(param->grav_field[pL],param->grav_field[pR]);
153 
154  /* Roe-Fixed upwinding scheme */
155 
156  _Euler1DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
157  _Euler1DEigenvalues_ (uavg,D,param,0);
158  _Euler1DLeftEigenvectors_ (uavg,L,param,0);
159  _Euler1DRightEigenvectors_(uavg,R,param,0);
160 
161  /* calculate characteristic fluxes and variables */
162  MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
163  MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
164  MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
165  MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
166 
167  for (k = 0; k < _MODEL_NVARS_; k++) {
168  double eigL,eigC,eigR;
169  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
170  eigL = D[k*_MODEL_NVARS_+k];
171  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
172  eigR = D[k*_MODEL_NVARS_+k];
173  _Euler1DEigenvalues_(uavg,D,param,0);
174  eigC = D[k*_MODEL_NVARS_+k];
175 
176  if ((eigL > 0) && (eigC > 0) && (eigR > 0)) fc[k] = fcL[k];
177  else if ((eigL < 0) && (eigC < 0) && (eigR < 0)) fc[k] = fcR[k];
178  else {
179  double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
180  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
181  }
182 
183  }
184 
185  /* calculate the interface flux from the characteristic flux */
186  MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
187  }
188  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
189  }
190 
191  return(0);
192 }
193 
212 int Euler1DUpwindLLF(
213  double *fI,
214  double *fL,
215  double *fR,
216  double *uL,
217  double *uR,
218  double *u,
219  int dir,
220  void *s,
221  double t
222  )
223 {
224  HyPar *solver = (HyPar*) s;
225  Euler1D *param = (Euler1D*) solver->physics;
226  int done,k;
227  _DECLARE_IERR_;
228 
229  int ndims = solver->ndims;
230  int *dim = solver->dim_local;
231  int ghosts = solver->ghosts;
232 
233  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
234  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
235  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
236  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
237 
238  done = 0; _ArraySetValue_(index_outer,ndims,0);
239  while (!done) {
240  _ArrayCopy1D_(index_outer,index_inter,ndims);
241  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
242  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
243  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
244  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
245  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
246  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
247  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
248  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
249  double kappa = max(param->grav_field[pL],param->grav_field[pR]);
250 
251  /* Local Lax-Friedrich upwinding scheme */
252 
253  _Euler1DRoeAverage_ (uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
254  _Euler1DEigenvalues_ (uavg,D,param,0);
255  _Euler1DLeftEigenvectors_ (uavg,L,param,0);
256  _Euler1DRightEigenvectors_(uavg,R,param,0);
257 
258  /* calculate characteristic fluxes and variables */
259  MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
260  MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
261  MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
262  MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
263 
264  for (k = 0; k < _MODEL_NVARS_; k++) {
265  double eigL,eigC,eigR;
266  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
267  eigL = D[k*_MODEL_NVARS_+k];
268  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
269  eigR = D[k*_MODEL_NVARS_+k];
270  _Euler1DEigenvalues_(uavg,D,param,0);
271  eigC = D[k*_MODEL_NVARS_+k];
272 
273  double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
274  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
275  }
276 
277  /* calculate the interface flux from the characteristic flux */
278  MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
279  }
280  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
281  }
282 
283  return(0);
284 }
285 
293 int Euler1DUpwindSWFS(
294  double *fI,
295  double *fL,
296  double *fR,
297  double *uL,
298  double *uR,
299  double *u,
300  int dir,
301  void *s,
302  double t
303  )
304 {
305  HyPar *solver = (HyPar*) s;
306  Euler1D *param = (Euler1D*) solver->physics;
307  int done;
308  _DECLARE_IERR_;
309 
310  int ndims = solver->ndims;
311  int *dim = solver->dim_local;
312 
313  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
314  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
315  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
316  static double fp[_MODEL_NVARS_], fm[_MODEL_NVARS_],uavg[_MODEL_NVARS_];
317 
318  done = 0; _ArraySetValue_(index_outer,ndims,0);
319  while (!done) {
320  _ArrayCopy1D_(index_outer,index_inter,ndims);
321  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
322  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
323  double rho,v,e,P,c,gamma=param->gamma,term,Mach;
324 
325  /* Steger Warming flux splitting */
326  _Euler1DRoeAverage_(uavg,(uL+_MODEL_NVARS_*p),(uR+_MODEL_NVARS_*p),param);
327  _Euler1DGetFlowVar_(uavg,rho,v,e,P,param);
328  Mach = v/sqrt(gamma*P/rho);
329 
330  if (Mach < -1.0) {
331 
332  _ArrayCopy1D3_((fR+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);
333 
334  } else if (Mach < 1.0) {
335 
336  _Euler1DGetFlowVar_((uL+_MODEL_NVARS_*p),rho,v,e,P,param);
337  c = sqrt(gamma*P/rho);
338  term = rho/(2.0*gamma);
339 
340  fp[0] = term * (2*gamma*v + c - v);
341  fp[1] = term * (2*(gamma-1.0)*v*v + (v+c)*(v+c));
342  fp[2] = term * ((gamma-1.0)*v*v*v + 0.5*(v+c)*(v+c)*(v+c) + ((3.0-gamma)*(v+c)*c*c)/(2.0*(gamma-1.0)));
343 
344  _Euler1DGetFlowVar_((uR+_MODEL_NVARS_*p),rho,v,e,P,param);
345  c = sqrt(gamma*P/rho);
346  term = rho/(2.0*gamma);
347 
348  fm[0] = term * (v - c);
349  fm[1] = term * (v-c) * (v-c);
350  fm[2] = term * (0.5*(v-c)*(v-c)*(v-c) + ((3.0-gamma)*(v-c)*c*c)/(2.0*(gamma-1.0)));
351 
352  _ArrayAdd1D_((fI+_MODEL_NVARS_*p),fp,fm,_MODEL_NVARS_);
353 
354  } else {
355 
356  _ArrayCopy1D3_((fL+_MODEL_NVARS_*p),(fI+_MODEL_NVARS_*p),_MODEL_NVARS_);
357 
358  }
359 
360  }
361  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
362  }
363 
364  return(0);
365 }
366 
376 int Euler1DUpwindRusanov(
377  double *fI,
378  double *fL,
379  double *fR,
380  double *uL,
381  double *uR,
382  double *u,
383  int dir,
384  void *s,
385  double t
386  )
387 {
388  HyPar *solver = (HyPar*) s;
389  Euler1D *param = (Euler1D*) solver->physics;
390  int done,k;
391  _DECLARE_IERR_;
392 
393  int ndims = solver->ndims;
394  int ghosts= solver->ghosts;
395  int *dim = solver->dim_local;
396 
397  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
398  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
399  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
400 
401  static double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_];
402 
403  done = 0; _ArraySetValue_(index_outer,ndims,0);
404  while (!done) {
405 
406  _ArrayCopy1D_(index_outer,index_inter,ndims);
407 
408  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
409 
410  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
411  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
412  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
413  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
414  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
415 
416  _Euler1DRoeAverage_(uavg,(u+_MODEL_NVARS_*pL),(u+_MODEL_NVARS_*pR),param);
417  for (k = 0; k < _MODEL_NVARS_; k++) udiff[k] = 0.5 * (uR[_MODEL_NVARS_*p+k] - uL[_MODEL_NVARS_*p+k]);
418 
419  double rho, uvel, E, P, c;
420 
421  _Euler1DGetFlowVar_((u+_MODEL_NVARS_*pL),rho,uvel,E,P,param);
422  c = param->gamma*P/rho;
423  double alphaL = c + absolute(uvel);
424 
425  _Euler1DGetFlowVar_((u+_MODEL_NVARS_*pR),rho,uvel,E,P,param);
426  c = param->gamma*P/rho;
427  double alphaR = c + absolute(uvel);
428 
429  _Euler1DGetFlowVar_(uavg,rho,uvel,E,P,param);
430  c = param->gamma*P/rho;
431  double alphaavg = c + absolute(uvel);
432 
433  double kappa = max(param->grav_field[pL],param->grav_field[pR]);
434  double alpha = kappa*max3(alphaL,alphaR,alphaavg);
435 
436  for (k = 0; k < _MODEL_NVARS_; k++) {
437  fI[_MODEL_NVARS_*p+k] = 0.5 * (fL[_MODEL_NVARS_*p+k]+fR[_MODEL_NVARS_*p+k]) - alpha*udiff[k];
438  }
439 
440  }
441 
442  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
443 
444  }
445 
446  return(0);
447 }
448 
453 int Euler1DUpwinddFRoe(
454  double *fI,
455  double *fL,
456  double *fR,
457  double *uL,
458  double *uR,
459  double *u,
460  int dir,
461  void *s,
462  double t
463  )
464 {
465  HyPar *solver = (HyPar*) s;
466  Euler1D *param = (Euler1D*) solver->physics;
467  int done,k;
468  _DECLARE_IERR_;
469 
470  int ndims = solver->ndims;
471  int ghosts= solver->ghosts;
472  int *dim = solver->dim_local;
473  double *uref = param->solution;
474 
475  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
476  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
477  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
478  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_],
479  DL[_MODEL_NVARS_*_MODEL_NVARS_], modA[_MODEL_NVARS_*_MODEL_NVARS_];
480 
481  done = 0; _ArraySetValue_(index_outer,ndims,0);
482  while (!done) {
483  _ArrayCopy1D_(index_outer,index_inter,ndims);
484  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
485  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
486  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
487  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
488  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
489  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
490  double udiff[_MODEL_NVARS_], uavg[_MODEL_NVARS_],udiss[_MODEL_NVARS_];
491 
492  /* Roe's upwinding scheme */
493 
494  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
495  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
496  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
497 
498  _Euler1DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param);
499  _Euler1DEigenvalues_ (uavg,D,param,0);
500  _Euler1DLeftEigenvectors_ (uavg,L,param,0);
501  _Euler1DRightEigenvectors_ (uavg,R,param,0);
502 
503  double kappa = max(param->grav_field[pL],param->grav_field[pR]);
504  k = 0; D[k] = 0.0;
505  k = 4; D[k] = kappa*absolute(D[k]);
506  k = 8; D[k] = kappa*absolute(D[k]);
507 
508  MatMult3(3,DL,D,L);
509  MatMult3(3,modA,R,DL);
510 
511  udiss[0] = modA[0*_MODEL_NVARS_+0]*udiff[0] + modA[0*_MODEL_NVARS_+1]*udiff[1] + modA[0*_MODEL_NVARS_+2]*udiff[2];
512  udiss[1] = modA[1*_MODEL_NVARS_+0]*udiff[0] + modA[1*_MODEL_NVARS_+1]*udiff[1] + modA[1*_MODEL_NVARS_+2]*udiff[2];
513  udiss[2] = modA[2*_MODEL_NVARS_+0]*udiff[0] + modA[2*_MODEL_NVARS_+1]*udiff[1] + modA[2*_MODEL_NVARS_+2]*udiff[2];
514 
515  fI[_MODEL_NVARS_*p+0] = 0.5 * (fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0]) - udiss[0];
516  fI[_MODEL_NVARS_*p+1] = 0.5 * (fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1]) - udiss[1];
517  fI[_MODEL_NVARS_*p+2] = 0.5 * (fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2]) - udiss[2];
518  }
519  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
520  }
521 
522  return(0);
523 }
524 
534 int Euler1DUpwinddFRF(
535  double *fI,
536  double *fL,
537  double *fR,
538  double *uL,
539  double *uR,
540  double *u,
541  int dir,
542  void *s,
543  double t
544  )
545 {
546  HyPar *solver = (HyPar*) s;
547  Euler1D *param = (Euler1D*) solver->physics;
548  int done,k;
549  _DECLARE_IERR_;
550 
551  int ndims = solver->ndims;
552  int *dim = solver->dim_local;
553  int ghosts = solver->ghosts;
554  double *uref = param->solution;
555 
556  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
557  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
558  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
559  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
560 
561  done = 0; _ArraySetValue_(index_outer,ndims,0);
562  while (!done) {
563  _ArrayCopy1D_(index_outer,index_inter,ndims);
564  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
565  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
566  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
567  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
568  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
569  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
570  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
571  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
572  double kappa = max(param->grav_field[pL],param->grav_field[pR]);
573 
574  /* Roe-Fixed upwinding scheme */
575 
576  _Euler1DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param);
577  _Euler1DEigenvalues_ (uavg,D,param,0);
578  _Euler1DLeftEigenvectors_ (uavg,L,param,0);
579  _Euler1DRightEigenvectors_(uavg,R,param,0);
580 
581  /* calculate characteristic fluxes and variables */
582  MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
583  MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
584  MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
585  MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
586 
587  for (k = 0; k < _MODEL_NVARS_; k++) {
588  double eigL,eigC,eigR;
589  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
590  eigL = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
591  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
592  eigR = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
593  _Euler1DEigenvalues_(uavg,D,param,0);
594  eigC = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
595 
596  if ((eigL > 0) && (eigC > 0) && (eigR > 0)) fc[k] = fcL[k];
597  else if ((eigL < 0) && (eigC < 0) && (eigR < 0)) fc[k] = fcR[k];
598  else {
599  double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
600  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
601  }
602 
603  }
604 
605  /* calculate the interface flux from the characteristic flux */
606  MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
607  }
608  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
609  }
610 
611  return(0);
612 }
613 
623 int Euler1DUpwinddFLLF(
624  double *fI,
625  double *fL,
626  double *fR,
627  double *uL,
628  double *uR,
629  double *u,
630  int dir,
631  void *s,
632  double t
633  )
634 {
635  HyPar *solver = (HyPar*) s;
636  Euler1D *param = (Euler1D*) solver->physics;
637  int done,k;
638  _DECLARE_IERR_;
639 
640  int ndims = solver->ndims;
641  int *dim = solver->dim_local;
642  int ghosts = solver->ghosts;
643  double *uref = param->solution;
644 
645  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
646  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
647  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
648  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], L[_MODEL_NVARS_*_MODEL_NVARS_];
649 
650  done = 0; _ArraySetValue_(index_outer,ndims,0);
651  while (!done) {
652  _ArrayCopy1D_(index_outer,index_inter,ndims);
653  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
654  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
655  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
656  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
657  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
658  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
659  double uavg[_MODEL_NVARS_], fcL[_MODEL_NVARS_], fcR[_MODEL_NVARS_],
660  ucL[_MODEL_NVARS_], ucR[_MODEL_NVARS_], fc[_MODEL_NVARS_];
661  double kappa = max(param->grav_field[pL],param->grav_field[pR]);
662 
663  /* Local Lax-Friedrich upwinding scheme */
664 
665  _Euler1DRoeAverage_ (uavg,(uref+_MODEL_NVARS_*pL),(uref+_MODEL_NVARS_*pR),param);
666  _Euler1DEigenvalues_ (uavg,D,param,0);
667  _Euler1DLeftEigenvectors_ (uavg,L,param,0);
668  _Euler1DRightEigenvectors_(uavg,R,param,0);
669 
670  /* calculate characteristic fluxes and variables */
671  MatVecMult3(_MODEL_NVARS_,ucL,L,(uL+_MODEL_NVARS_*p));
672  MatVecMult3(_MODEL_NVARS_,ucR,L,(uR+_MODEL_NVARS_*p));
673  MatVecMult3(_MODEL_NVARS_,fcL,L,(fL+_MODEL_NVARS_*p));
674  MatVecMult3(_MODEL_NVARS_,fcR,L,(fR+_MODEL_NVARS_*p));
675 
676  for (k = 0; k < _MODEL_NVARS_; k++) {
677  double eigL,eigC,eigR;
678  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pL),D,param,0);
679  eigL = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
680  _Euler1DEigenvalues_((u+_MODEL_NVARS_*pR),D,param,0);
681  eigR = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
682  _Euler1DEigenvalues_(uavg,D,param,0);
683  eigC = (k == 0? 0.0 : D[k*_MODEL_NVARS_+k]);
684 
685  double alpha = kappa * max3(absolute(eigL),absolute(eigC),absolute(eigR));
686  fc[k] = 0.5 * (fcL[k] + fcR[k] + alpha * (ucL[k]-ucR[k]));
687  }
688 
689  /* calculate the interface flux from the characteristic flux */
690  MatVecMult3(_MODEL_NVARS_,(fI+_MODEL_NVARS_*p),R,fc);
691  }
692  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
693  }
694 
695  return(0);
696 }
max3
#define max3(a, b, c)
Definition: math_ops.h:27
Euler1DUpwindLLF
int Euler1DUpwindLLF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:212
Euler1D::gamma
double gamma
Definition: euler1d.h:274
absolute
#define absolute(a)
Definition: math_ops.h:32
mathfunctions.h
Contains function definitions for common mathematical functions.
Euler1D
Structure containing variables and parameters specific to the 1D Euler equations. This structure cont...
Definition: euler1d.h:273
Euler1DUpwinddFLLF
int Euler1DUpwinddFLLF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:623
Euler1DUpwindRusanov
int Euler1DUpwindRusanov(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:376
Euler1DUpwinddFRoe
int Euler1DUpwinddFRoe(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:453
basic.h
Some basic definitions and macros.
HyPar::ndims
int ndims
Definition: hypar.h:26
_Euler1DLeftEigenvectors_
#define _Euler1DLeftEigenvectors_(u, L, p, dir)
Definition: euler1d.h:225
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
_ArrayAdd1D_
#define _ArrayAdd1D_(x, a, b, size)
Definition: arrayfunctions.h:212
hypar.h
Contains structure definition for hypar.
_ArrayCopy1D3_
#define _ArrayCopy1D3_(x, y, size)
Definition: arrayfunctions.h:343
MatVecMult3
#define MatVecMult3(N, y, A, x)
Definition: matmult_native.h:137
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
Euler1DUpwindRF
int Euler1DUpwindRF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:115
_Euler1DEigenvalues_
#define _Euler1DEigenvalues_(u, D, p, dir)
Definition: euler1d.h:207
_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
Euler1DUpwindSWFS
int Euler1DUpwindSWFS(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:293
HyPar::physics
void * physics
Definition: hypar.h:266
_Euler1DRightEigenvectors_
#define _Euler1DRightEigenvectors_(u, R, p, dir)
Definition: euler1d.h:249
euler1d.h
1D Euler Equations (inviscid, compressible flows)
Euler1D::solution
double * solution
Definition: euler1d.h:290
Euler1DUpwindRoe
int Euler1DUpwindRoe(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:30
HyPar::ghosts
int ghosts
Definition: hypar.h:52
matmult_native.h
Contains macros and function definitions for common matrix multiplication.
MatMult3
#define MatMult3(N, A, X, Y)
Definition: matmult_native.h:43
_MODEL_NVARS_
#define _MODEL_NVARS_
Definition: euler1d.h:58
_Euler1DGetFlowVar_
#define _Euler1DGetFlowVar_(u, rho, v, e, P, p)
Definition: euler1d.h:81
_DECLARE_IERR_
#define _DECLARE_IERR_
Definition: basic.h:17
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
max
#define max(a, b)
Definition: math_ops.h:18
arrayfunctions.h
Contains macros and function definitions for common array operations.
_Euler1DRoeAverage_
#define _Euler1DRoeAverage_(uavg, uL, uR, p)
Definition: euler1d.h:169
Euler1D::grav_field
double * grav_field
Definition: euler1d.h:288
Euler1DUpwinddFRF
int Euler1DUpwinddFRF(double *fI, double *fL, double *fR, double *uL, double *uR, double *u, int dir, void *s, double t)
Definition: Euler1DUpwind.c:534