HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Interp1PrimThirdOrderMUSCL.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 <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_ 2
25 
76  double *fI,
77  double *fC,
78  double *u,
79  double *x,
80  int upw,
81  int dir,
82  void *s,
83  void *m,
84  int uflag
85  )
86 {
87  HyPar *solver = (HyPar*) s;
88  MUSCLParameters *muscl = (MUSCLParameters*) solver->interp;
89 
90  int ghosts = solver->ghosts;
91  int ndims = solver->ndims;
92  int nvars = solver->nvars;
93  int *dim = solver->dim_local;
94 
95  /* define some constants */
96  double one_third = 1.0/3.0;
97  double one_sixth = 1.0/6.0;
98 
99  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
100  dimension "dir" */
101  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
102  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
103  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
104  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
105 
106  int i;
107  if (upw > 0) {
108 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
109  for (i=0; i<N_outer; i++) {
110  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
111  _ArrayCopy1D_(index_outer,indexC,ndims);
112  _ArrayCopy1D_(index_outer,indexI,ndims);
113  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
114  int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
115  int qm1,qm2,qp1,v;
116  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
117  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
118  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
119  for (v=0; v<nvars; v++) {
120  /* Defining stencil points */
121  double m2, m1, p1;
122  m2 = fC[qm2*nvars+v];
123  m1 = fC[qm1*nvars+v];
124  p1 = fC[qp1*nvars+v];
125 
126  double fdiff = p1 - m1;
127  double bdiff = m1 - m2;
128  double limit = (3*fdiff*bdiff + muscl->eps)
129  / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
130 
131  fI[p*nvars+v] = m1 + limit * (one_third*fdiff + one_sixth*bdiff);
132  }
133  }
134  }
135  } else {
136 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
137  for (i=0; i<N_outer; i++) {
138  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
139  _ArrayCopy1D_(index_outer,indexC,ndims);
140  _ArrayCopy1D_(index_outer,indexI,ndims);
141  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
142  int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
143  int qm1,qp1,qp2,v;
144  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
145  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
146  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
147  for (v=0; v<nvars; v++) {
148  /* Defining stencil points */
149  double m1, p1, p2;
150  m1 = fC[qm1*nvars+v];
151  p1 = fC[qp1*nvars+v];
152  p2 = fC[qp2*nvars+v];
153 
154  double fdiff = p2 - p1;
155  double bdiff = p1 - m1;
156  double limit = (3*fdiff*bdiff + muscl->eps)
157  / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
158 
159  fI[p*nvars+v] = p1 - limit * (one_third*fdiff + one_sixth*bdiff);
160  }
161  }
162  }
163  }
164 
165  return(0);
166 }
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.
int Interp1PrimThirdOrderMUSCL(double *, double *, double *, double *, int, int, void *, void *, int)
3rd order MUSCL scheme with Koren&#39;s limiter (component-wise) on a uniform grid
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Contains function definitions for common mathematical functions.
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
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 _ArrayProduct1D_(x, size, p)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23