HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Interp1PrimFifthOrderWENO.c
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <time.h>
8 
9 #include <basic.h>
10 #include <arrayfunctions.h>
11 #include <mathfunctions.h>
12 #include <interpolation.h>
13 #include <mpivars.h>
14 #include <hypar.h>
15 
16 #ifdef with_omp
17 #include <omp.h>
18 #endif
19 
20 #undef _MINIMUM_GHOSTS_
21 
25 #define _MINIMUM_GHOSTS_ 3
26 
75  double *fI,
76  double *fC,
77  double *u,
78  double *x,
79  int upw,
80  int dir,
81  void *s,
82  void *m,
83  int uflag
84  )
85 {
86  HyPar *solver = (HyPar*) s;
87  WENOParameters *weno = (WENOParameters*) solver->interp;
88 
89  int ghosts = solver->ghosts;
90  int ndims = solver->ndims;
91  int nvars = solver->nvars;
92  int *dim = solver->dim_local;
93  int *stride= solver->stride_with_ghosts;
94 
95  /* define some constants */
96  static const double one_sixth = 1.0/6.0;
97 
98  double *ww1, *ww2, *ww3;
99  ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
100  ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
101  ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
102 
103  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
104  dimension "dir" */
105  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
106  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
107  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
108  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
109 
110 #if defined(CPU_STAT)
111  clock_t cpu_start, cpu_end;
112  cpu_start = clock();
113 #endif
114 
115  int i;
116 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
117  for (i=0; i<N_outer; i++) {
118  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
119  _ArrayCopy1D_(index_outer,indexC,ndims);
120  _ArrayCopy1D_(index_outer,indexI,ndims);
121  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
122  int qm1,qm2,qm3,qp1,qp2,p;
123  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
124  if (upw > 0) {
125  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
126  qm3 = qm1 - 2*stride[dir];
127  qm2 = qm1 - stride[dir];
128  qp1 = qm1 + stride[dir];
129  qp2 = qm1 + 2*stride[dir];
130  } else {
131  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
132  qm3 = qm1 + 2*stride[dir];
133  qm2 = qm1 + stride[dir];
134  qp1 = qm1 - stride[dir];
135  qp2 = qm1 - 2*stride[dir];
136  }
137 
138  /* Defining stencil points */
139  double *fm3, *fm2, *fm1, *fp1, *fp2;
140  fm3 = (fC+qm3*nvars);
141  fm2 = (fC+qm2*nvars);
142  fm1 = (fC+qm1*nvars);
143  fp1 = (fC+qp1*nvars);
144  fp2 = (fC+qp2*nvars);
145 
146  /* Candidate stencils and their optimal weights*/
147  double f1[nvars], f2[nvars], f3[nvars];
148  _ArrayAXBYCZ_(f1,(2*one_sixth),fm3,(-7*one_sixth) ,fm2,(11*one_sixth) ,fm1,nvars);
149  _ArrayAXBYCZ_(f2,(-one_sixth) ,fm2,(5*one_sixth) ,fm1,(2*one_sixth) ,fp1,nvars);
150  _ArrayAXBYCZ_(f3,(2*one_sixth),fm1,(5*one_sixth) ,fp1,(-one_sixth) ,fp2,nvars);
151 
152  /* calculate WENO weights */
153  double *w1,*w2,*w3;
154  w1 = (ww1+p*nvars);
155  w2 = (ww2+p*nvars);
156  w3 = (ww3+p*nvars);
157 
158  _ArrayMultiply3Add1D_((fI+p*nvars),w1,f1,w2,f2,w3,f3,nvars);
159  }
160  }
161 
162 #if defined(CPU_STAT)
163  cpu_end = clock();
164  printf("Interp1PrimFifthOrderWENO CPU time = %8.6lf dir = %d\n", (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC, dir);
165 #endif
166 
167  return(0);
168 }
void * interp
Definition: hypar.h:362
Definitions for the functions computing the interpolated value of the primitive at the cell interface...
#define _ArrayIndexnD_(N, index, imax, i, ghost)
int * dim_local
Definition: hypar.h:37
MPI related function definitions.
Structure of variables/parameters needed by the WENO-type scheme.
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayMultiply3Add1D_(x, a, b, c, d, e, f, size)
int Interp1PrimFifthOrderWENO(double *, double *, double *, double *, int, int, void *, void *, int)
5th order WENO reconstruction (component-wise) on a uniform grid
Contains function definitions for common mathematical functions.
int * stride_with_ghosts
Definition: hypar.h:414
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
#define _ArrayAXBYCZ_(w, a, x, b, y, c, z, size)
Contains structure definition for hypar.
Some basic definitions and macros.
int ndims
Definition: hypar.h:26
Contains macros and function definitions for common array operations.
#define _ArrayProduct1D_(x, size, p)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23