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

Characteristic-based 2nd-order MUSCL 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 Interp1PrimSecondOrderMUSCLChar (double *fI, double *fC, double *u, double *x, int upw, int dir, void *s, void *m, int uflag)
 2nd order MUSCL scheme (characteristic-based) on a uniform grid More...
 

Detailed Description

Characteristic-based 2nd-order MUSCL scheme.

Author
Debojyoti Ghosh

Definition in file Interp1PrimSecondOrderMUSCLChar.c.

Macro Definition Documentation

◆ _MINIMUM_GHOSTS_

#define _MINIMUM_GHOSTS_   2

Minimum number of ghost points required for this interpolation method.

Definition at line 25 of file Interp1PrimSecondOrderMUSCLChar.c.

Function Documentation

◆ Interp1PrimSecondOrderMUSCLChar()

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

2nd order MUSCL scheme (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 2nd order MUSCL 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 numerical approximation \(\hat{\bf f}_{j+1/2} \approx {\bf h}_{j+1/2}\) as: using the 2nd order MUSCL scheme as follows:

\begin{equation} \hat{\alpha}^k_{j+1/2} = {\alpha}^k_{j} + \frac{1}{2} \phi\left(r_j\right) \left[{\alpha}^k_{j+1}-{\alpha}^k_{j}\right] \end{equation}

where

\begin{equation} r_j = \left( \alpha^k_j - \alpha^k_{j-1} \right) / \left( \alpha^k_{j+1} - \alpha^k_j \right), \end{equation}

\(\phi\left(r\right)\) is a limiter (minmod, mc, generalized minmod, etc.), and

\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.
  • 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:

  • van Leer, B. (1979), Towards the Ultimate Conservative Difference Scheme, V. A Second Order Sequel to Godunov's Method, J. Com. Phys., 32, 101–136.
  • Nessyahu, H. and E. Tadmor (1990), Non-oscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87, 408–463.
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 92 of file Interp1PrimSecondOrderMUSCLChar.c.

103 {
104  HyPar *solver = (HyPar*) s;
105  MUSCLParameters *muscl = (MUSCLParameters*) solver->interp;
106  int i, k, v;
108 
109  int ghosts = solver->ghosts;
110  int ndims = solver->ndims;
111  int nvars = solver->nvars;
112  int *dim = solver->dim_local;
113 
114  /* define some constants */
115  double one_third = 1.0/3.0;
116  double one_sixth = 1.0/6.0;
117 
118  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
119  dimension "dir" */
120  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
121  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
122  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
123  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
124 
125  /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
126  double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
127 
128  if (upw > 0) {
129 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
130  for (i=0; i<N_outer; i++) {
131  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
132  _ArrayCopy1D_(index_outer,indexC,ndims);
133  _ArrayCopy1D_(index_outer,indexI,ndims);
134 
135  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
136 
137  /* 1D indices of the stencil grid points */
138  int qm1,qm2,qp1;
139  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
140  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
141  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
142 
143  int p; /* 1D index of the interface */
144  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
145 
146  /* find averaged state at this interface */
147  IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
148  CHECKERR(ierr);
149 
150  /* Get the left and right eigenvectors */
151  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
152  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
153 
154  /* For each characteristic field */
155  for (v = 0; v < nvars; v++) {
156  /* calculate the characteristic flux components along this characteristic */
157  double m2, m1, p1;
158  m2 = m1 = p1 = 0;
159  for (k = 0; k < nvars; k++) {
160  m2 += L[v*nvars+k] * fC[qm2*nvars+k];
161  m1 += L[v*nvars+k] * fC[qm1*nvars+k];
162  p1 += L[v*nvars+k] * fC[qp1*nvars+k];
163  }
164 
165  double slope_ratio = (m1 - m2) / ((p1 - m1) + 1e-40);
166  double phi = muscl->LimiterFunction(slope_ratio);
167  fchar[v] = m1 + 0.5 * phi * (p1-m1);
168  }
169 
170  /* calculate the interface u from the characteristic u */
171  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
172  }
173  }
174  } else {
175 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
176  for (i=0; i<N_outer; i++) {
177  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
178  _ArrayCopy1D_(index_outer,indexC,ndims);
179  _ArrayCopy1D_(index_outer,indexI,ndims);
180 
181  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
182 
183  /* 1D indices of the stencil grid points */
184  int qm1,qp1,qp2;
185  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
186  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
187  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
188 
189  int p; /* 1D index of the interface */
190  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
191 
192  /* find averaged state at this interface */
193  IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
194  CHECKERR(ierr);
195 
196  /* Get the left and right eigenvectors */
197  IERR solver->GetLeftEigenvectors (uavg,L,solver->physics,dir); CHECKERR(ierr);
198  IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
199 
200  /* For each characteristic field */
201  for (v = 0; v < nvars; v++) {
202  /* calculate the characteristic flux components along this characteristic */
203  double m1, p1, p2;
204  m1 = p1 = p2 = 0;
205  for (k = 0; k < nvars; k++) {
206  m1 += L[v*nvars+k] * fC[qm1*nvars+k];
207  p1 += L[v*nvars+k] * fC[qp1*nvars+k];
208  p2 += L[v*nvars+k] * fC[qp2*nvars+k];
209  }
210 
211  double slope_ratio = (p1 - m1) / ((p2 - p1) + 1e-40);
212  double phi = muscl->LimiterFunction(slope_ratio);
213  fchar[v] = p1 + 0.5 * phi * (p1-p2);
214  }
215 
216  /* calculate the interface u from the characteristic u */
217  IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
218  }
219  }
220  }
221 
222  return(0);
223 }
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
double(* LimiterFunction)(double)
#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)
int * dim_local
Definition: hypar.h:37
void * interp
Definition: hypar.h:362
void * physics
Definition: hypar.h:266
Structure of variables/parameters needed by the MUSCL scheme.
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