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

Compute the hyperbolic term of the governing equations. More...

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

Go to the source code of this file.

Functions

static int ReconstructHyperbolic (double *, double *, double *, double *, int, void *, void *, double, int, int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
 
static int DefaultUpwinding (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int HyperbolicFunction (double *hyp, double *u, void *s, void *m, double t, int LimFlag, int(*FluxFunction)(double *, double *, int, void *, double), int(*UpwindFunction)(double *, double *, double *, double *, double *, double *, int, void *, double))
 

Detailed Description

Compute the hyperbolic term of the governing equations.

Author
Debojyoti Ghosh

Definition in file HyperbolicFunction.c.

Function Documentation

int ReconstructHyperbolic ( double *  fluxI,
double *  fluxC,
double *  u,
double *  x,
int  dir,
void *  s,
void *  m,
double  t,
int  LimFlag,
int(*)(double *, double *, double *, double *, double *, double *, int, void *, double)  UpwindFunction 
)
static

This function computes the numerical flux \(\hat{\bf f}_{j+1/2}\) at the interface from the cell-centered flux function \({\bf f}_j\). This happens in two steps:-

Interpolation: High-order accurate approximations to the flux at the interface \(j+1/2\) are computed from the cell-centered flux with left- and right-biased interpolation methods. This is done by the HyPar::InterpolateInterfacesHyp function. This can be expressed as follows:

\begin{align} \hat{\bf f}^L_{j+1/2} &= \mathcal{I}\left({\bf f}_j,+1\right), \\ \hat{\bf f}^R_{j+1/2} &= \mathcal{I}\left({\bf f}_j,-1\right), \end{align}

where the \(\pm 1\) indicates the interpolation bias, and \(\mathcal{I}\) is the interpolation operator pointed to by HyPar::InterpolateInterfacesHyp (see src/InterpolationFunctions for all the available operators). The interface values of the solution are similarly computed:

\begin{align} \hat{\bf u}^L_{j+1/2} &= \mathcal{I}\left({\bf u}_j,+1\right), \\ \hat{\bf u}^R_{j+1/2} &= \mathcal{I}\left({\bf u}_j,-1\right), \end{align}

The specific choice of \(\mathcal{I}\) is set based on HyPar::spatial_scheme_hyp.

Upwinding: The final flux at the interface is computed as

\begin{equation} \hat{\bf f}_{j+1/2} = \mathcal{U}\left( \hat{\bf f}^L_{j+1/2}, \hat{\bf f}^R_{j+1/2}, \hat{\bf u}^L_{j+1/2}, \hat{\bf u}^R_{j+1/2} \right), \end{equation}

where \(\mathcal{U}\) denotes the upwinding function UpwindFunction() passed as an argument (if NULL, DefaultUpwinding() is used). The upwinding function is specified by the physical model.

Note: Solution-dependent, nonlinear interpolation methods (such as WENO, CRWENO) are implemented in a way that separates the calculation of the nonlinear interpolation weights (based on, say, the smoothness of the flux function), and the actual evaluation of the interpolant, into different functions. This allows the flexibility to choose if and when the nonlinear coefficients are evaluated (or previously computed values are reused). Some possible scenarios are:

  • For explicit time integration, they are computed every time the hyperbolic flux term is being computed.
  • For implicit time integration, consistency or linearization may require that they be computed and "frozen" at the beginning of a time step or stage.

The argument LimFlag controls this behavior:

  • LimFlag = 1 means recompute the nonlinear coefficients.
  • LimFlag = 0 means reuse the the previously computed coefficients.
Parameters
fluxIArray to hold the computed interface fluxes. This array does not have ghost points. The dimensions are the same as those of u without ghost points in all dimensions, except along dir, where it is one more
fluxCArray of the flux function computed at the cell centers (same layout as u)
uSolution array
xArray of spatial coordinates
dirSpatial dimension along which to reconstruct the interface fluxes
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent solution time
LimFlagFlag to indicate if the nonlinear coefficients for solution-dependent interpolation method should be recomputed
UpwindFunctionFunction pointer to the upwinding function for the interface flux computation. If NULL, DefaultUpwinding() will be used.

Definition at line 167 of file HyperbolicFunction.c.

185 {
186  HyPar *solver = (HyPar*) s;
187  MPIVariables *mpi = (MPIVariables*) m;
189 
190  double *uC = NULL;
191  double *uL = solver->uL;
192  double *uR = solver->uR;
193  double *fluxL = solver->fL;
194  double *fluxR = solver->fR;
195 
196  /*
197  precalculate the non-linear interpolation coefficients if required
198  else reuse the weights previously calculated
199  */
200 
201  if (LimFlag) IERR solver->SetInterpLimiterVar(fluxC,u,x,dir,solver,mpi);
202 
203  /* if defined, calculate the modified u-function to be used for upwinding
204  e.g.: used in well-balanced schemes for Euler/Navier-Stokes with gravity
205  otherwise, just copy u to uC */
206  if (solver->UFunction) {
207  uC = solver->uC;
208  IERR solver->UFunction(uC,u,dir,solver,mpi,t); CHECKERR(ierr);
209  } else uC = u;
210 
211  /* Interpolation -> to calculate left and right-biased interface flux and state variable*/
212  IERR solver->InterpolateInterfacesHyp(uL ,uC ,u,x, 1,dir,solver,mpi,1); CHECKERR(ierr);
213  IERR solver->InterpolateInterfacesHyp(uR ,uC ,u,x,-1,dir,solver,mpi,1); CHECKERR(ierr);
214  IERR solver->InterpolateInterfacesHyp(fluxL,fluxC,u,x, 1,dir,solver,mpi,0); CHECKERR(ierr);
215  IERR solver->InterpolateInterfacesHyp(fluxR,fluxC,u,x,-1,dir,solver,mpi,0); CHECKERR(ierr);
216 
217  /* Upwind -> to calculate the final interface flux */
218  if (UpwindFunction) { IERR UpwindFunction (fluxI,fluxL,fluxR,uL ,uR ,u ,dir,solver,t); CHECKERR(ierr); }
219  else { IERR DefaultUpwinding (fluxI,fluxL,fluxR,NULL,NULL,NULL,dir,solver,t); CHECKERR(ierr); }
220 
221  return(0);
222 }
double * uR
Definition: hypar.h:139
double * uL
Definition: hypar.h:139
static int DefaultUpwinding(double *, double *, double *, double *, double *, double *, int, void *, double)
double * fL
Definition: hypar.h:139
int(* SetInterpLimiterVar)(double *, double *, double *, int, void *, void *)
Definition: hypar.h:234
double * uC
Definition: hypar.h:131
int(* InterpolateInterfacesHyp)(double *, double *, double *, double *, int, int, void *, void *, int)
Definition: hypar.h:224
#define CHECKERR(ierr)
Definition: basic.h:18
int(* UFunction)(double *, double *, int, void *, void *, double)
Definition: hypar.h:321
#define IERR
Definition: basic.h:16
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * fR
Definition: hypar.h:139
#define _DECLARE_IERR_
Definition: basic.h:17
int DefaultUpwinding ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)
static

If no upwinding scheme is specified, this function defines the "upwind" flux as the arithmetic mean of the left- and right-biased fluxes.

Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension
sSolver object of type HyPar
tCurrent solution time

Definition at line 226 of file HyperbolicFunction.c.

237 {
238  HyPar *solver = (HyPar*) s;
239  int done;
240 
241  int *dim = solver->dim_local;
242  int ndims = solver->ndims;
243  int nvars = solver->nvars;
244 
245  int bounds_outer[ndims], bounds_inter[ndims];
246  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
247  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
248 
249  done = 0; int index_outer[ndims], index_inter[ndims];
250  _ArraySetValue_(index_outer,ndims,0);
251  while (!done) {
252  _ArrayCopy1D_(index_outer,index_inter,ndims);
253  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
254  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
255  int v; for (v=0; v<nvars; v++) fI[nvars*p+v] = 0.5 * (fL[nvars*p+v]+fR[nvars*p+v]);
256  }
257  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
258  }
259 
260  return(0);
261 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
int * dim_local
Definition: hypar.h:37
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int HyperbolicFunction ( double *  hyp,
double *  u,
void *  s,
void *  m,
double  t,
int  LimFlag,
int(*)(double *, double *, int, void *, double)  FluxFunction,
int(*)(double *, double *, double *, double *, double *, double *, int, void *, double)  UpwindFunction 
)

This function computes the hyperbolic term of the governing equations, expressed as follows:-

\begin{equation} {\bf F} \left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {\partial {\bf f}_d\left({\bf u}\right)} {\partial x_d} \end{equation}

using a conservative finite-difference discretization is space:

\begin{equation} \hat{\bf F} \left({\bf u}\right) = \sum_{d=0}^{D-1} \frac {1} {\Delta x_d} \left[ \hat{\bf f}_{d,j+1/2} - \hat{\bf f}_{d,j-1/2} \right] \end{equation}

where \(d\) denotes the spatial dimension, \(D\) denotes the total number of spatial dimensions, the hat denotes the discretized quantity, and \(j\) is the grid coordinate along dimension \(d\). The approximation to the flux function \({\bf f}_d\) at the interfaces \(j\pm1/2\), denoted by \(\hat{\bf f}_{d,j\pm 1/2}\), are computed using the function ReconstructHyperbolic().

Parameters
hypArray to hold the computed hyperbolic term (shares the same layout as u
uSolution array
sSolver object of type HyPar
mMPI object of type MPIVariables
tCurrent simulation time
LimFlagFlag to indicate if the nonlinear coefficients for solution-dependent interpolation method should be recomputed (see ReconstructHyperbolic() for an explanation on why this is needed)
FluxFunctionFunction pointer to the flux function for the hyperbolic term
UpwindFunctionFunction pointer to the upwinding function for the hyperbolic term

Definition at line 31 of file HyperbolicFunction.c.

45 {
46  HyPar *solver = (HyPar*) s;
47  MPIVariables *mpi = (MPIVariables*) m;
48  int d, v, i, done;
49  double *FluxI = solver->fluxI; /* interface flux */
50  double *FluxC = solver->fluxC; /* cell centered flux */
52 
53  int ndims = solver->ndims;
54  int nvars = solver->nvars;
55  int ghosts = solver->ghosts;
56  int *dim = solver->dim_local;
57  int size = solver->npoints_local_wghosts;
58  double *x = solver->x;
59  double *dxinv = solver->dxinv;
60  int index[ndims], index1[ndims], index2[ndims], dim_interface[ndims];
61 
62  LimFlag = (LimFlag && solver->flag_nonlinearinterp && solver->SetInterpLimiterVar);
63 
64  _ArraySetValue_(hyp,size*nvars,0.0);
65  _ArraySetValue_(solver->StageBoundaryIntegral,2*ndims*nvars,0.0);
66  if (!FluxFunction) return(0); /* zero hyperbolic term */
67  solver->count_hyp++;
68 
69 #if defined(CPU_STAT)
70  double cpu_time = 0.0;
71  clock_t cpu_start, cpu_end;
72 #endif
73 
74  int offset = 0;
75  for (d = 0; d < ndims; d++) {
76  _ArrayCopy1D_(dim,dim_interface,ndims); dim_interface[d]++;
77  int size_cellcenter = 1; for (i = 0; i < ndims; i++) size_cellcenter *= (dim[i] + 2*ghosts);
78  int size_interface = 1; for (i = 0; i < ndims; i++) size_interface *= dim_interface[i];
79 
80  /* evaluate cell-centered flux */
81  IERR FluxFunction(FluxC,u,d,solver,t); CHECKERR(ierr);
82  /* compute interface fluxes */
83  IERR ReconstructHyperbolic(FluxI,FluxC,u,x+offset,d,solver,mpi,t,LimFlag,UpwindFunction);
84  CHECKERR(ierr);
85 
86  /* calculate the first derivative */
87  done = 0; _ArraySetValue_(index,ndims,0);
88  int p, p1, p2;
89 
90 #if defined(CPU_STAT)
91  cpu_start = clock();
92 #endif
93 
94  while (!done) {
95  _ArrayCopy1D_(index,index1,ndims);
96  _ArrayCopy1D_(index,index2,ndims); index2[d]++;
97  _ArrayIndex1D_(ndims,dim ,index ,ghosts,p);
98  _ArrayIndex1D_(ndims,dim_interface,index1,0 ,p1);
99  _ArrayIndex1D_(ndims,dim_interface,index2,0 ,p2);
100  for (v=0; v<nvars; v++) hyp[nvars*p+v] += dxinv[offset+ghosts+index[d]]
101  * (FluxI[nvars*p2+v]-FluxI[nvars*p1+v]);
102  /* boundary flux integral */
103  if (index[d] == 0)
104  for (v=0; v<nvars; v++) solver->StageBoundaryIntegral[(2*d+0)*nvars+v] -= FluxI[nvars*p1+v];
105  if (index[d] == dim[d]-1)
106  for (v=0; v<nvars; v++) solver->StageBoundaryIntegral[(2*d+1)*nvars+v] += FluxI[nvars*p2+v];
107 
108  _ArrayIncrementIndex_(ndims,dim,index,done);
109  }
110 
111 #if defined(CPU_STAT)
112  cpu_end = clock();
113  cpu_time += (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC;
114 #endif
115 
116  offset += dim[d] + 2*ghosts;
117  }
118 
119 #if defined(CPU_STAT)
120  printf("HyperbolicFunction CPU time = %8.6lf\n", cpu_time);
121 #endif
122 
123  if (solver->flag_ib) _ArrayBlockMultiply_(hyp,solver->iblank,size,nvars);
124 
125  return(0);
126 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
double * iblank
Definition: hypar.h:436
#define _ArrayIncrementIndex_(N, imax, i, done)
int flag_ib
Definition: hypar.h:441
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
int(* SetInterpLimiterVar)(double *, double *, double *, int, void *, void *)
Definition: hypar.h:234
static int ReconstructHyperbolic(double *, double *, double *, double *, int, void *, void *, double, int, int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
double * fluxI
Definition: hypar.h:136
#define _ArrayBlockMultiply_(x, a, n, bs)
#define _ArrayCopy1D_(x, y, size)
double * fluxC
Definition: hypar.h:128
int nvars
Definition: hypar.h:29
#define CHECKERR(ierr)
Definition: basic.h:18
double * StageBoundaryIntegral
Definition: hypar.h:382
int count_hyp
Definition: hypar.h:418
double * dxinv
Definition: hypar.h:110
int ndims
Definition: hypar.h:26
#define IERR
Definition: basic.h:16
Structure of MPI-related variables.
double * x
Definition: hypar.h:107
int flag_nonlinearinterp
Definition: hypar.h:411
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _DECLARE_IERR_
Definition: basic.h:17