HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Interp1PrimThirdOrderMUSCLChar.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_ 3
26 
92  double *fI,
93  double *fC,
94  double *u,
95  double *x,
96  int upw,
97  int dir,
98  void *s,
99  void *m,
100  int uflag
101  )
102 {
103  HyPar *solver = (HyPar*) s;
104  MUSCLParameters *muscl = (MUSCLParameters*) solver->interp;
105  int i, k, v;
107 
108  int ghosts = solver->ghosts;
109  int ndims = solver->ndims;
110  int nvars = solver->nvars;
111  int *dim = solver->dim_local;
112 
113  /* define some constants */
114  double one_third = 1.0/3.0;
115  double one_sixth = 1.0/6.0;
116 
117  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
118  dimension "dir" */
119  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
120  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
121  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
122  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
123 
124  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
125  double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
126 
127  if (upw > 0) {
128 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
129  for (i=0; i<N_outer; i++) {
130  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
131  _ArrayCopy1D_(index_outer,indexC,ndims);
132  _ArrayCopy1D_(index_outer,indexI,ndims);
133 
134  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
135 
136  /* 1D indices of the stencil grid points */
137  int qm1,qm2,qp1;
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 
142  int p; /* 1D index of the interface */
143  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
144 
145  /* find averaged state at this interface */
146  IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
147  CHECKERR(ierr);
148 
149  /* Get the left and right eigenvectors */
150  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
151  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
152 
153  /* For each characteristic field */
154  for (v = 0; v < nvars; v++) {
155  /* calculate the characteristic flux components along this characteristic */
156  double m2, m1, p1;
157  m2 = m1 = p1 = 0;
158  for (k = 0; k < nvars; k++) {
159  m2 += L[v*nvars+k] * fC[qm2*nvars+k];
160  m1 += L[v*nvars+k] * fC[qm1*nvars+k];
161  p1 += L[v*nvars+k] * fC[qp1*nvars+k];
162  }
163  double fdiff = p1 - m1;
164  double bdiff = m1 - m2;
165  double limit = (3*fdiff*bdiff + muscl->eps)
166  / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
167  fchar[v] = m1 + limit * (one_third*fdiff + one_sixth*bdiff);
168  }
169 
170  /* calculate the interface u from the characteristic u */
171  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
172  }
173  }
174  } else {
175 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
176  for (i=0; i<N_outer; i++) {
177  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
178  _ArrayCopy1D_(index_outer,indexC,ndims);
179  _ArrayCopy1D_(index_outer,indexI,ndims);
180 
181  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
182 
183  /* 1D indices of the stencil grid points */
184  int qm1,qp1,qp2;
185  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
186  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
187  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
188 
189  int p; /* 1D index of the interface */
190  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
191 
192  /* find averaged state at this interface */
193  IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
194  CHECKERR(ierr);
195 
196  /* Get the left and right eigenvectors */
197  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
198  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
199 
200  /* For each characteristic field */
201  for (v = 0; v < nvars; v++) {
202  /* calculate the characteristic flux components along this characteristic */
203  double m1, p1, p2;
204  m1 = p1 = p2 = 0;
205  for (k = 0; k < nvars; k++) {
206  m1 += L[v*nvars+k] * fC[qm1*nvars+k];
207  p1 += L[v*nvars+k] * fC[qp1*nvars+k];
208  p2 += L[v*nvars+k] * fC[qp2*nvars+k];
209  }
210  double fdiff = p2 - p1;
211  double bdiff = p1 - m1;
212  double limit = (3*fdiff*bdiff + muscl->eps)
213  / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
214  fchar[v] = p1 - limit * (one_third*fdiff + one_sixth*bdiff);
215  }
216 
217  /* calculate the interface u from the characteristic u */
218  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
219  }
220  }
221  }
222 
223  return(0);
224 }
void * interp
Definition: hypar.h:362
int Interp1PrimThirdOrderMUSCLChar(double *, double *, double *, double *, int, int, void *, void *, int)
3rd order MUSCL scheme with Koren&#39;s limiter (characteristic-based) on a uniform grid ...
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.
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
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.
Structure of variables/parameters needed by the MUSCL scheme.
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