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

Characteristic-based WENO5 Scheme. More...

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <matmult_native.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 Interp1PrimFifthOrderWENOChar (double *fI, double *fC, double *u, double *x, int upw, int dir, void *s, void *m, int uflag)
 5th order WENO reconstruction (characteristic-based) on a uniform grid More...
 

Detailed Description

Characteristic-based WENO5 Scheme.

Author
Debojyoti Ghosh

Definition in file Interp1PrimFifthOrderWENOChar.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 Interp1PrimFifthOrderWENOChar.c.

Function Documentation

◆ Interp1PrimFifthOrderWENOChar()

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

5th order WENO reconstruction (characteristic-based) 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{\alpha}^k_{j+1/2} = \frac{1}{3} {\alpha}^k_{j-2} - \frac{7}{6} {\alpha}^k_{j-1} + \frac{11}{6} {\alpha}^k_j \right]\\ + &\ \omega_2\ \times\ \left[ \hat{\alpha}^k_{j+1/2} = -\frac{1}{6} {\alpha}^k_{j-1} + \frac{5}{6} {\alpha}^k_j + \frac{1}{3} {\alpha}^k_{j+1} \right]\\ + &\ \omega_3\ \times\ \left[ \hat{\alpha}^k_{j+1/2} = \frac{1}{3} {\alpha}^k_j + \frac{5}{6} {\alpha}^k_{j+1} - \frac{1}{6} {\alpha}^k_{j+2} \right]\\ \Rightarrow &\ \hat{\alpha}^k_{j+1/2} = \frac{\omega_1}{3} {\alpha}^k_{j-2} - \frac{1}{6}(7\omega_1+\omega_2){\alpha}^k_{j-1} + \frac{1}{6}(11\omega_1+5\omega_2+2\omega_3){\alpha}^k_j + \frac{1}{6}(2\omega_2+5\omega_3){\alpha}^k_{j+1} - \frac{\omega_3}{6}{\alpha}^k_{j+2}, \end{align}

where

\begin{equation} \alpha^k = {\bf l}_k \cdot {\bf f},\ k=1,\cdots,n \end{equation}

is the \(k\)-th characteristic quantity, and \({\bf l}_k\) is the \(k\)-th left eigenvector, \({\bf r}_k\) is the \(k\)-th right eigenvector, and \(n\) is HyPar::nvars. The nonlinear weights \(\omega_k; k=1,2,3\) are the WENO weights computed in WENOFifthOrderCalculateWeightsChar(). The final interpolated function is computed from the interpolated characteristic quantities as:

\begin{equation} \hat{\bf f}_{j+1/2} = \sum_{k=1}^n \alpha^k_{j+1/2} {\bf r}_k \end{equation}

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 left and right eigenvectors are computed at an averaged quantity at j+1/2. Thus, this function requires functions to compute the average state, and the left and right eigenvectors. These are provided by the physical model through

    If these functions are not provided by the physical model, then a characteristic-based interpolation cannot be used.

  • 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 86 of file Interp1PrimFifthOrderWENOChar.c.

97 {
98  HyPar *solver = (HyPar*) s;
99  WENOParameters *weno = (WENOParameters*) solver->interp;
100  int i, k, v;
102 
103  int ghosts = solver->ghosts;
104  int ndims = solver->ndims;
105  int nvars = solver->nvars;
106  int *dim = solver->dim_local;
107 
108  /* define some constants */
109  static const double one_sixth = 1.0/6.0;
110 
111  double *ww1, *ww2, *ww3;
112  ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
113  ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
114  ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
115 
116  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
117  dimension "dir" */
118  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
119  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
120  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
121  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
122 
123  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
124  double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
125 
126 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
127  for (i=0; i<N_outer; i++) {
128  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
129  _ArrayCopy1D_(index_outer,indexC,ndims);
130  _ArrayCopy1D_(index_outer,indexI,ndims);
131 
132  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
133 
134  /* 1D indices of the stencil grid points */
135  int qm1,qm2,qm3,qp1,qp2;
136  if (upw > 0) {
137  indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
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  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
142  } else {
143  indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
144  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
145  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
146  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
147  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
148  }
149 
150  int p; /* 1D index of the interface */
151  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
152 
153  /* find averaged state at this interface */
154  IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics); CHECKERR(ierr);
155 
156  /* Get the left and right eigenvectors */
157  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
158  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
159 
160  /* For each characteristic field */
161  for (v = 0; v < nvars; v++) {
162 
163  /* calculate the characteristic flux components along this characteristic */
164  double fm3, fm2, fm1, fp1, fp2;
165  fm3 = fm2 = fm1 = fp1 = fp2 = 0;
166  for (k = 0; k < nvars; k++) {
167  fm3 += L[v*nvars+k] * fC[qm3*nvars+k];
168  fm2 += L[v*nvars+k] * fC[qm2*nvars+k];
169  fm1 += L[v*nvars+k] * fC[qm1*nvars+k];
170  fp1 += L[v*nvars+k] * fC[qp1*nvars+k];
171  fp2 += L[v*nvars+k] * fC[qp2*nvars+k];
172  }
173 
174  /* Candidate stencils and their optimal weights*/
175  double f1, f2, f3;
176  f1 = (2*one_sixth)*fm3 - (7.0*one_sixth)*fm2 + (11.0*one_sixth)*fm1;
177  f2 = (-one_sixth)*fm2 + (5.0*one_sixth)*fm1 + (2*one_sixth)*fp1;
178  f3 = (2*one_sixth)*fm1 + (5*one_sixth)*fp1 - (one_sixth)*fp2;
179 
180  /* calculate WENO weights */
181  double w1,w2,w3;
182  w1 = *(ww1+p*nvars+v);
183  w2 = *(ww2+p*nvars+v);
184  w3 = *(ww3+p*nvars+v);
185 
186  /* fifth order WENO approximation of the characteristic flux */
187  fchar[v] = w1*f1 + w2*f2 + w3*f3;
188 
189  }
190 
191  /* calculate the interface u from the characteristic u */
192  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
193 
194  }
195  }
196 
197  return(0);
198 }
int nvars
Definition: hypar.h:29
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
int(* AveragingFunction)(double *, double *, double *, void *)
Definition: hypar.h:354
#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
#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
void * physics
Definition: hypar.h:266
int ghosts
Definition: hypar.h:52
#define MatVecMult(N, y, A, x)
#define _DECLARE_IERR_
Definition: basic.h:17
#define _ArrayCopy1D_(x, y, size)
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