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

2nd order MUSCL scheme (component-wise application to vectors) More...

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

Functions

int Interp1PrimSecondOrderMUSCL (double *fI, double *fC, double *u, double *x, int upw, int dir, void *s, void *m, int uflag)
 2nd order MUSCL scheme (component-wise) on a uniform grid More...
 

Detailed Description

2nd order MUSCL scheme (component-wise application to vectors)

Author
Debojyoti Ghosh

Definition in file Interp1PrimSecondOrderMUSCL.c.

Macro Definition Documentation

#define _MINIMUM_GHOSTS_   2

Minimum number of ghost points required for this interpolation method.

Definition at line 24 of file Interp1PrimSecondOrderMUSCL.c.

Function Documentation

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

2nd order MUSCL scheme (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 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:

\begin{equation} \hat{\bf f}_{j+1/2} = {\bf f}_{j} + \frac{1}{2} \phi\left(r_j\right) \left[{\bf f}_{j+1}-{\bf f}_{j}\right] \end{equation}

where

\begin{equation} r_j = \left( f_j - f_{j-1} \right) / \left( f_{j+1} - f_j \right) \end{equation}

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

Implementation Notes:

  • The scalar interpolation method is applied to the vector function in a component-wise manner.
  • 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 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 76 of file Interp1PrimSecondOrderMUSCL.c.

87 {
88  HyPar *solver = (HyPar*) s;
89  MUSCLParameters *muscl = (MUSCLParameters*) solver->interp;
90 
91  int ghosts = solver->ghosts;
92  int ndims = solver->ndims;
93  int nvars = solver->nvars;
94  int *dim = solver->dim_local;
95 
96  /* define some constants */
97  double half = 1.0/2.0;
98 
99  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
100  dimension "dir" */
101  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
102  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
103  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
104  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
105 
106  int i;
107  if (upw > 0) {
108 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
109  for (i=0; i<N_outer; i++) {
110  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
111  _ArrayCopy1D_(index_outer,indexC,ndims);
112  _ArrayCopy1D_(index_outer,indexI,ndims);
113  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
114  int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
115  int qm2,qm1,qp1,v;
116  indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
117  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
118  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
119  for (v=0; v<nvars; v++) {
120  /* Defining stencil points */
121  double m2, m1, p1;
122  m2 = fC[qm2*nvars+v];
123  m1 = fC[qm1*nvars+v];
124  p1 = fC[qp1*nvars+v];
125 
126  double slope_ratio = (m1 - m2) / ((p1 - m1) + 1e-40);
127  double phi = muscl->LimiterFunction(slope_ratio);
128  fI[p*nvars+v] = m1 + 0.5 * phi * (p1-m1);
129  }
130  }
131  }
132  } else {
133 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
134  for (i=0; i<N_outer; i++) {
135  _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
136  _ArrayCopy1D_(index_outer,indexC,ndims);
137  _ArrayCopy1D_(index_outer,indexI,ndims);
138  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
139  int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
140  int qm1,qp1,qp2,v;
141  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
142  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
143  indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
144  for (v=0; v<nvars; v++) {
145  /* Defining stencil points */
146  double m1, p1, p2;
147  m1 = fC[qm1*nvars+v];
148  p1 = fC[qp1*nvars+v];
149  p2 = fC[qp2*nvars+v];
150 
151  double slope_ratio = (p1 - m1) / ((p2 - p1) + 1e-40);
152  double phi = muscl->LimiterFunction(slope_ratio);
153  fI[p*nvars+v] = p1 + 0.5 * phi * (p1-p2);
154  }
155  }
156  }
157  }
158 
159  return(0);
160 }
void * interp
Definition: hypar.h:362
#define _ArrayIndexnD_(N, index, imax, i, ghost)
int * dim_local
Definition: hypar.h:37
double(* LimiterFunction)(double)
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
Structure of variables/parameters needed by the MUSCL scheme.
int ndims
Definition: hypar.h:26
#define _ArrayProduct1D_(x, size, p)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23