HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
WENOFifthOrderCalculateWeights.c
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <string.h>
8 #include <time.h>
9 
10 #include <basic.h>
11 #include <arrayfunctions.h>
12 #include <matmult_native.h>
13 #include <mathfunctions.h>
14 #include <interpolation.h>
15 #include <mpivars.h>
16 #include <hypar.h>
17 
18 #include <arrayfunctions_gpu.h>
19 
20 #undef _MINIMUM_GHOSTS_
21 
25 #define _MINIMUM_GHOSTS_ 3
26 
27 static int WENOFifthOrderCalculateWeightsJS(double*,double*,double*,int,void*,void*);
28 static int WENOFifthOrderCalculateWeightsM (double*,double*,double*,int,void*,void*);
29 static int WENOFifthOrderCalculateWeightsZ (double*,double*,double*,int,void*,void*);
30 static int WENOFifthOrderCalculateWeightsYC(double*,double*,double*,int,void*,void*);
31 
32 static int WENOFifthOrderCalculateWeightsCharJS(double*,double*,double*,int,void*,void*);
33 static int WENOFifthOrderCalculateWeightsCharM (double*,double*,double*,int,void*,void*);
34 static int WENOFifthOrderCalculateWeightsCharZ (double*,double*,double*,int,void*,void*);
35 static int WENOFifthOrderCalculateWeightsCharYC(double*,double*,double*,int,void*,void*);
36 
40 int WENOFifthOrderCalculateWeights(
41  double *fC,
42  double *uC,
43  double *x,
44  int dir,
45  void *s,
46  void *m
47  )
48 {
49  HyPar *solver = (HyPar*) s;
50  WENOParameters *weno = (WENOParameters*) solver->interp;
51  MPIVariables *mpi = (MPIVariables*) m;
52 
53  int ret;
54 
55  if (weno->yc) ret = WENOFifthOrderCalculateWeightsYC (fC,uC,x,dir,solver,mpi);
56  else if (weno->borges) ret = WENOFifthOrderCalculateWeightsZ (fC,uC,x,dir,solver,mpi);
57  else if (weno->mapped) ret = WENOFifthOrderCalculateWeightsM (fC,uC,x,dir,solver,mpi);
58  else ret = WENOFifthOrderCalculateWeightsJS (fC,uC,x,dir,solver,mpi);
59 
60  return ret;
61 }
62 
66 int WENOFifthOrderCalculateWeightsChar(
67  double *fC,
68  double *uC,
69  double *x,
70  int dir,
71  void *s,
72  void *m
73  )
74 {
75  HyPar *solver = (HyPar*) s;
76  WENOParameters *weno = (WENOParameters*) solver->interp;
77  MPIVariables *mpi = (MPIVariables*) m;
78 
79  if (weno->yc) return(WENOFifthOrderCalculateWeightsCharYC (fC,uC,x,dir,solver,mpi));
80  else if (weno->borges) return(WENOFifthOrderCalculateWeightsCharZ (fC,uC,x,dir,solver,mpi));
81  else if (weno->mapped) return(WENOFifthOrderCalculateWeightsCharM (fC,uC,x,dir,solver,mpi));
82  else return(WENOFifthOrderCalculateWeightsCharJS (fC,uC,x,dir,solver,mpi));
83 }
84 
111 int WENOFifthOrderCalculateWeightsJS(
112  double *fC,
113  double *uC,
114  double *x,
115  int dir,
116  void *s,
117  void *m
118  )
119 {
120  HyPar *solver = (HyPar*) s;
121  WENOParameters *weno = (WENOParameters*) solver->interp;
122  MPIVariables *mpi = (MPIVariables*) m;
123  int i;
124  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
125  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
126 
127  int ghosts = solver->ghosts;
128  int ndims = solver->ndims;
129  int nvars = solver->nvars;
130  int *dim = solver->dim_local;
131  int *stride= solver->stride_with_ghosts;
132 
133  /* define some constants */
134  static const double thirteen_by_twelve = 13.0/12.0;
135  static const double one_fourth = 1.0/4.0;
136 
137  /* calculate dimension offset */
138  int offset = weno->offset[dir];
139 
140  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
141  dimension "dir" */
142  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
143  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
144  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
145  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
146 
147  ww1LF = weno->w1 + offset;
148  ww2LF = weno->w2 + offset;
149  ww3LF = weno->w3 + offset;
150  ww1RF = weno->w1 + 2*weno->size + offset;
151  ww2RF = weno->w2 + 2*weno->size + offset;
152  ww3RF = weno->w3 + 2*weno->size + offset;
153  ww1LU = weno->w1 + weno->size + offset;
154  ww2LU = weno->w2 + weno->size + offset;
155  ww3LU = weno->w3 + weno->size + offset;
156  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
157  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
158  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
159 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
160  for (i=0; i<N_outer; i++) {
161  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
162  _ArrayCopy1D_(index_outer,indexC,ndims);
163  _ArrayCopy1D_(index_outer,indexI,ndims);
164  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
165  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
166  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
167  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
168  qm3L = qm1L - 2*stride[dir];
169  qm2L = qm1L - stride[dir];
170  qp1L = qm1L + stride[dir];
171  qp2L = qm1L + 2*stride[dir];
172 
173  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
174  qm3R = qm1R + 2*stride[dir];
175  qm2R = qm1R + stride[dir];
176  qp1R = qm1R - stride[dir];
177  qp2R = qm1R - 2*stride[dir];
178 
179  /* Defining stencil points */
180  double *m3LF, *m2LF, *m1LF, *p1LF, *p2LF;
181  m3LF = (fC+qm3L*nvars);
182  m2LF = (fC+qm2L*nvars);
183  m1LF = (fC+qm1L*nvars);
184  p1LF = (fC+qp1L*nvars);
185  p2LF = (fC+qp2L*nvars);
186  double *m3RF, *m2RF, *m1RF, *p1RF, *p2RF;
187  m3RF = (fC+qm3R*nvars);
188  m2RF = (fC+qm2R*nvars);
189  m1RF = (fC+qm1R*nvars);
190  p1RF = (fC+qp1R*nvars);
191  p2RF = (fC+qp2R*nvars);
192  double *m3LU, *m2LU, *m1LU, *p1LU, *p2LU;
193  m3LU = (uC+qm3L*nvars);
194  m2LU = (uC+qm2L*nvars);
195  m1LU = (uC+qm1L*nvars);
196  p1LU = (uC+qp1L*nvars);
197  p2LU = (uC+qp2L*nvars);
198  double *m3RU, *m2RU, *m1RU, *p1RU, *p2RU;
199  m3RU = (uC+qm3R*nvars);
200  m2RU = (uC+qm2R*nvars);
201  m1RU = (uC+qm1R*nvars);
202  p1RU = (uC+qp1R*nvars);
203  p2RU = (uC+qp2R*nvars);
204 
205  /* optimal weights*/
206  double c1, c2, c3;
207  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
208  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
209  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
210  /* Use WENO5 at the physical boundaries */
211  c1 = _WENO_OPTIMAL_WEIGHT_1_;
212  c2 = _WENO_OPTIMAL_WEIGHT_2_;
213  c3 = _WENO_OPTIMAL_WEIGHT_3_;
214  } else {
215  /* CRWENO5 at the interior points */
216  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
217  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
218  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
219  }
220  } else {
221  /* WENO5 and HCWENO5 */
222  c1 = _WENO_OPTIMAL_WEIGHT_1_;
223  c2 = _WENO_OPTIMAL_WEIGHT_2_;
224  c3 = _WENO_OPTIMAL_WEIGHT_3_;
225  }
226 
227  /* calculate WENO weights */
228  _WENOWeights_v_JS_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
229  _WENOWeights_v_JS_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
230  _WENOWeights_v_JS_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
231  _WENOWeights_v_JS_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
232  }
233  }
234 
235  return(0);
236 }
237 
265 int WENOFifthOrderCalculateWeightsM(
266  double *fC,
267  double *uC,
268  double *x,
269  int dir,
270  void *s,
271  void *m
272  )
273 {
274  HyPar *solver = (HyPar*) s;
275  WENOParameters *weno = (WENOParameters*) solver->interp;
276  MPIVariables *mpi = (MPIVariables*) m;
277  int i;
278  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
279  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
280 
281  int ghosts = solver->ghosts;
282  int ndims = solver->ndims;
283  int nvars = solver->nvars;
284  int *dim = solver->dim_local;
285  int *stride= solver->stride_with_ghosts;
286 
287  /* define some constants */
288  static const double thirteen_by_twelve = 13.0/12.0;
289  static const double one_fourth = 1.0/4.0;
290 
291  /* calculate dimension offset */
292  int offset = weno->offset[dir];
293 
294  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
295  dimension "dir" */
296  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
297  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
298  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
299  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
300 
301  ww1LF = weno->w1 + offset;
302  ww2LF = weno->w2 + offset;
303  ww3LF = weno->w3 + offset;
304  ww1RF = weno->w1 + 2*weno->size + offset;
305  ww2RF = weno->w2 + 2*weno->size + offset;
306  ww3RF = weno->w3 + 2*weno->size + offset;
307  ww1LU = weno->w1 + weno->size + offset;
308  ww2LU = weno->w2 + weno->size + offset;
309  ww3LU = weno->w3 + weno->size + offset;
310  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
311  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
312  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
313 
314 #if defined(CPU_STAT)
315  clock_t startTime, endTime;
316  double copyTime = 0.0;
317  startTime = clock();
318 #endif
319 
320 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
321  for (i=0; i<N_outer; i++) {
322  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
323  _ArrayCopy1D_(index_outer,indexC,ndims);
324  _ArrayCopy1D_(index_outer,indexI,ndims);
325  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
326  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
327  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
328  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
329  qm3L = qm1L - 2*stride[dir];
330  qm2L = qm1L - stride[dir];
331  qp1L = qm1L + stride[dir];
332  qp2L = qm1L + 2*stride[dir];
333 
334  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
335  qm3R = qm1R + 2*stride[dir];
336  qm2R = qm1R + stride[dir];
337  qp1R = qm1R - stride[dir];
338  qp2R = qm1R - 2*stride[dir];
339 
340  /* Defining stencil points */
341  double *m3LF, *m2LF, *m1LF, *p1LF, *p2LF;
342  m3LF = (fC+qm3L*nvars);
343  m2LF = (fC+qm2L*nvars);
344  m1LF = (fC+qm1L*nvars);
345  p1LF = (fC+qp1L*nvars);
346  p2LF = (fC+qp2L*nvars);
347  double *m3RF, *m2RF, *m1RF, *p1RF, *p2RF;
348  m3RF = (fC+qm3R*nvars);
349  m2RF = (fC+qm2R*nvars);
350  m1RF = (fC+qm1R*nvars);
351  p1RF = (fC+qp1R*nvars);
352  p2RF = (fC+qp2R*nvars);
353  double *m3LU, *m2LU, *m1LU, *p1LU, *p2LU;
354  m3LU = (uC+qm3L*nvars);
355  m2LU = (uC+qm2L*nvars);
356  m1LU = (uC+qm1L*nvars);
357  p1LU = (uC+qp1L*nvars);
358  p2LU = (uC+qp2L*nvars);
359  double *m3RU, *m2RU, *m1RU, *p1RU, *p2RU;
360  m3RU = (uC+qm3R*nvars);
361  m2RU = (uC+qm2R*nvars);
362  m1RU = (uC+qm1R*nvars);
363  p1RU = (uC+qp1R*nvars);
364  p2RU = (uC+qp2R*nvars);
365 
366  /* optimal weights*/
367  double c1, c2, c3;
368  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
369  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
370  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
371  /* Use WENO5 at the physical boundaries */
372  c1 = _WENO_OPTIMAL_WEIGHT_1_;
373  c2 = _WENO_OPTIMAL_WEIGHT_2_;
374  c3 = _WENO_OPTIMAL_WEIGHT_3_;
375  } else {
376  /* CRWENO5 at the interior points */
377  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
378  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
379  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
380  }
381  } else {
382  /* WENO5 and HCWENO5 */
383  c1 = _WENO_OPTIMAL_WEIGHT_1_;
384  c2 = _WENO_OPTIMAL_WEIGHT_2_;
385  c3 = _WENO_OPTIMAL_WEIGHT_3_;
386  }
387 
388  /* calculate WENO weights */
389  _WENOWeights_v_M_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
390  _WENOWeights_v_M_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
391  _WENOWeights_v_M_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
392  _WENOWeights_v_M_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
393  }
394  }
395 
396 #if defined(CPU_STAT)
397  endTime = clock();
398  printf("WENOFifthOrderCalculateWeightsM CPU time = %8.6f dir = %d\n",
399  (double)(endTime - startTime) / CLOCKS_PER_SEC, dir);
400 #endif
401 
402  return(0);
403 }
404 
434 int WENOFifthOrderCalculateWeightsZ(
435  double *fC,
436  double *uC,
437  double *x,
438  int dir,
439  void *s,
440  void *m
441  )
442 {
443  HyPar *solver = (HyPar*) s;
444  WENOParameters *weno = (WENOParameters*) solver->interp;
445  MPIVariables *mpi = (MPIVariables*) m;
446  int i;
447  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
448  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
449 
450  int ghosts = solver->ghosts;
451  int ndims = solver->ndims;
452  int nvars = solver->nvars;
453  int *dim = solver->dim_local;
454  int *stride= solver->stride_with_ghosts;
455 
456  /* define some constants */
457  static const double thirteen_by_twelve = 13.0/12.0;
458  static const double one_fourth = 1.0/4.0;
459 
460  /* calculate dimension offset */
461  int offset = weno->offset[dir];
462 
463  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
464  dimension "dir" */
465  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
466  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
467  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
468  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
469 
470  ww1LF = weno->w1 + offset;
471  ww2LF = weno->w2 + offset;
472  ww3LF = weno->w3 + offset;
473  ww1RF = weno->w1 + 2*weno->size + offset;
474  ww2RF = weno->w2 + 2*weno->size + offset;
475  ww3RF = weno->w3 + 2*weno->size + offset;
476  ww1LU = weno->w1 + weno->size + offset;
477  ww2LU = weno->w2 + weno->size + offset;
478  ww3LU = weno->w3 + weno->size + offset;
479  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
480  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
481  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
482 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
483  for (i=0; i<N_outer; i++) {
484  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
485  _ArrayCopy1D_(index_outer,indexC,ndims);
486  _ArrayCopy1D_(index_outer,indexI,ndims);
487  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
488  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
489  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
490  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
491  qm3L = qm1L - 2*stride[dir];
492  qm2L = qm1L - stride[dir];
493  qp1L = qm1L + stride[dir];
494  qp2L = qm1L + 2*stride[dir];
495 
496  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
497  qm3R = qm1R + 2*stride[dir];
498  qm2R = qm1R + stride[dir];
499  qp1R = qm1R - stride[dir];
500  qp2R = qm1R - 2*stride[dir];
501 
502  /* Defining stencil points */
503  double *m3LF, *m2LF, *m1LF, *p1LF, *p2LF;
504  m3LF = (fC+qm3L*nvars);
505  m2LF = (fC+qm2L*nvars);
506  m1LF = (fC+qm1L*nvars);
507  p1LF = (fC+qp1L*nvars);
508  p2LF = (fC+qp2L*nvars);
509  double *m3RF, *m2RF, *m1RF, *p1RF, *p2RF;
510  m3RF = (fC+qm3R*nvars);
511  m2RF = (fC+qm2R*nvars);
512  m1RF = (fC+qm1R*nvars);
513  p1RF = (fC+qp1R*nvars);
514  p2RF = (fC+qp2R*nvars);
515  double *m3LU, *m2LU, *m1LU, *p1LU, *p2LU;
516  m3LU = (uC+qm3L*nvars);
517  m2LU = (uC+qm2L*nvars);
518  m1LU = (uC+qm1L*nvars);
519  p1LU = (uC+qp1L*nvars);
520  p2LU = (uC+qp2L*nvars);
521  double *m3RU, *m2RU, *m1RU, *p1RU, *p2RU;
522  m3RU = (uC+qm3R*nvars);
523  m2RU = (uC+qm2R*nvars);
524  m1RU = (uC+qm1R*nvars);
525  p1RU = (uC+qp1R*nvars);
526  p2RU = (uC+qp2R*nvars);
527 
528  /* optimal weights*/
529  double c1, c2, c3;
530  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
531  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
532  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
533  /* Use WENO5 at the physical boundaries */
534  c1 = _WENO_OPTIMAL_WEIGHT_1_;
535  c2 = _WENO_OPTIMAL_WEIGHT_2_;
536  c3 = _WENO_OPTIMAL_WEIGHT_3_;
537  } else {
538  /* CRWENO5 at the interior points */
539  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
540  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
541  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
542  }
543  } else {
544  /* WENO5 and HCWENO5 */
545  c1 = _WENO_OPTIMAL_WEIGHT_1_;
546  c2 = _WENO_OPTIMAL_WEIGHT_2_;
547  c3 = _WENO_OPTIMAL_WEIGHT_3_;
548  }
549 
550  /* calculate WENO weights */
551  _WENOWeights_v_Z_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
552  _WENOWeights_v_Z_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
553  _WENOWeights_v_Z_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
554  _WENOWeights_v_Z_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
555  }
556  }
557 
558  return(0);
559 }
560 
590 int WENOFifthOrderCalculateWeightsYC(
591  double *fC,
592  double *uC,
593  double *x,
594  int dir,
595  void *s,
596  void *m
597  )
598 {
599  HyPar *solver = (HyPar*) s;
600  WENOParameters *weno = (WENOParameters*) solver->interp;
601  MPIVariables *mpi = (MPIVariables*) m;
602  int i;
603  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
604  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
605 
606  int ghosts = solver->ghosts;
607  int ndims = solver->ndims;
608  int nvars = solver->nvars;
609  int *dim = solver->dim_local;
610  int *stride= solver->stride_with_ghosts;
611 
612  /* define some constants */
613  static const double thirteen_by_twelve = 13.0/12.0;
614  static const double one_fourth = 1.0/4.0;
615 
616  /* calculate dimension offset */
617  int offset = weno->offset[dir];
618 
619  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
620  dimension "dir" */
621  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
622  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
623  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
624  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
625 
626  ww1LF = weno->w1 + offset;
627  ww2LF = weno->w2 + offset;
628  ww3LF = weno->w3 + offset;
629  ww1RF = weno->w1 + 2*weno->size + offset;
630  ww2RF = weno->w2 + 2*weno->size + offset;
631  ww3RF = weno->w3 + 2*weno->size + offset;
632  ww1LU = weno->w1 + weno->size + offset;
633  ww2LU = weno->w2 + weno->size + offset;
634  ww3LU = weno->w3 + weno->size + offset;
635  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
636  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
637  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
638 
639 #if defined(CPU_STAT)
640  clock_t cpu_start, cpu_end;
641  cpu_start = clock();
642 #endif
643 
644 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
645  for (i=0; i<N_outer; i++) {
646  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
647  _ArrayCopy1D_(index_outer,indexC,ndims);
648  _ArrayCopy1D_(index_outer,indexI,ndims);
649  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
650  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
651  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
652  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
653  qm3L = qm1L - 2*stride[dir];
654  qm2L = qm1L - stride[dir];
655  qp1L = qm1L + stride[dir];
656  qp2L = qm1L + 2*stride[dir];
657 
658  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
659  qm3R = qm1R + 2*stride[dir];
660  qm2R = qm1R + stride[dir];
661  qp1R = qm1R - stride[dir];
662  qp2R = qm1R - 2*stride[dir];
663 
664  /* Defining stencil points */
665  double *m3LF, *m2LF, *m1LF, *p1LF, *p2LF;
666  m3LF = (fC+qm3L*nvars);
667  m2LF = (fC+qm2L*nvars);
668  m1LF = (fC+qm1L*nvars);
669  p1LF = (fC+qp1L*nvars);
670  p2LF = (fC+qp2L*nvars);
671  double *m3RF, *m2RF, *m1RF, *p1RF, *p2RF;
672  m3RF = (fC+qm3R*nvars);
673  m2RF = (fC+qm2R*nvars);
674  m1RF = (fC+qm1R*nvars);
675  p1RF = (fC+qp1R*nvars);
676  p2RF = (fC+qp2R*nvars);
677  double *m3LU, *m2LU, *m1LU, *p1LU, *p2LU;
678  m3LU = (uC+qm3L*nvars);
679  m2LU = (uC+qm2L*nvars);
680  m1LU = (uC+qm1L*nvars);
681  p1LU = (uC+qp1L*nvars);
682  p2LU = (uC+qp2L*nvars);
683  double *m3RU, *m2RU, *m1RU, *p1RU, *p2RU;
684  m3RU = (uC+qm3R*nvars);
685  m2RU = (uC+qm2R*nvars);
686  m1RU = (uC+qm1R*nvars);
687  p1RU = (uC+qp1R*nvars);
688  p2RU = (uC+qp2R*nvars);
689 
690  /* optimal weights*/
691  double c1, c2, c3;
692  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
693  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
694  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
695  /* Use WENO5 at the physical boundaries */
696  c1 = _WENO_OPTIMAL_WEIGHT_1_;
697  c2 = _WENO_OPTIMAL_WEIGHT_2_;
698  c3 = _WENO_OPTIMAL_WEIGHT_3_;
699  } else {
700  /* CRWENO5 at the interior points */
701  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
702  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
703  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
704  }
705  } else {
706  /* WENO5 and HCWENO5 */
707  c1 = _WENO_OPTIMAL_WEIGHT_1_;
708  c2 = _WENO_OPTIMAL_WEIGHT_2_;
709  c3 = _WENO_OPTIMAL_WEIGHT_3_;
710  }
711 
712  /* calculate WENO weights */
713  _WENOWeights_v_YC_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
714  _WENOWeights_v_YC_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
715  _WENOWeights_v_YC_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
716  _WENOWeights_v_YC_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
717  }
718  }
719 
720 #if defined(CPU_STAT)
721  cpu_end = clock();
722  printf("WENOFifthOrderCalculateWeightsYC:\n");
723  printf(" CPU time = %8.6lf dir = %d\n", (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC, dir);
724 #endif
725 
726  return(0);
727 }
728 
760 int WENOFifthOrderCalculateWeightsCharJS(
761  double *fC,
762  double *uC,
763  double *x,
764  int dir,
765  void *s,
766  void *m
767  )
768 {
769  HyPar *solver = (HyPar*) s;
770  WENOParameters *weno = (WENOParameters*) solver->interp;
771  MPIVariables *mpi = (MPIVariables*) m;
772  int i;
773  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
774  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
775 
776  int ghosts = solver->ghosts;
777  int ndims = solver->ndims;
778  int nvars = solver->nvars;
779  int *dim = solver->dim_local;
780  int *stride= solver->stride_with_ghosts;
781 
782  /* define some constants */
783  static const double thirteen_by_twelve = 13.0/12.0;
784  static const double one_fourth = 1.0/4.0;
785 
786  /* calculate dimension offset */
787  int offset = weno->offset[dir];
788 
789  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
790  dimension "dir" */
791  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
792  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
793  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
794  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
795 
796  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
797  double L[nvars*nvars], uavg[nvars];
798 
799  ww1LF = weno->w1 + offset;
800  ww2LF = weno->w2 + offset;
801  ww3LF = weno->w3 + offset;
802  ww1RF = weno->w1 + 2*weno->size + offset;
803  ww2RF = weno->w2 + 2*weno->size + offset;
804  ww3RF = weno->w3 + 2*weno->size + offset;
805  ww1LU = weno->w1 + weno->size + offset;
806  ww2LU = weno->w2 + weno->size + offset;
807  ww3LU = weno->w3 + weno->size + offset;
808  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
809  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
810  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
811 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
812  for (i=0; i<N_outer; i++) {
813  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
814  _ArrayCopy1D_(index_outer,indexC,ndims);
815  _ArrayCopy1D_(index_outer,indexI,ndims);
816  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
817  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
818  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
819  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
820  qm3L = qm1L - 2*stride[dir];
821  qm2L = qm1L - stride[dir];
822  qp1L = qm1L + stride[dir];
823  qp2L = qm1L + 2*stride[dir];
824 
825  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
826  qm3R = qm1R + 2*stride[dir];
827  qm2R = qm1R + stride[dir];
828  qp1R = qm1R - stride[dir];
829  qp2R = qm1R - 2*stride[dir];
830 
831  /* find averaged state and left eigenvectors at this interface */
832  IERR solver->AveragingFunction(uavg,(uC+nvars*qm1L),(uC+nvars*qp1L),solver->physics); CHECKERR(ierr);
833  IERR solver->GetLeftEigenvectors(uavg,L,solver->physics,dir); CHECKERR(ierr);
834 
835  /* Defining stencil points */
836  double m3LF[nvars], m2LF[nvars], m1LF[nvars], p1LF[nvars], p2LF[nvars];
837  double m3RF[nvars], m2RF[nvars], m1RF[nvars], p1RF[nvars], p2RF[nvars];
838  double m3LU[nvars], m2LU[nvars], m1LU[nvars], p1LU[nvars], p2LU[nvars];
839  double m3RU[nvars], m2RU[nvars], m1RU[nvars], p1RU[nvars], p2RU[nvars];
840 
841  MatVecMult(nvars,m3LF,L,(fC+nvars*qm3L));
842  MatVecMult(nvars,m2LF,L,(fC+nvars*qm2L));
843  MatVecMult(nvars,m1LF,L,(fC+nvars*qm1L));
844  MatVecMult(nvars,p1LF,L,(fC+nvars*qp1L));
845  MatVecMult(nvars,p2LF,L,(fC+nvars*qp2L));
846 
847  MatVecMult(nvars,m3RF,L,(fC+nvars*qm3R));
848  MatVecMult(nvars,m2RF,L,(fC+nvars*qm2R));
849  MatVecMult(nvars,m1RF,L,(fC+nvars*qm1R));
850  MatVecMult(nvars,p1RF,L,(fC+nvars*qp1R));
851  MatVecMult(nvars,p2RF,L,(fC+nvars*qp2R));
852 
853  MatVecMult(nvars,m3LU,L,(uC+nvars*qm3L));
854  MatVecMult(nvars,m2LU,L,(uC+nvars*qm2L));
855  MatVecMult(nvars,m1LU,L,(uC+nvars*qm1L));
856  MatVecMult(nvars,p1LU,L,(uC+nvars*qp1L));
857  MatVecMult(nvars,p2LU,L,(uC+nvars*qp2L));
858 
859  MatVecMult(nvars,m3RU,L,(uC+nvars*qm3R));
860  MatVecMult(nvars,m2RU,L,(uC+nvars*qm2R));
861  MatVecMult(nvars,m1RU,L,(uC+nvars*qm1R));
862  MatVecMult(nvars,p1RU,L,(uC+nvars*qp1R));
863  MatVecMult(nvars,p2RU,L,(uC+nvars*qp2R));
864 
865  /* optimal weights*/
866  double c1, c2, c3;
867  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
868  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
869  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
870  /* Use WENO5 at the physical boundaries */
871  c1 = _WENO_OPTIMAL_WEIGHT_1_;
872  c2 = _WENO_OPTIMAL_WEIGHT_2_;
873  c3 = _WENO_OPTIMAL_WEIGHT_3_;
874  } else {
875  /* CRWENO5 at the interior points */
876  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
877  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
878  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
879  }
880  } else {
881  /* WENO5 and HCWENO5 */
882  c1 = _WENO_OPTIMAL_WEIGHT_1_;
883  c2 = _WENO_OPTIMAL_WEIGHT_2_;
884  c3 = _WENO_OPTIMAL_WEIGHT_3_;
885  }
886 
887  /* calculate WENO weights */
888  _WENOWeights_v_JS_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
889  _WENOWeights_v_JS_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
890  _WENOWeights_v_JS_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
891  _WENOWeights_v_JS_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
892  }
893  }
894 
895  return(0);
896 }
897 
930 int WENOFifthOrderCalculateWeightsCharM(
931  double *fC,
932  double *uC,
933  double *x,
934  int dir,
935  void *s,
936  void *m
937  )
938 {
939  HyPar *solver = (HyPar*) s;
940  WENOParameters *weno = (WENOParameters*) solver->interp;
941  MPIVariables *mpi = (MPIVariables*) m;
942  int i;
943  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
944  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
945 
946  int ghosts = solver->ghosts;
947  int ndims = solver->ndims;
948  int nvars = solver->nvars;
949  int *dim = solver->dim_local;
950  int *stride= solver->stride_with_ghosts;
951 
952  /* define some constants */
953  static const double thirteen_by_twelve = 13.0/12.0;
954  static const double one_fourth = 1.0/4.0;
955 
956  /* calculate dimension offset */
957  int offset = weno->offset[dir];
958 
959  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
960  dimension "dir" */
961  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
962  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
963  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
964  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
965 
966  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
967  double L[nvars*nvars], uavg[nvars];
968 
969  ww1LF = weno->w1 + offset;
970  ww2LF = weno->w2 + offset;
971  ww3LF = weno->w3 + offset;
972  ww1RF = weno->w1 + 2*weno->size + offset;
973  ww2RF = weno->w2 + 2*weno->size + offset;
974  ww3RF = weno->w3 + 2*weno->size + offset;
975  ww1LU = weno->w1 + weno->size + offset;
976  ww2LU = weno->w2 + weno->size + offset;
977  ww3LU = weno->w3 + weno->size + offset;
978  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
979  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
980  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
981 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
982  for (i=0; i<N_outer; i++) {
983  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
984  _ArrayCopy1D_(index_outer,indexC,ndims);
985  _ArrayCopy1D_(index_outer,indexI,ndims);
986  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
987  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
988  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
989  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
990  qm3L = qm1L - 2*stride[dir];
991  qm2L = qm1L - stride[dir];
992  qp1L = qm1L + stride[dir];
993  qp2L = qm1L + 2*stride[dir];
994 
995  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
996  qm3R = qm1R + 2*stride[dir];
997  qm2R = qm1R + stride[dir];
998  qp1R = qm1R - stride[dir];
999  qp2R = qm1R - 2*stride[dir];
1000 
1001  /* find averaged state and left eigenvectors at this interface */
1002  IERR solver->AveragingFunction(uavg,(uC+nvars*qm1L),(uC+nvars*qp1L),solver->physics); CHECKERR(ierr);
1003  IERR solver->GetLeftEigenvectors(uavg,L,solver->physics,dir); CHECKERR(ierr);
1004 
1005  /* Defining stencil points */
1006  double m3LF[nvars], m2LF[nvars], m1LF[nvars], p1LF[nvars], p2LF[nvars];
1007  double m3RF[nvars], m2RF[nvars], m1RF[nvars], p1RF[nvars], p2RF[nvars];
1008  double m3LU[nvars], m2LU[nvars], m1LU[nvars], p1LU[nvars], p2LU[nvars];
1009  double m3RU[nvars], m2RU[nvars], m1RU[nvars], p1RU[nvars], p2RU[nvars];
1010 
1011  MatVecMult(nvars,m3LF,L,(fC+nvars*qm3L));
1012  MatVecMult(nvars,m2LF,L,(fC+nvars*qm2L));
1013  MatVecMult(nvars,m1LF,L,(fC+nvars*qm1L));
1014  MatVecMult(nvars,p1LF,L,(fC+nvars*qp1L));
1015  MatVecMult(nvars,p2LF,L,(fC+nvars*qp2L));
1016 
1017  MatVecMult(nvars,m3RF,L,(fC+nvars*qm3R));
1018  MatVecMult(nvars,m2RF,L,(fC+nvars*qm2R));
1019  MatVecMult(nvars,m1RF,L,(fC+nvars*qm1R));
1020  MatVecMult(nvars,p1RF,L,(fC+nvars*qp1R));
1021  MatVecMult(nvars,p2RF,L,(fC+nvars*qp2R));
1022 
1023  MatVecMult(nvars,m3LU,L,(uC+nvars*qm3L));
1024  MatVecMult(nvars,m2LU,L,(uC+nvars*qm2L));
1025  MatVecMult(nvars,m1LU,L,(uC+nvars*qm1L));
1026  MatVecMult(nvars,p1LU,L,(uC+nvars*qp1L));
1027  MatVecMult(nvars,p2LU,L,(uC+nvars*qp2L));
1028 
1029  MatVecMult(nvars,m3RU,L,(uC+nvars*qm3R));
1030  MatVecMult(nvars,m2RU,L,(uC+nvars*qm2R));
1031  MatVecMult(nvars,m1RU,L,(uC+nvars*qm1R));
1032  MatVecMult(nvars,p1RU,L,(uC+nvars*qp1R));
1033  MatVecMult(nvars,p2RU,L,(uC+nvars*qp2R));
1034 
1035  /* optimal weights*/
1036  double c1, c2, c3;
1037  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
1038  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
1039  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
1040  /* Use WENO5 at the physical boundaries */
1041  c1 = _WENO_OPTIMAL_WEIGHT_1_;
1042  c2 = _WENO_OPTIMAL_WEIGHT_2_;
1043  c3 = _WENO_OPTIMAL_WEIGHT_3_;
1044  } else {
1045  /* CRWENO5 at the interior points */
1046  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
1047  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
1048  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
1049  }
1050  } else {
1051  /* WENO5 and HCWENO5 */
1052  c1 = _WENO_OPTIMAL_WEIGHT_1_;
1053  c2 = _WENO_OPTIMAL_WEIGHT_2_;
1054  c3 = _WENO_OPTIMAL_WEIGHT_3_;
1055  }
1056 
1057  /* calculate WENO weights */
1058  _WENOWeights_v_M_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
1059  _WENOWeights_v_M_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
1060  _WENOWeights_v_M_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
1061  _WENOWeights_v_M_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
1062  }
1063  }
1064 
1065  return(0);
1066 }
1067 
1101 int WENOFifthOrderCalculateWeightsCharZ(
1102  double *fC,
1103  double *uC,
1104  double *x,
1105  int dir,
1106  void *s,
1107  void *m
1108  )
1109 {
1110  HyPar *solver = (HyPar*) s;
1111  WENOParameters *weno = (WENOParameters*) solver->interp;
1112  MPIVariables *mpi = (MPIVariables*) m;
1113  int i;
1114  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
1115  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
1116 
1117  int ghosts = solver->ghosts;
1118  int ndims = solver->ndims;
1119  int nvars = solver->nvars;
1120  int *dim = solver->dim_local;
1121  int *stride= solver->stride_with_ghosts;
1122 
1123  /* define some constants */
1124  static const double thirteen_by_twelve = 13.0/12.0;
1125  static const double one_fourth = 1.0/4.0;
1126 
1127  /* calculate dimension offset */
1128  int offset = weno->offset[dir];
1129 
1130  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
1131  dimension "dir" */
1132  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
1133  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
1134  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
1135  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
1136 
1137  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
1138  double L[nvars*nvars], uavg[nvars];
1139 
1140  ww1LF = weno->w1 + offset;
1141  ww2LF = weno->w2 + offset;
1142  ww3LF = weno->w3 + offset;
1143  ww1RF = weno->w1 + 2*weno->size + offset;
1144  ww2RF = weno->w2 + 2*weno->size + offset;
1145  ww3RF = weno->w3 + 2*weno->size + offset;
1146  ww1LU = weno->w1 + weno->size + offset;
1147  ww2LU = weno->w2 + weno->size + offset;
1148  ww3LU = weno->w3 + weno->size + offset;
1149  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
1150  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
1151  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
1152 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
1153  for (i=0; i<N_outer; i++) {
1154  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
1155  _ArrayCopy1D_(index_outer,indexC,ndims);
1156  _ArrayCopy1D_(index_outer,indexI,ndims);
1157  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
1158  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
1159  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
1160  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
1161  qm3L = qm1L - 2*stride[dir];
1162  qm2L = qm1L - stride[dir];
1163  qp1L = qm1L + stride[dir];
1164  qp2L = qm1L + 2*stride[dir];
1165 
1166  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
1167  qm3R = qm1R + 2*stride[dir];
1168  qm2R = qm1R + stride[dir];
1169  qp1R = qm1R - stride[dir];
1170  qp2R = qm1R - 2*stride[dir];
1171 
1172  /* find averaged state and left eigenvectors at this interface */
1173  IERR solver->AveragingFunction(uavg,(uC+nvars*qm1L),(uC+nvars*qp1L),solver->physics); CHECKERR(ierr);
1174  IERR solver->GetLeftEigenvectors(uavg,L,solver->physics,dir); CHECKERR(ierr);
1175 
1176  /* Defining stencil points */
1177  double m3LF[nvars], m2LF[nvars], m1LF[nvars], p1LF[nvars], p2LF[nvars];
1178  double m3RF[nvars], m2RF[nvars], m1RF[nvars], p1RF[nvars], p2RF[nvars];
1179  double m3LU[nvars], m2LU[nvars], m1LU[nvars], p1LU[nvars], p2LU[nvars];
1180  double m3RU[nvars], m2RU[nvars], m1RU[nvars], p1RU[nvars], p2RU[nvars];
1181 
1182  MatVecMult(nvars,m3LF,L,(fC+nvars*qm3L));
1183  MatVecMult(nvars,m2LF,L,(fC+nvars*qm2L));
1184  MatVecMult(nvars,m1LF,L,(fC+nvars*qm1L));
1185  MatVecMult(nvars,p1LF,L,(fC+nvars*qp1L));
1186  MatVecMult(nvars,p2LF,L,(fC+nvars*qp2L));
1187 
1188  MatVecMult(nvars,m3RF,L,(fC+nvars*qm3R));
1189  MatVecMult(nvars,m2RF,L,(fC+nvars*qm2R));
1190  MatVecMult(nvars,m1RF,L,(fC+nvars*qm1R));
1191  MatVecMult(nvars,p1RF,L,(fC+nvars*qp1R));
1192  MatVecMult(nvars,p2RF,L,(fC+nvars*qp2R));
1193 
1194  MatVecMult(nvars,m3LU,L,(uC+nvars*qm3L));
1195  MatVecMult(nvars,m2LU,L,(uC+nvars*qm2L));
1196  MatVecMult(nvars,m1LU,L,(uC+nvars*qm1L));
1197  MatVecMult(nvars,p1LU,L,(uC+nvars*qp1L));
1198  MatVecMult(nvars,p2LU,L,(uC+nvars*qp2L));
1199 
1200  MatVecMult(nvars,m3RU,L,(uC+nvars*qm3R));
1201  MatVecMult(nvars,m2RU,L,(uC+nvars*qm2R));
1202  MatVecMult(nvars,m1RU,L,(uC+nvars*qm1R));
1203  MatVecMult(nvars,p1RU,L,(uC+nvars*qp1R));
1204  MatVecMult(nvars,p2RU,L,(uC+nvars*qp2R));
1205 
1206  /* optimal weights*/
1207  double c1, c2, c3;
1208  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
1209  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
1210  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
1211  /* Use WENO5 at the physical boundaries */
1212  c1 = _WENO_OPTIMAL_WEIGHT_1_;
1213  c2 = _WENO_OPTIMAL_WEIGHT_2_;
1214  c3 = _WENO_OPTIMAL_WEIGHT_3_;
1215  } else {
1216  /* CRWENO5 at the interior points */
1217  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
1218  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
1219  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
1220  }
1221  } else {
1222  /* WENO5 and HCWENO5 */
1223  c1 = _WENO_OPTIMAL_WEIGHT_1_;
1224  c2 = _WENO_OPTIMAL_WEIGHT_2_;
1225  c3 = _WENO_OPTIMAL_WEIGHT_3_;
1226  }
1227 
1228  /* calculate WENO weights */
1229  _WENOWeights_v_Z_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
1230  _WENOWeights_v_Z_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
1231  _WENOWeights_v_Z_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
1232  _WENOWeights_v_Z_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
1233  }
1234  }
1235 
1236  return(0);
1237 }
1238 
1272 int WENOFifthOrderCalculateWeightsCharYC(
1273  double *fC,
1274  double *uC,
1275  double *x,
1276  int dir,
1277  void *s,
1278  void *m
1279  )
1280 {
1281  HyPar *solver = (HyPar*) s;
1282  WENOParameters *weno = (WENOParameters*) solver->interp;
1283  MPIVariables *mpi = (MPIVariables*) m;
1284  int i;
1285  double *ww1LF, *ww2LF, *ww3LF, *ww1RF, *ww2RF, *ww3RF;
1286  double *ww1LU, *ww2LU, *ww3LU, *ww1RU, *ww2RU, *ww3RU;
1287 
1288  int ghosts = solver->ghosts;
1289  int ndims = solver->ndims;
1290  int nvars = solver->nvars;
1291  int *dim = solver->dim_local;
1292  int *stride= solver->stride_with_ghosts;
1293 
1294  /* define some constants */
1295  static const double thirteen_by_twelve = 13.0/12.0;
1296  static const double one_fourth = 1.0/4.0;
1297 
1298  /* calculate dimension offset */
1299  int offset = weno->offset[dir];
1300 
1301  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
1302  dimension "dir" */
1303  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
1304  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
1305  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
1306  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
1307 
1308  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
1309  double L[nvars*nvars], uavg[nvars];
1310 
1311  ww1LF = weno->w1 + offset;
1312  ww2LF = weno->w2 + offset;
1313  ww3LF = weno->w3 + offset;
1314  ww1RF = weno->w1 + 2*weno->size + offset;
1315  ww2RF = weno->w2 + 2*weno->size + offset;
1316  ww3RF = weno->w3 + 2*weno->size + offset;
1317  ww1LU = weno->w1 + weno->size + offset;
1318  ww2LU = weno->w2 + weno->size + offset;
1319  ww3LU = weno->w3 + weno->size + offset;
1320  ww1RU = weno->w1 + 2*weno->size + weno->size + offset;
1321  ww2RU = weno->w2 + 2*weno->size + weno->size + offset;
1322  ww3RU = weno->w3 + 2*weno->size + weno->size + offset;
1323 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
1324  for (i=0; i<N_outer; i++) {
1325  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
1326  _ArrayCopy1D_(index_outer,indexC,ndims);
1327  _ArrayCopy1D_(index_outer,indexI,ndims);
1328  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
1329  int qm1L,qm2L,qm3L,qp1L,qp2L,p,qm1R,qm2R,qm3R,qp1R,qp2R;
1330  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
1331  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1L);
1332  qm3L = qm1L - 2*stride[dir];
1333  qm2L = qm1L - stride[dir];
1334  qp1L = qm1L + stride[dir];
1335  qp2L = qm1L + 2*stride[dir];
1336 
1337  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1R);
1338  qm3R = qm1R + 2*stride[dir];
1339  qm2R = qm1R + stride[dir];
1340  qp1R = qm1R - stride[dir];
1341  qp2R = qm1R - 2*stride[dir];
1342 
1343  /* find averaged state and left eigenvectors at this interface */
1344  IERR solver->AveragingFunction(uavg,(uC+nvars*qm1L),(uC+nvars*qp1L),solver->physics); CHECKERR(ierr);
1345  IERR solver->GetLeftEigenvectors(uavg,L,solver->physics,dir); CHECKERR(ierr);
1346 
1347  /* Defining stencil points */
1348  double m3LF[nvars], m2LF[nvars], m1LF[nvars], p1LF[nvars], p2LF[nvars];
1349  double m3RF[nvars], m2RF[nvars], m1RF[nvars], p1RF[nvars], p2RF[nvars];
1350  double m3LU[nvars], m2LU[nvars], m1LU[nvars], p1LU[nvars], p2LU[nvars];
1351  double m3RU[nvars], m2RU[nvars], m1RU[nvars], p1RU[nvars], p2RU[nvars];
1352 
1353  MatVecMult(nvars,m3LF,L,(fC+nvars*qm3L));
1354  MatVecMult(nvars,m2LF,L,(fC+nvars*qm2L));
1355  MatVecMult(nvars,m1LF,L,(fC+nvars*qm1L));
1356  MatVecMult(nvars,p1LF,L,(fC+nvars*qp1L));
1357  MatVecMult(nvars,p2LF,L,(fC+nvars*qp2L));
1358 
1359  MatVecMult(nvars,m3RF,L,(fC+nvars*qm3R));
1360  MatVecMult(nvars,m2RF,L,(fC+nvars*qm2R));
1361  MatVecMult(nvars,m1RF,L,(fC+nvars*qm1R));
1362  MatVecMult(nvars,p1RF,L,(fC+nvars*qp1R));
1363  MatVecMult(nvars,p2RF,L,(fC+nvars*qp2R));
1364 
1365  MatVecMult(nvars,m3LU,L,(uC+nvars*qm3L));
1366  MatVecMult(nvars,m2LU,L,(uC+nvars*qm2L));
1367  MatVecMult(nvars,m1LU,L,(uC+nvars*qm1L));
1368  MatVecMult(nvars,p1LU,L,(uC+nvars*qp1L));
1369  MatVecMult(nvars,p2LU,L,(uC+nvars*qp2L));
1370 
1371  MatVecMult(nvars,m3RU,L,(uC+nvars*qm3R));
1372  MatVecMult(nvars,m2RU,L,(uC+nvars*qm2R));
1373  MatVecMult(nvars,m1RU,L,(uC+nvars*qm1R));
1374  MatVecMult(nvars,p1RU,L,(uC+nvars*qp1R));
1375  MatVecMult(nvars,p2RU,L,(uC+nvars*qp2R));
1376 
1377  /* optimal weights*/
1378  double c1, c2, c3;
1379  if (!strcmp(solver->spatial_scheme_hyp,_FIFTH_ORDER_CRWENO_)) {
1380  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
1381  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
1382  /* Use WENO5 at the physical boundaries */
1383  c1 = _WENO_OPTIMAL_WEIGHT_1_;
1384  c2 = _WENO_OPTIMAL_WEIGHT_2_;
1385  c3 = _WENO_OPTIMAL_WEIGHT_3_;
1386  } else {
1387  /* CRWENO5 at the interior points */
1388  c1 = _CRWENO_OPTIMAL_WEIGHT_1_;
1389  c2 = _CRWENO_OPTIMAL_WEIGHT_2_;
1390  c3 = _CRWENO_OPTIMAL_WEIGHT_3_;
1391  }
1392  } else {
1393  /* WENO5 and HCWENO5 */
1394  c1 = _WENO_OPTIMAL_WEIGHT_1_;
1395  c2 = _WENO_OPTIMAL_WEIGHT_2_;
1396  c3 = _WENO_OPTIMAL_WEIGHT_3_;
1397  }
1398 
1399  /* calculate WENO weights */
1400  _WENOWeights_v_YC_((ww1LF+p*nvars),(ww2LF+p*nvars),(ww3LF+p*nvars),c1,c2,c3,m3LF,m2LF,m1LF,p1LF,p2LF,weno->eps,nvars);
1401  _WENOWeights_v_YC_((ww1RF+p*nvars),(ww2RF+p*nvars),(ww3RF+p*nvars),c1,c2,c3,m3RF,m2RF,m1RF,p1RF,p2RF,weno->eps,nvars);
1402  _WENOWeights_v_YC_((ww1LU+p*nvars),(ww2LU+p*nvars),(ww3LU+p*nvars),c1,c2,c3,m3LU,m2LU,m1LU,p1LU,p2LU,weno->eps,nvars);
1403  _WENOWeights_v_YC_((ww1RU+p*nvars),(ww2RU+p*nvars),(ww3RU+p*nvars),c1,c2,c3,m3RU,m2RU,m1RU,p1RU,p2RU,weno->eps,nvars);
1404  }
1405  }
1406 
1407  return(0);
1408 }
HyPar::spatial_scheme_hyp
char spatial_scheme_hyp[_MAX_STRING_SIZE_]
Definition: hypar.h:84
HyPar::interp
void * interp
Definition: hypar.h:362
_WENOWeights_v_M_
#define _WENOWeights_v_M_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)
Definition: interpolation.h:314
WENOFifthOrderCalculateWeightsChar
int WENOFifthOrderCalculateWeightsChar(double *fC, double *uC, double *x, int dir, void *s, void *m)
Definition: WENOFifthOrderCalculateWeights.c:66
interpolation.h
Definitions for the functions computing the interpolated value of the primitive at the cell interface...
WENOFifthOrderCalculateWeights
int WENOFifthOrderCalculateWeights(double *fC, double *uC, double *x, int dir, void *s, void *m)
Definition: WENOFifthOrderCalculateWeights.c:40
_ArrayIndexnD_
#define _ArrayIndexnD_(N, index, imax, i, ghost)
Definition: arrayfunctions.h:21
WENOFifthOrderCalculateWeightsCharJS
static int WENOFifthOrderCalculateWeightsCharJS(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:760
HyPar::physics
void * physics
Definition: hypar.h:266
HyPar::dim_local
int * dim_local
Definition: hypar.h:37
mpivars.h
MPI related function definitions.
_WENO_OPTIMAL_WEIGHT_2_
#define _WENO_OPTIMAL_WEIGHT_2_
Definition: interpolation.h:230
WENOParameters
Structure of variables/parameters needed by the WENO-type scheme.
Definition: interpolation.h:194
HyPar::ghosts
int ghosts
Definition: hypar.h:52
_ArrayIndex1D_
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Definition: arrayfunctions.h:40
_CRWENO_OPTIMAL_WEIGHT_1_
#define _CRWENO_OPTIMAL_WEIGHT_1_
Definition: interpolation.h:234
_WENO_OPTIMAL_WEIGHT_3_
#define _WENO_OPTIMAL_WEIGHT_3_
Definition: interpolation.h:232
WENOFifthOrderCalculateWeightsM
static int WENOFifthOrderCalculateWeightsM(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:265
HyPar::AveragingFunction
int(* AveragingFunction)(double *, double *, double *, void *)
Definition: hypar.h:354
MatVecMult
#define MatVecMult(N, y, A, x)
Definition: matmult_native.h:116
mathfunctions.h
Contains function definitions for common mathematical functions.
WENOFifthOrderCalculateWeightsJS
static int WENOFifthOrderCalculateWeightsJS(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:111
_WENOWeights_v_JS_
#define _WENOWeights_v_JS_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)
Definition: interpolation.h:265
HyPar::GetLeftEigenvectors
int(* GetLeftEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:357
_WENO_OPTIMAL_WEIGHT_1_
#define _WENO_OPTIMAL_WEIGHT_1_
Definition: interpolation.h:228
arrayfunctions_gpu.h
Contains function definitions for common array operations on GPU.
HyPar::stride_with_ghosts
int * stride_with_ghosts
Definition: hypar.h:414
WENOFifthOrderCalculateWeightsCharZ
static int WENOFifthOrderCalculateWeightsCharZ(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:1101
WENOFifthOrderCalculateWeightsYC
static int WENOFifthOrderCalculateWeightsYC(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:590
_WENOWeights_v_Z_
#define _WENOWeights_v_Z_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)
Definition: interpolation.h:426
_ArrayCopy1D_
#define _ArrayCopy1D_(x, y, size)
Definition: arrayfunctions.h:319
HyPar::nvars
int nvars
Definition: hypar.h:29
CHECKERR
#define CHECKERR(ierr)
Definition: basic.h:18
hypar.h
Contains structure definition for hypar.
WENOFifthOrderCalculateWeightsCharYC
static int WENOFifthOrderCalculateWeightsCharYC(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:1272
_CRWENO_OPTIMAL_WEIGHT_2_
#define _CRWENO_OPTIMAL_WEIGHT_2_
Definition: interpolation.h:236
_CRWENO_OPTIMAL_WEIGHT_3_
#define _CRWENO_OPTIMAL_WEIGHT_3_
Definition: interpolation.h:238
WENOFifthOrderCalculateWeightsCharM
static int WENOFifthOrderCalculateWeightsCharM(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:930
basic.h
Some basic definitions and macros.
HyPar::ndims
int ndims
Definition: hypar.h:26
_WENOWeights_v_YC_
#define _WENOWeights_v_YC_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)
Definition: interpolation.h:479
arrayfunctions.h
Contains macros and function definitions for common array operations.
IERR
#define IERR
Definition: basic.h:16
_FIFTH_ORDER_CRWENO_
#define _FIFTH_ORDER_CRWENO_
Definition: interpolation.h:28
_ArrayProduct1D_
#define _ArrayProduct1D_(x, size, p)
Definition: arrayfunctions.h:377
MPIVariables
Structure of MPI-related variables.
Definition: mpivars_struct.h:24
HyPar
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
WENOFifthOrderCalculateWeightsZ
static int WENOFifthOrderCalculateWeightsZ(double *, double *, double *, int, void *, void *)
Definition: WENOFifthOrderCalculateWeights.c:434
matmult_native.h
Contains macros and function definitions for common matrix multiplication.