HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Interp1PrimFourthOrderCentralChar.c
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <basic.h>
9 #include <arrayfunctions.h>
10 #include <matmult_native.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_ 2
26 
80  double *fI,
81  double *fC,
82  double *u,
83  double *x,
84  int upw,
85  int dir,
86  void *s,
87  void *m,
88  int uflag
89  )
90 {
91  HyPar *solver = (HyPar*) s;
92  int i, k, v;
94 
95  int ghosts = solver->ghosts;
96  int ndims = solver->ndims;
97  int nvars = solver->nvars;
98  int *dim = solver->dim_local;
99 
100  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
101  dimension "dir" */
102  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
103  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
104  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
105  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
106 
107  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
108  double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
109 
110  static const double c1 = 7.0 / 12.0;
111  static const double c2 = -1.0 / 12.0;
112 
113 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
114  for (i=0; i<N_outer; i++) {
115  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
116  _ArrayCopy1D_(index_outer,indexC,ndims);
117  _ArrayCopy1D_(index_outer,indexI,ndims);
118 
119  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
120 
121  int pLL, pL, pR, pRR; /* 1D index of the left and right cells */
122  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pLL);
123  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pL );
124  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pR );
125  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pRR);
126  int p; /* 1D index of the interface */
127  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
128 
129  /* find averaged state at this interface */
130  IERR solver->AveragingFunction(uavg,&u[nvars*pL],&u[nvars*pR],solver->physics); CHECKERR(ierr);
131 
132  /* Get the left and right eigenvectors */
133  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
134  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
135 
136  /* For each characteristic field */
137  for (v = 0; v < nvars; v++) {
138 
139  /* calculate the characteristic flux components along this characteristic */
140  double fcLL = 0, fcL = 0, fcR = 0, fcRR = 0;
141  for (k = 0; k < nvars; k++) {
142  fcLL += L[v*nvars+k] * fC[pLL*nvars+k];
143  fcL += L[v*nvars+k] * fC[pL *nvars+k];
144  fcR += L[v*nvars+k] * fC[pR *nvars+k];
145  fcRR += L[v*nvars+k] * fC[pRR*nvars+k];
146  }
147 
148  /* first order upwind approximation of the characteristic flux */
149  fchar[v] = c2*fcLL + c1*fcL + c1*fcR + c2*fcRR;
150 
151  }
152 
153  /* calculate the interface u from the characteristic u */
154  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
155 
156  }
157  }
158 
159  return(0);
160 }
int nvars
Definition: hypar.h:29
#define IERR
Definition: basic.h:16
MPI related function definitions.
#define CHECKERR(ierr)
Definition: basic.h:18
Contains function definitions for common mathematical functions.
int(* AveragingFunction)(double *, double *, double *, void *)
Definition: hypar.h:354
Some basic definitions and macros.
#define _ArrayIndexnD_(N, index, imax, i, ghost)
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Contains structure definition for hypar.
int Interp1PrimFourthOrderCentralChar(double *fI, double *fC, double *u, double *x, int upw, int dir, void *s, void *m, int uflag)
4th order central reconstruction (characteristic-based) on a uniform grid
#define _ArrayIndex1D_(N, imax, i, ghost, index)
int * dim_local
Definition: hypar.h:37
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
Contains macros and function definitions for common matrix multiplication.
Definitions for the functions computing the interpolated value of the primitive at the cell interface...
#define MatVecMult(N, y, A, x)
#define _DECLARE_IERR_
Definition: basic.h:17
#define _ArrayCopy1D_(x, y, size)
Contains macros and function definitions for common array operations.
int(* GetRightEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:359
#define _ArrayProduct1D_(x, size, p)
int(* GetLeftEigenvectors)(double *, double *, void *, int)
Definition: hypar.h:357