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

CRWENO5 Scheme (Component-wise application to vectors). More...

#include <stdio.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <interpolation.h>
#include <tridiagLU.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Macros

#define _MINIMUM_GHOSTS_   3
 

Functions

int Interp1PrimFifthOrderCRWENO (double *fI, double *fC, double *u, double *x, int upw, int dir, void *s, void *m, int uflag)
 5th order CRWENO reconstruction (component-wise) on a uniform grid More...
 

Detailed Description

CRWENO5 Scheme (Component-wise application to vectors).

Author
Debojyoti Ghosh

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

Function Documentation

◆ Interp1PrimFifthOrderCRWENO()

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

5th order CRWENO reconstruction (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 fifth order CRWENO 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 CRWENO 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[ \frac{2}{3}\hat{\bf f}_{j-1/2} + \frac{1}{3}\hat{\bf f}_{j+1/2} = \frac{1}{6} \left( f_{j-1} + 5f_j \right) \right]\\ + &\ \omega_2\ \times\ \left[ \frac{1}{3}\hat{\bf f}_{j-1/2}+\frac{2}{3}\hat{\bf f}_{j+1/2} = \frac{1}{6} \left( 5f_j + f_{j+1} \right) \right] \\ + &\ \omega_3\ \times\ \left[ \frac{2}{3}\hat{\bf f}_{j+1/2} + \frac{1}{3}\hat{\bf f}_{j+3/2} = \frac{1}{6} \left( f_j + 5f_{j+1} \right) \right] \\ \Rightarrow &\ \left(\frac{2}{3}\omega_1+\frac{1}{3}\omega_2\right)\hat{\bf f}_{j-1/2} + \left[\frac{1}{3}\omega_1+\frac{2}{3}(\omega_2+\omega_3)\right]\hat{\bf f}_{j+1/2} + \frac{1}{3}\omega_3\hat{\bf f}_{j+3/2} = \frac{\omega_1}{6}{\bf f}_{j-1} + \frac{5(\omega_1+\omega_2)+\omega_3}{6}{\bf f}_j + \frac{\omega_2+5\omega_3}{6}{\bf f}_{j+1}, \end{align}

where \(\omega_k; k=1,2,3\) are the nonlinear WENO weights computed in WENOFifthOrderCalculateWeights() (note that the \(\omega\) are different for each component of the vector \(\hat{\bf f}\)). The resulting tridiagonal system is solved using tridiagLU() (see also TridiagLU, tridiagLU.h).

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 scalar interpolation method is applied to the vector function in a component-wise manner.
  • 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:

  • Ghosh, D., Baeder, J. D., Compact Reconstruction Schemes with Weighted ENO Limiting for Hyperbolic Conservation Laws, SIAM Journal on Scientific Computing, 34 (3), 2012, A1678–A1706, http://dx.doi.org/10.1137/110857659
  • Ghosh, D., Constantinescu, E. M., Brown, J., Efficient Implementation of Nonlinear Compact Schemes on Massively Parallel Platforms, SIAM Journal on Scientific Computing, 37 (3), 2015, C354–C383, http://dx.doi.org/10.1137/140989261
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 78 of file Interp1PrimFifthOrderCRWENO.c.

89 {
90  HyPar *solver = (HyPar*) s;
91  MPIVariables *mpi = (MPIVariables*) m;
92  CompactScheme *compact= (CompactScheme*) solver->compact;
93  WENOParameters *weno = (WENOParameters*) solver->interp;
94  TridiagLU *lu = (TridiagLU*) solver->lusolver;
95  int sys,Nsys,d;
97 
98  int ghosts = solver->ghosts;
99  int ndims = solver->ndims;
100  int nvars = solver->nvars;
101  int *dim = solver->dim_local;
102  int *stride= solver->stride_with_ghosts;
103 
104  /* define some constants */
105  static const double one_third = 1.0/3.0;
106  static const double one_sixth = 1.0/6.0;
107 
108  double *ww1, *ww2, *ww3;
109  ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
110  ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
111  ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
112 
113  /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
114  dimension "dir" */
115  int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
116  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
117  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
118  int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
119 
120  /* calculate total number of tridiagonal systems to solve */
121  _ArrayProduct1D_(bounds_outer,ndims,Nsys); Nsys *= nvars;
122 
123  /* Allocate arrays for tridiagonal system */
124  double *A = compact->A;
125  double *B = compact->B;
126  double *C = compact->C;
127  double *R = compact->R;
128 
129 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
130  for (sys=0; sys < N_outer; sys++) {
131  _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
132  _ArrayCopy1D_(index_outer,indexC,ndims);
133  _ArrayCopy1D_(index_outer,indexI,ndims);
134  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
135  int qm1,qm2,qm3,qp1,qp2,p;
136  _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
137  if (upw > 0) {
138  indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
139  qm3 = qm1 - 2*stride[dir];
140  qm2 = qm1 - stride[dir];
141  qp1 = qm1 + stride[dir];
142  qp2 = qm1 + 2*stride[dir];
143  } else {
144  indexC[dir] = indexI[dir] ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
145  qm3 = qm1 + 2*stride[dir];
146  qm2 = qm1 + stride[dir];
147  qp1 = qm1 - stride[dir];
148  qp2 = qm1 - 2*stride[dir];
149  }
150 
151  /* Defining stencil points */
152  double *fm3, *fm2, *fm1, *fp1, *fp2;
153  fm3 = fC+qm3*nvars;
154  fm2 = fC+qm2*nvars;
155  fm1 = fC+qm1*nvars;
156  fp1 = fC+qp1*nvars;
157  fp2 = fC+qp2*nvars;
158 
159  /* Candidate stencils and their optimal weights*/
160  double f1[nvars], f2[nvars], f3[nvars];
161  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
162  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
163  /* Use WENO5 at the physical boundaries */
164  _ArrayAXBYCZ_(f1,(2*one_sixth),fm3,(-7*one_sixth) ,fm2,(11*one_sixth) ,fm1,nvars);
165  _ArrayAXBYCZ_(f2,(-one_sixth) ,fm2,(5*one_sixth) ,fm1,(2*one_sixth) ,fp1,nvars);
166  _ArrayAXBYCZ_(f3,(2*one_sixth),fm1,(5*one_sixth) ,fp1,(-one_sixth) ,fp2,nvars);
167  } else {
168  /* CRWENO5 at the interior points */
169  _ArrayAXBY_(f1,(one_sixth) ,fm2,(5*one_sixth),fm1,nvars);
170  _ArrayAXBY_(f2,(5*one_sixth),fm1,(one_sixth) ,fp1,nvars);
171  _ArrayAXBY_(f3,(one_sixth) ,fm1,(5*one_sixth),fp1,nvars);
172  }
173 
174  /* retrieve the WENO weights */
175  double *w1, *w2, *w3;
176  w1 = (ww1+p*nvars);
177  w2 = (ww2+p*nvars);
178  w3 = (ww3+p*nvars);
179 
180  if ( ((mpi->ip[dir] == 0 ) && (indexI[dir] == 0 ))
181  || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
182  _ArraySetValue_ ((A+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
183  _ArraySetValue_ ((B+Nsys*indexI[dir]+sys*nvars),nvars,1.0)
184  _ArraySetValue_ ((C+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
185  } else {
186  if (upw > 0) {
187  _ArrayAXBY_ ((A+Nsys*indexI[dir]+sys*nvars),(2*one_third) ,w1,(one_third) ,w2,nvars);
188  _ArrayAXBYCZ_ ((B+Nsys*indexI[dir]+sys*nvars),(one_third) ,w1,(2*one_third),w2,(2*one_third),w3,nvars);
189  _ArrayScaleCopy1D_(w3,(one_third),(C+Nsys*indexI[dir]+sys*nvars),nvars);
190  } else {
191  _ArrayAXBY_ ((C+Nsys*indexI[dir]+sys*nvars),(2*one_third) ,w1,(one_third) ,w2,nvars);
192  _ArrayAXBYCZ_ ((B+Nsys*indexI[dir]+sys*nvars),(one_third) ,w1,(2*one_third),w2,(2*one_third),w3,nvars);
193  _ArrayScaleCopy1D_(w3,(one_third),(A+Nsys*indexI[dir]+sys*nvars),nvars);
194  }
195  }
196  _ArrayMultiply3Add1D_ ((R+Nsys*indexI[dir]+sys*nvars),w1,f1,w2,f2,w3,f3,nvars);
197  }
198  }
199 
200 #ifdef serial
201 
202  /* Solve the tridiagonal system */
203  IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,NULL); CHECKERR(ierr);
204 
205 #else
206 
207  /* Solve the tridiagonal system */
208  /* all processes except the last will solve without the last interface to avoid overlap */
209  if (mpi->ip[dir] != mpi->iproc[dir]-1) { IERR tridiagLU(A,B,C,R,dim[dir] ,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
210  else { IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
211 
212  /* Now get the solution to the last interface from the next proc */
213  double *sendbuf = compact->sendbuf;
214  double *recvbuf = compact->recvbuf;
215  MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
216  if (mpi->ip[dir]) for (d=0; d<Nsys; d++) sendbuf[d] = R[d];
217  if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
218  if (mpi->ip[dir]) MPI_Isend(sendbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
219  MPI_Status status_arr[2];
220  MPI_Waitall(2,&req[0],status_arr);
221  if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys; d++) R[d+Nsys*dim[dir]] = recvbuf[d];
222 
223 #endif
224 
225  /* save the solution to fI */
226 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
227  for (sys=0; sys < N_outer; sys++) {
228  _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
229  _ArrayCopy1D_(index_outer,indexI,ndims);
230  for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
231  int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
232  _ArrayCopy1D_((R+sys*nvars+Nsys*indexI[dir]),(fI+nvars*p),nvars);
233  }
234  }
235 
236  return(0);
237 }
int nvars
Definition: hypar.h:29
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
#define _ArrayAXBY_(z, a, x, b, y, size)
#define _ArrayIndexnD_(N, index, imax, i, ghost)
int * stride_with_ghosts
Definition: hypar.h:414
int ndims
Definition: hypar.h:26
#define _ArrayScaleCopy1D_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _ArrayAXBYCZ_(w, a, x, b, y, c, z, size)
MPI_Comm * comm
#define _ArrayIndex1D_(N, imax, i, ghost, index)
void * lusolver
Definition: hypar.h:368
double * recvbuf
Structure of variables/parameters needed by the WENO-type scheme.
Structure of variables/parameters needed by the compact schemes.
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
void * interp
Definition: hypar.h:362
int tridiagLU(double *, double *, double *, double *, int, int, void *, void *)
Definition: tridiagLU.c:83
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
#define _DECLARE_IERR_
Definition: basic.h:17
#define _ArrayCopy1D_(x, y, size)
#define _ArrayProduct1D_(x, size, p)
void * compact
Definition: hypar.h:364
double * sendbuf
#define _ArrayMultiply3Add1D_(x, a, b, c, d, e, f, size)