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

WENO5 Scheme (Component-wise application to vectors). More...

#include <stdlib.h>
#include <time.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <interpolation.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Macros

#define _MINIMUM_GHOSTS_   3
 

Functions

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

Detailed Description

WENO5 Scheme (Component-wise application to vectors).

Author
Debojyoti Ghosh

Definition in file Interp1PrimFifthOrderWENO.c.

Macro Definition Documentation

◆ _MINIMUM_GHOSTS_

#define _MINIMUM_GHOSTS_   3

Minimum number of ghost points required for this interpolation method.

Definition at line 25 of file Interp1PrimFifthOrderWENO.c.

Function Documentation

◆ Interp1PrimFifthOrderWENO()

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

5th order WENO 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 WENO 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 WENO numerical approximation \(\hat{\bf f}_{j+1/2} \approx {\bf h}_{j+1/2}\) as the convex combination of three 3rd order methods:

\begin{align} &\ \omega_1\ \times\ \left[ \hat{\bf f}_{j+1/2}^1 = \frac{1}{3} {\bf f}_{j-2} - \frac{7}{6} {\bf f}_{j-1} + \frac{11}{6} {\bf f}_j \right]\\ + &\ \omega_2\ \times\ \left[ \hat{\bf f}_{j+1/2}^2 = -\frac{1}{6} {\bf f}_{j-1} + \frac{5}{6} {\bf f}_j + \frac{1}{3} {\bf f}_{j+1} \right]\\ + &\ \omega_3\ \times\ \left[ \hat{\bf f}_{j+1/2}^3 = \frac{1}{3} {\bf f}_j + \frac{5}{6} {\bf f}_{j+1} - \frac{1}{6} {\bf f}_{j+2} \right]\\ \Rightarrow &\ \hat{\bf f}_{j+1/2} = \frac{\omega_1}{3} {\bf f}_{j-2} - \frac{1}{6}(7\omega_1+\omega_2){\bf f}_{j-1} + \frac{1}{6}(11\omega_1+5\omega_2+2\omega_3){\bf f}_j + \frac{1}{6}(2\omega_2+5\omega_3){\bf f}_{j+1} - \frac{\omega_3}{6}{\bf f}_{j+2}, \end{align}

where \(\omega_k; k=1,2,3\) are the nonlinear WENO weights computed in WENOFifthOrderCalculateWeights() (note that the \(\omega\) are different for each component of the vector \(\hat{\bf f}\)).

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:

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 74 of file Interp1PrimFifthOrderWENO.c.

85 {
86  HyPar *solver = (HyPar*) s;
87  WENOParameters *weno = (WENOParameters*) solver->interp;
88 
89  int ghosts = solver->ghosts;
90  int ndims = solver->ndims;
91  int nvars = solver->nvars;
92  int *dim = solver->dim_local;
93  int *stride= solver->stride_with_ghosts;
94 
95  /* define some constants */
96  static const double one_sixth = 1.0/6.0;
97 
98  double *ww1, *ww2, *ww3;
99  ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
100  ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
101  ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
102 
103  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
104  dimension "dir" */
105  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
106  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
107  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
108  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
109 
110 #if defined(CPU_STAT)
111  clock_t cpu_start, cpu_end;
112  cpu_start = clock();
113 #endif
114 
115  int i;
116 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
117  for (i=0; i<N_outer; i++) {
118  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
119  _ArrayCopy1D_(index_outer,indexC,ndims);
120  _ArrayCopy1D_(index_outer,indexI,ndims);
121  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
122  int qm1,qm2,qm3,qp1,qp2,p;
123  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
124  if (upw > 0) {
125  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
126  qm3 = qm1 - 2*stride[dir];
127  qm2 = qm1 - stride[dir];
128  qp1 = qm1 + stride[dir];
129  qp2 = qm1 + 2*stride[dir];
130  } else {
131  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
132  qm3 = qm1 + 2*stride[dir];
133  qm2 = qm1 + stride[dir];
134  qp1 = qm1 - stride[dir];
135  qp2 = qm1 - 2*stride[dir];
136  }
137 
138  /* Defining stencil points */
139  double *fm3, *fm2, *fm1, *fp1, *fp2;
140  fm3 = (fC+qm3*nvars);
141  fm2 = (fC+qm2*nvars);
142  fm1 = (fC+qm1*nvars);
143  fp1 = (fC+qp1*nvars);
144  fp2 = (fC+qp2*nvars);
145 
146  /* Candidate stencils and their optimal weights*/
147  double f1[nvars], f2[nvars], f3[nvars];
148  _ArrayAXBYCZ_(f1,(2*one_sixth),fm3,(-7*one_sixth) ,fm2,(11*one_sixth) ,fm1,nvars);
149  _ArrayAXBYCZ_(f2,(-one_sixth) ,fm2,(5*one_sixth) ,fm1,(2*one_sixth) ,fp1,nvars);
150  _ArrayAXBYCZ_(f3,(2*one_sixth),fm1,(5*one_sixth) ,fp1,(-one_sixth) ,fp2,nvars);
151 
152  /* calculate WENO weights */
153  double *w1,*w2,*w3;
154  w1 = (ww1+p*nvars);
155  w2 = (ww2+p*nvars);
156  w3 = (ww3+p*nvars);
157 
158  _ArrayMultiply3Add1D_((fI+p*nvars),w1,f1,w2,f2,w3,f3,nvars);
159  }
160  }
161 
162 #if defined(CPU_STAT)
163  cpu_end = clock();
164  printf("Interp1PrimFifthOrderWENO CPU time = %8.6lf dir = %d\n", (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC, dir);
165 #endif
166 
167  return(0);
168 }
int nvars
Definition: hypar.h:29
#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
#define _ArrayAXBYCZ_(w, a, x, b, y, c, z, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Structure of variables/parameters needed by the WENO-type scheme.
int * dim_local
Definition: hypar.h:37
void * interp
Definition: hypar.h:362
int ghosts
Definition: hypar.h:52
#define _ArrayCopy1D_(x, y, size)
#define _ArrayProduct1D_(x, size, p)
#define _ArrayMultiply3Add1D_(x, a, b, c, d, e, f, size)