HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Interp1PrimFifthOrderWENOChar.c
Go to the documentation of this file.
1 
6 #include <stdlib.h>
7 #include <basic.h>
8 #include <arrayfunctions.h>
9 #include <matmult_native.h>
10 #include <mathfunctions.h>
11 #include <interpolation.h>
12 #include <mpivars.h>
13 #include <hypar.h>
14 
15 #ifdef with_omp
16 #include <omp.h>
17 #endif
18 
19 #undef _MINIMUM_GHOSTS_
20 
24 #define _MINIMUM_GHOSTS_ 3
25 
87  double *fI,
88  double *fC,
89  double *u,
90  double *x,
91  int upw,
92  int dir,
93  void *s,
94  void *m,
95  int uflag
96  )
97 {
98  HyPar *solver = (HyPar*) s;
99  WENOParameters *weno = (WENOParameters*) solver->interp;
100  int i, k, v;
102 
103  int ghosts = solver->ghosts;
104  int ndims = solver->ndims;
105  int nvars = solver->nvars;
106  int *dim = solver->dim_local;
107 
108  /* define some constants */
109  static const double one_sixth = 1.0/6.0;
110 
111  double *ww1, *ww2, *ww3;
112  ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
113  ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
114  ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
115 
116  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
117  dimension "dir" */
118  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
119  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
120  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
121  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
122 
123  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
124  double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
125 
126 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
127  for (i=0; i<N_outer; i++) {
128  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
129  _ArrayCopy1D_(index_outer,indexC,ndims);
130  _ArrayCopy1D_(index_outer,indexI,ndims);
131 
132  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
133 
134  /* 1D indices of the stencil grid points */
135  int qm1,qm2,qm3,qp1,qp2;
136  if (upw > 0) {
137  indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
138  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
139  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
140  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
141  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
142  } else {
143  indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
144  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
145  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
146  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
147  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
148  }
149 
150  int p; /* 1D index of the interface */
151  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
152 
153  /* find averaged state at this interface */
154  IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics); CHECKERR(ierr);
155 
156  /* Get the left and right eigenvectors */
157  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
158  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
159 
160  /* For each characteristic field */
161  for (v = 0; v < nvars; v++) {
162 
163  /* calculate the characteristic flux components along this characteristic */
164  double fm3, fm2, fm1, fp1, fp2;
165  fm3 = fm2 = fm1 = fp1 = fp2 = 0;
166  for (k = 0; k < nvars; k++) {
167  fm3 += L[v*nvars+k] * fC[qm3*nvars+k];
168  fm2 += L[v*nvars+k] * fC[qm2*nvars+k];
169  fm1 += L[v*nvars+k] * fC[qm1*nvars+k];
170  fp1 += L[v*nvars+k] * fC[qp1*nvars+k];
171  fp2 += L[v*nvars+k] * fC[qp2*nvars+k];
172  }
173 
174  /* Candidate stencils and their optimal weights*/
175  double f1, f2, f3;
176  f1 = (2*one_sixth)*fm3 - (7.0*one_sixth)*fm2 + (11.0*one_sixth)*fm1;
177  f2 = (-one_sixth)*fm2 + (5.0*one_sixth)*fm1 + (2*one_sixth)*fp1;
178  f3 = (2*one_sixth)*fm1 + (5*one_sixth)*fp1 - (one_sixth)*fp2;
179 
180  /* calculate WENO weights */
181  double w1,w2,w3;
182  w1 = *(ww1+p*nvars+v);
183  w2 = *(ww2+p*nvars+v);
184  w3 = *(ww3+p*nvars+v);
185 
186  /* fifth order WENO approximation of the characteristic flux */
187  fchar[v] = w1*f1 + w2*f2 + w3*f3;
188 
189  }
190 
191  /* calculate the interface u from the characteristic u */
192  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
193 
194  }
195  }
196 
197  return(0);
198 }
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)
void * physics
Definition: hypar.h:266
int(* GetRightEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:359
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)
int Interp1PrimFifthOrderWENOChar(double *, double *, double *, double *, int, int, void *, void *, int)
5th order WENO reconstruction (characteristic-based) on a uniform grid
int(* AveragingFunction)(double *, double *, double *, void *)
Definition: hypar.h:354
#define MatVecMult(N, y, A, x)
Contains function definitions for common mathematical functions.
int(* GetLeftEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:357
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
#define CHECKERR(ierr)
Definition: basic.h:18
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 IERR
Definition: basic.h:16
#define _ArrayProduct1D_(x, size, p)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Contains macros and function definitions for common matrix multiplication.
#define _DECLARE_IERR_
Definition: basic.h:17