HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Interp1PrimFourthOrderCentralChar.c File Reference

Characteristic-based fourth order central scheme. More...

#include <stdio.h>
#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_   2
 

Functions

int Interp1PrimFourthOrderCentralChar (double *fI, double *fC, double *u, double *x, int upw, int dir, void *s, void *m, int uflag)
 4th order central reconstruction (characteristic-based) on a uniform grid More...
 

Detailed Description

Characteristic-based fourth order central scheme.

Author
Debojyoti Ghosh

Definition in file Interp1PrimFourthOrderCentralChar.c.

Macro Definition Documentation

#define _MINIMUM_GHOSTS_   2

Minimum number of ghost points required for this interpolation method.

Definition at line 25 of file Interp1PrimFourthOrderCentralChar.c.

Function Documentation

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

4th order central 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 4th order central 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 4th order central numerical approximation \(\hat{\bf f}_{j+1/2} \approx {\bf h}_{j+1/2}\) as:

\begin{equation} \hat{\alpha}^k_{j+1/2} = -\frac{1}{12} {\alpha}^k_{j-1} + \frac{7}{12} {\alpha}^k_{j} + \frac{7}{12} {\alpha}^k_{j+1} - \frac{1}{12} {\alpha}^k_{j+2}, \end{equation}

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 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.
  • Since this is a central scheme, the input argument upw has no effect.
  • 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}\)).
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 79 of file Interp1PrimFourthOrderCentralChar.c.

90 {
91  HyPar *solver = (HyPar*) s;
92  int i, k, v;
94 
95  int ghosts = solver->ghosts;
96  int ndims = solver->ndims;
97  int nvars = solver->nvars;
98  int *dim = solver->dim_local;
99 
100  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
101  dimension "dir" */
102  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
103  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
104  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
105  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
106 
107  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
108  double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
109 
110  static const double c1 = 7.0 / 12.0;
111  static const double c2 = -1.0 / 12.0;
112 
113 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
114  for (i=0; i<N_outer; i++) {
115  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
116  _ArrayCopy1D_(index_outer,indexC,ndims);
117  _ArrayCopy1D_(index_outer,indexI,ndims);
118 
119  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
120 
121  int pLL, pL, pR, pRR; /* 1D index of the left and right cells */
122  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pLL);
123  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pL );
124  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pR );
125  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pRR);
126  int p; /* 1D index of the interface */
127  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
128 
129  /* find averaged state at this interface */
130  IERR solver->AveragingFunction(uavg,&u[nvars*pL],&u[nvars*pR],solver->physics); CHECKERR(ierr);
131 
132  /* Get the left and right eigenvectors */
133  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
134  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
135 
136  /* For each characteristic field */
137  for (v = 0; v < nvars; v++) {
138 
139  /* calculate the characteristic flux components along this characteristic */
140  double fcLL = 0, fcL = 0, fcR = 0, fcRR = 0;
141  for (k = 0; k < nvars; k++) {
142  fcLL += L[v*nvars+k] * fC[pLL*nvars+k];
143  fcL += L[v*nvars+k] * fC[pL *nvars+k];
144  fcR += L[v*nvars+k] * fC[pR *nvars+k];
145  fcRR += L[v*nvars+k] * fC[pRR*nvars+k];
146  }
147 
148  /* first order upwind approximation of the characteristic flux */
149  fchar[v] = c2*fcLL + c1*fcL + c1*fcR + c2*fcRR;
150 
151  }
152 
153  /* calculate the interface u from the characteristic u */
154  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
155 
156  }
157  }
158 
159  return(0);
160 }
#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
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)
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
int ndims
Definition: hypar.h:26
#define IERR
Definition: basic.h:16
#define _ArrayProduct1D_(x, size, p)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _DECLARE_IERR_
Definition: basic.h:17