HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Interp1PrimFifthOrderCompactUpwind.c File Reference

5th order compact upwind scheme (component-wise application to vectors). More...

#include <stdio.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <interpolation.h>
#include <tridiagLU.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Macros

#define _MINIMUM_GHOSTS_   3
 

Functions

int Interp1PrimFifthOrderCompactUpwind (double *fI, double *fC, double *u, double *x, int upw, int dir, void *s, void *m, int uflag)
 5th order compact upwind reconstruction (component-wise) on a uniform grid More...
 

Detailed Description

5th order compact upwind scheme (component-wise application to vectors).

Author
Debojyoti Ghosh

Definition in file Interp1PrimFifthOrderCompactUpwind.c.

Macro Definition Documentation

◆ _MINIMUM_GHOSTS_

#define _MINIMUM_GHOSTS_   3

Minimum number of ghost points required for this interpolation method.

Definition at line 24 of file Interp1PrimFifthOrderCompactUpwind.c.

Function Documentation

◆ Interp1PrimFifthOrderCompactUpwind()

int Interp1PrimFifthOrderCompactUpwind ( double *  fI,
double *  fC,
double *  u,
double *  x,
int  upw,
int  dir,
void *  s,
void *  m,
int  uflag 
)

5th order compact upwind reconstruction (component-wise) on a uniform grid

Computes the interpolated values of the first primitive of a function \({\bf f}\left({\bf u}\right)\) at the interfaces from the cell-centered values of the function using the fifth order compact upwind scheme on a uniform grid. The first primitive is defined as a function \({\bf h}\left({\bf u}\right)\) that satisfies:

\begin{equation} {\bf f}\left({\bf u}\left(x\right)\right) = \frac{1}{\Delta x} \int_{x-\Delta x/2}^{x+\Delta x/2} {\bf h}\left({\bf u}\left(\zeta\right)\right)d\zeta, \end{equation}

where \(x\) is the spatial coordinate along the dimension of the interpolation. This function computes the 5th order compact upwind numerical approximation \(\hat{\bf f}_{j+1/2} \approx {\bf h}_{j+1/2}\) as:

\begin{align} \frac{3}{10}\hat{\bf f}_{j-1/2} + \frac{6}{10}\hat{\bf f}_{j+1/2} + \frac{1}{10}\hat{\bf f}_{j+3/2} = \frac{1}{30}{\bf f}_{j-1} + \frac{19}{30}{\bf f}_j + \frac{1}{3}{\bf f}_{j+1}. \end{align}

The resulting tridiagonal system is solved using tridiagLU() (see also TridiagLU, tridiagLU.h).

Implementation Notes:

  • This method assumes a uniform grid in the spatial dimension corresponding to the interpolation.
  • The method described above corresponds to a left-biased interpolation. The corresponding right-biased interpolation can be obtained by reflecting the equations about interface j+1/2.
  • The scalar interpolation method is applied to the vector function in a component-wise manner.
  • The function computes the interpolant for the entire grid in one call. It loops over all the grid lines along the interpolation direction and carries out the 1D interpolation along these grid lines.
  • Location of cell-centers and cell interfaces along the spatial dimension of the interpolation is shown in the following figure:
    chap1_1Ddomain.png

Function arguments:

Argument Type Explanation ------—
fI double* Array to hold the computed interpolant at the grid interfaces. This array must have the same layout as the solution, but with no ghost points. Its size should be the same as u in all dimensions, except dir (the dimension along which to interpolate) along which it should be larger by 1 (number of interfaces is 1 more than the number of interior cell centers).
fC double* Array with the cell-centered values of the flux function \({\bf f}\left({\bf u}\right)\). This array must have the same layout and size as the solution, with ghost points.
u double* The solution array \({\bf u}\) (with ghost points). If the interpolation is characteristic based, this is needed to compute the eigendecomposition. For a multidimensional problem, the layout is as follows: u is a contiguous 1D array of size (nvars*dim[0]*dim[1]*...*dim[D-1]) corresponding to the multi-dimensional solution, with the following ordering - nvars, dim[0], dim[1], ..., dim[D-1], where nvars is the number of solution components (HyPar::nvars), dim is the local size (HyPar::dim_local), D is the number of spatial dimensions.
x double* The grid array (with ghost points). This is used only by non-uniform-grid interpolation methods. For multidimensional problems, the layout is as follows: x is a contiguous 1D array of size (dim[0]+dim[1]+...+dim[D-1]), with the spatial coordinates along dim[0] stored from 0,...,dim[0]-1, the spatial coordinates along dim[1] stored along dim[0],...,dim[0]+dim[1]-1, and so forth.
upw int Upwinding direction: if positive, a left-biased interpolant will be computed; if negative, a right-biased interpolant will be computed. If the interpolation method is central, then this has no effect.
dir int Spatial dimension along which to interpolate (eg: 0 for 1D; 0 or 1 for 2D; 0,1 or 2 for 3D)
s void* Solver object of type HyPar: the following variables are needed - HyPar::ghosts, HyPar::ndims, HyPar::nvars, HyPar::dim_local.
m void* MPI object of type MPIVariables: this is needed only by compact interpolation method that need to solve a global implicit system across MPI ranks.
uflag int A flag indicating if the function being interpolated \({\bf f}\) is the solution itself \({\bf u}\) (if 1, \({\bf f}\left({\bf u}\right) \equiv {\bf u}\)).

Reference:

  • Ghosh, D., Baeder, J. D., Compact Reconstruction Schemes with Weighted ENO Limiting for Hyperbolic Conservation Laws, SIAM Journal on Scientific Computing, 34 (3), 2012, A1678–A1706, http://dx.doi.org/10.1137/110857659
Parameters
fIArray of interpolated function values at the interfaces
fCArray of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
uArray of cell-centered values of the solution \({\bf u}\)
xGrid coordinates
upwUpwind direction (left or right biased)
dirSpatial dimension along which to interpolation
sObject of type HyPar containing solver-related variables
mObject of type MPIVariables containing MPI-related variables
uflagFlag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 72 of file Interp1PrimFifthOrderCompactUpwind.c.

83 {
84  HyPar *solver = (HyPar*) s;
85  MPIVariables *mpi = (MPIVariables*) m;
86  CompactScheme *compact= (CompactScheme*) solver->compact;
87  TridiagLU *lu = (TridiagLU*) solver->lusolver;
88  int sys,Nsys,d,v;
90 
91  int ghosts = solver->ghosts;
92  int ndims = solver->ndims;
93  int nvars = solver->nvars;
94  int *dim = solver->dim_local;
95  int *stride= solver->stride_with_ghosts;
96 
97  /* define some constants */
98  static const double one_third = 1.0/3.0,
99  thirteen_by_sixty = 13.0/60.0,
100  fortyseven_by_sixty = 47.0/60.0,
101  twentyseven_by_sixty = 27.0/60.0,
102  one_by_twenty = 1.0/20.0,
103  one_by_thirty = 1.0/30.0,
104  nineteen_by_thirty = 19.0/30.0,
105  three_by_ten = 3.0/10.0,
106  six_by_ten = 6.0/10.0,
107  one_by_ten = 1.0/10.0;
108 
109  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
110  dimension "dir" */
111  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
112  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
113  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
114  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
115 
116  /* calculate total number of tridiagonal systems to solve */
117  _ArrayProduct1D_(bounds_outer,ndims,Nsys); Nsys *= nvars;
118 
119  /* Allocate arrays for tridiagonal system */
120  double *A = compact->A;
121  double *B = compact->B;
122  double *C = compact->C;
123  double *R = compact->R;
124 
125 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
126  for (sys=0; sys < N_outer; sys++) {
127  _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
128  _ArrayCopy1D_(index_outer,indexC,ndims);
129  _ArrayCopy1D_(index_outer,indexI,ndims);
130  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
131  int qm1,qm2,qm3,qp1,qp2,p;
132  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
133  if (upw > 0) {
134  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
135  qm3 = qm1 - 2*stride[dir];
136  qm2 = qm1 - stride[dir];
137  qp1 = qm1 + stride[dir];
138  qp2 = qm1 + 2*stride[dir];
139  } else {
140  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
141  qm3 = qm1 + 2*stride[dir];
142  qm2 = qm1 + stride[dir];
143  qp1 = qm1 - stride[dir];
144  qp2 = qm1 - 2*stride[dir];
145  }
146 
147  /* Defining stencil points */
148  double *fm3, *fm2, *fm1, *fp1, *fp2;
149  fm3 = fC+qm3*nvars;
150  fm2 = fC+qm2*nvars;
151  fm1 = fC+qm1*nvars;
152  fp1 = fC+qp1*nvars;
153  fp2 = fC+qp2*nvars;
154 
155  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
156  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
157 
158  /* Use 5th order upwind at the physical boundaries */
159  _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
160  _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,1.0)
161  _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
162  for (v=0; v<nvars; v++) {
163  (R+Nsys*indexI[dir]+sys*nvars)[v] = one_by_thirty * fm3[v]
164  - thirteen_by_sixty * fm2[v]
165  + fortyseven_by_sixty * fm1[v]
166  + twentyseven_by_sixty * fp1[v]
167  - one_by_twenty * fp2[v];
168  }
169 
170  } else {
171 
172  /* Use 5th order upwind at the physical boundaries */
173  if (upw > 0) {
174  _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,three_by_ten);
175  _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,six_by_ten );
176  _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,one_by_ten );
177  } else {
178  _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,three_by_ten);
179  _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,six_by_ten );
180  _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,one_by_ten );
181  }
182  for (v=0; v<nvars; v++) {
183  (R+Nsys*indexI[dir]+sys*nvars)[v] = one_by_thirty * fm2[v]
184  + nineteen_by_thirty * fm1[v]
185  + one_third * fp1[v];
186  }
187  }
188  }
189  }
190 
191 #ifdef serial
192 
193  /* Solve the tridiagonal system */
194  IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,NULL); CHECKERR(ierr);
195 
196 #else
197 
198  /* Solve the tridiagonal system */
199  /* all processes except the last will solve without the last interface to avoid overlap */
200  if (mpi->ip[dir] != mpi->iproc[dir]-1) { IERR tridiagLU(A,B,C,R,dim[dir] ,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
201  else { IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
202 
203  /* Now get the solution to the last interface from the next proc */
204  double *sendbuf = compact->sendbuf;
205  double *recvbuf = compact->recvbuf;
206  MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
207  if (mpi->ip[dir]) for (d=0; d<Nsys; d++) sendbuf[d] = R[d];
208  if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
209  if (mpi->ip[dir]) MPI_Isend(sendbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
210  MPI_Status status_arr[2];
211  MPI_Waitall(2,&req[0],status_arr);
212  if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys; d++) R[d+Nsys*dim[dir]] = recvbuf[d];
213 
214 #endif
215 
216  /* save the solution to fI */
217 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
218  for (sys=0; sys < N_outer; sys++) {
219  _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
220  _ArrayCopy1D_(index_outer,indexI,ndims);
221  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
222  int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
223  _ArrayCopy1D_((R+sys*nvars+Nsys*indexI[dir]),(fI+nvars*p),nvars);
224  }
225  }
226 
227  return(0);
228 }
int nvars
Definition: hypar.h:29
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
#define _ArrayIndexnD_(N, index, imax, i, ghost)
int * stride_with_ghosts
Definition: hypar.h:414
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
MPI_Comm * comm
#define _ArrayIndex1D_(N, imax, i, ghost, index)
void * lusolver
Definition: hypar.h:368
double * recvbuf
Structure of variables/parameters needed by the compact schemes.
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
int tridiagLU(double *, double *, double *, double *, int, int, void *, void *)
Definition: tridiagLU.c:83
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
#define _DECLARE_IERR_
Definition: basic.h:17
#define _ArrayCopy1D_(x, y, size)
#define _ArrayProduct1D_(x, size, p)
void * compact
Definition: hypar.h:364
double * sendbuf