Definitions for the functions computing the interpolated value of the primitive at the cell interfaces from the cell-centered values. More...

#include "basic.h"

Data Structures
struct	MUSCLParameters
	Structure of variables/parameters needed by the MUSCL scheme. More...

struct	WENOParameters
	Structure of variables/parameters needed by the WENO-type scheme. More...

struct	CompactScheme
	Structure of variables/parameters needed by the compact schemes. More...

Macros
#define	_FIRST_ORDER_UPWIND_ "1"

#define	_SECOND_ORDER_CENTRAL_ "2"

#define	_SECOND_ORDER_MUSCL_ "muscl2"

#define	_THIRD_ORDER_MUSCL_ "muscl3"

#define	_FOURTH_ORDER_CENTRAL_ "4"

#define	_FIFTH_ORDER_UPWIND_ "upw5"

#define	_FIFTH_ORDER_COMPACT_UPWIND_ "cupw5"

#define	_FIFTH_ORDER_WENO_ "weno5"

#define	_FIFTH_ORDER_CRWENO_ "crweno5"

#define	_FIFTH_ORDER_HCWENO_ "hcweno5"

#define	_CHARACTERISTIC_ "characteristic"

#define	_COMPONENTS_ "components"

#define	_WENO_OPTIMAL_WEIGHT_1_ 0.1

#define	_WENO_OPTIMAL_WEIGHT_2_ 0.6

#define	_WENO_OPTIMAL_WEIGHT_3_ 0.3

#define	_CRWENO_OPTIMAL_WEIGHT_1_ 0.2

#define	_CRWENO_OPTIMAL_WEIGHT_2_ 0.5

#define	_CRWENO_OPTIMAL_WEIGHT_3_ 0.3

#define	_WENOWeights_v_JS_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)

#define	_WENOWeights_v_M_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)

#define	_WENOWeights_v_M_Scalar_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, idx)

#define	_WENOWeights_v_Z_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)

#define	_WENOWeights_v_YC_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, N)

#define	_WENOWeights_v_YC_Scalar_(w1, w2, w3, c1, c2, c3, m3, m2, m1, p1, p2, eps, idx)

Functions
int	Interp1PrimFirstOrderUpwind (double , double , double , double , int, int, void , void , int)
	1st order upwind reconstruction (component-wise) on a uniform grid More...

int	Interp1PrimSecondOrderCentral (double , double , double , double , int, int, void , void , int)
	2nd order central reconstruction (component-wise) on a uniform grid More...

int	Interp1PrimSecondOrderMUSCL (double , double , double , double , int, int, void , void , int)
	2nd order MUSCL scheme (component-wise) on a uniform grid More...

int	Interp1PrimThirdOrderMUSCL (double , double , double , double , int, int, void , void , int)
	3rd order MUSCL scheme with Koren's limiter (component-wise) on a uniform grid More...

int	Interp1PrimFourthOrderCentral (double , double , double , double , int, int, void , void , int)
	4th order central reconstruction (component-wise) on a uniform grid More...

int	Interp1PrimFifthOrderUpwind (double , double , double , double , int, int, void , void , int)
	5th order upwind reconstruction (component-wise) on a uniform grid More...

int	Interp1PrimFifthOrderCompactUpwind (double , double , double , double , int, int, void , void , int)
	5th order compact upwind reconstruction (component-wise) on a uniform grid More...

int	Interp1PrimFifthOrderWENO (double , double , double , double , int, int, void , void , int)
	5th order WENO reconstruction (component-wise) on a uniform grid More...

int	Interp1PrimFifthOrderCRWENO (double , double , double , double , int, int, void , void , int)
	5th order CRWENO reconstruction (component-wise) on a uniform grid More...

int	Interp1PrimFifthOrderHCWENO (double , double , double , double , int, int, void , void , int)
	5th order hybrid compact-WENO reconstruction (component-wise) on a uniform grid More...

int	gpuInterp1PrimFifthOrderWENO (double , double , double , double , int, int, void , void , int)

int	Interp1PrimFirstOrderUpwindChar (double , double , double , double , int, int, void , void , int)
	1st order upwind reconstruction (characteristic-based) on a uniform grid More...

int	Interp1PrimSecondOrderCentralChar (double , double , double , double , int, int, void , void , int)
	2nd order central reconstruction (characteristic-based) on a uniform grid More...

int	Interp1PrimSecondOrderMUSCLChar (double , double , double , double , int, int, void , void , int)
	2nd order MUSCL scheme (characteristic-based) on a uniform grid More...

int	Interp1PrimThirdOrderMUSCLChar (double , double , double , double , int, int, void , void , int)
	3rd order MUSCL scheme with Koren's limiter (characteristic-based) on a uniform grid More...

int	Interp1PrimFourthOrderCentralChar (double , double , double , double , int, int, void , void , int)
	4th order central reconstruction (characteristic-based) on a uniform grid More...

int	Interp1PrimFifthOrderUpwindChar (double , double , double , double , int, int, void , void , int)
	5th order upwind reconstruction (characteristic-based) on a uniform grid More...

int	Interp1PrimFifthOrderCompactUpwindChar (double , double , double , double , int, int, void , void , int)
	5th order compact upwind reconstruction (characteristic-based) on a uniform grid More...

int	Interp1PrimFifthOrderWENOChar (double , double , double , double , int, int, void , void , int)
	5th order WENO reconstruction (characteristic-based) on a uniform grid More...

int	Interp1PrimFifthOrderCRWENOChar (double , double , double , double , int, int, void , void , int)
	5th order CRWENO reconstruction (characteristic-based) on a uniform grid More...

int	Interp1PrimFifthOrderHCWENOChar (double , double , double , double , int, int, void , void , int)
	5th order hybrid compact-WENO reconstruction (characteristic-based) on a uniform grid More...

int	Interp2PrimSecondOrder (double , double , int, void , void )
	2nd order component-wise interpolation of the 2nd primitive on a uniform grid More...

int	InterpSetLimiterVar (double , double , double , int, void , void *)

int	MUSCLInitialize (void , void )

int	WENOInitialize (void , void , char , char )

int	WENOCleanup (void *, int)

int	CompactSchemeInitialize (void , void , char *)

int	CompactSchemeCleanup (void *)

Detailed Description

Definitions for the functions computing the interpolated value of the primitive at the cell interfaces from the cell-centered values.

Author: Debojyoti Ghosh

Definition in file interpolation.h.

Data Structure Documentation

◆ WENOParameters

struct WENOParameters

Structure of variables/parameters needed by the WENO-type scheme.

This structure contains the variables/parameters needed by the WENO-type scheme (_FIFTH_ORDER_WENO_, _FIFTH_ORDER_CRWENO_, _FIFTH_ORDER_HCWENO_).

Definition at line 194 of file interpolation.h.

Data Fields
int	mapped	Use mapped weights? (Henrick, Aslam, J. Comput. Phys., 2005)
int	borges	Use Borges' implementation of weights? (Borges, et. al, J. Comput. Phys., 2008)
int	yc	Use Yamaleev-Carpenter implementation of weights? (Yamaleev, Carpenter, J. Comput. Phys., 2009)
int	no_limiting	Remove limiting -> 5th order polynomial interpolation (freeze the WENO weights to the optimal coefficients)
double	eps	epsilon parameter
double	p	p parameter
double	tol	a general tolerance parameter
double	rc	Parameter for the hybrid compact-WENO scheme
double	xi	Parameter for the hybrid compact-WENO scheme
double *	w1	Array to save the first WENO weight
double *	w2	Array to save the second WENO weight
double *	w3	Array to save the third WENO weight
int *	offset	Array containing the offset information for the WENO weights
int	size	Size of the WENO weights array

◆ CompactScheme

struct CompactScheme

Structure of variables/parameters needed by the compact schemes.

This structure contains the variables/parameters needed by a compact scheme (_FIFTH_ORDER_COMPACT_UPWIND_, _FIFTH_ORDER_CRWENO_, _FIFTH_ORDER_HCWENO_).

Definition at line 560 of file interpolation.h.

Data Fields
double *	A	Array to save the sub-diagonal of the tridiagonal system resulting from the fifth-order CRWENO scheme
double *	B	Array to save the diagonal of the tridiagonal system resulting from the fifth-order CRWENO scheme
double *	C	Array to save the super-diagonal of the tridiagonal system resulting from the fifth-order CRWENO scheme
double *	R	Array to save the right-hand-side of the tridiagonal system resulting from the fifth-order CRWENO scheme
double *	sendbuf	Buffer array to send data across processors
double *	recvbuf	Buffer array to receive data across processors

Macro Definition Documentation

◆ _FIRST_ORDER_UPWIND_

#define _FIRST_ORDER_UPWIND_ "1"

First order upwind scheme: Interp1PrimFirstOrderUpwind(), Interp1PrimFirstOrderUpwindChar()

Definition at line 12 of file interpolation.h.

◆ _SECOND_ORDER_CENTRAL_

#define _SECOND_ORDER_CENTRAL_ "2"

Second order central scheme: Interp1PrimSecondOrderCentral(), Interp1PrimSecondOrderCentralChar()

Definition at line 14 of file interpolation.h.

◆ _SECOND_ORDER_MUSCL_

#define _SECOND_ORDER_MUSCL_ "muscl2"

Second order MUSCL scheme: Interp1PrimSecondOrderMUSCL(), Interp1PrimSecondOrderMUSCLChar()

Definition at line 16 of file interpolation.h.

◆ _THIRD_ORDER_MUSCL_

#define _THIRD_ORDER_MUSCL_ "muscl3"

Third order MUSCL scheme with Koren's limiter: Interp1PrimThirdOrderMUSCL(), Interp1PrimThirdOrderMUSCLChar()

Definition at line 18 of file interpolation.h.

◆ _FOURTH_ORDER_CENTRAL_

#define _FOURTH_ORDER_CENTRAL_ "4"

Fourth order central scheme: Interp1PrimFourthOrderCentral(), Interp1PrimFourthOrderCentralChar()

Definition at line 20 of file interpolation.h.

◆ _FIFTH_ORDER_UPWIND_

#define _FIFTH_ORDER_UPWIND_ "upw5"

Fifth order upwind scheme: Interp1PrimFifthOrderUpwind(), Interp1PrimFifthOrderUpwindChar()

Definition at line 22 of file interpolation.h.

◆ _FIFTH_ORDER_COMPACT_UPWIND_

#define _FIFTH_ORDER_COMPACT_UPWIND_ "cupw5"

Fifth order compact upwind scheme: Interp1PrimFifthOrderCompactUpwind(), Interp1PrimFifthOrderCompactUpwindChar()

Definition at line 24 of file interpolation.h.

◆ _FIFTH_ORDER_WENO_

#define _FIFTH_ORDER_WENO_ "weno5"

Fifth order Weighted Essentially Non-Oscillatory (WENO) scheme: Interp1PrimFifthOrderWENO(), Interp1PrimFifthOrderWENOChar()

Definition at line 26 of file interpolation.h.

◆ _FIFTH_ORDER_CRWENO_

#define _FIFTH_ORDER_CRWENO_ "crweno5"

Fifth order Compact Reconstruction Weighted Essentially Non-Oscillatory (CRWENO) scheme: Interp1PrimFifthOrderCRWENO(), Interp1PrimFifthOrderCRWENOChar()

Definition at line 28 of file interpolation.h.

◆ _FIFTH_ORDER_HCWENO_

#define _FIFTH_ORDER_HCWENO_ "hcweno5"

Fifth order hybrid compact-WENO scheme: Interp1PrimFifthOrderHCWENO(), Interp1PrimFifthOrderHCWENOChar()

Definition at line 30 of file interpolation.h.

◆ _CHARACTERISTIC_

#define _CHARACTERISTIC_ "characteristic"

Characteristic-based interpolation of vectors (Physical model must define left and right eigenvectors)

Definition at line 33 of file interpolation.h.

◆ _COMPONENTS_

#define _COMPONENTS_ "components"

Component-wise interpolation of vectors

Definition at line 34 of file interpolation.h.

◆ _WENO_OPTIMAL_WEIGHT_1_

#define _WENO_OPTIMAL_WEIGHT_1_ 0.1

Optimal value for the first fifth-order WENO weight

Definition at line 228 of file interpolation.h.

◆ _WENO_OPTIMAL_WEIGHT_2_

#define _WENO_OPTIMAL_WEIGHT_2_ 0.6

Optimal value for the second fifth-order WENO weight

Definition at line 230 of file interpolation.h.

◆ _WENO_OPTIMAL_WEIGHT_3_

#define _WENO_OPTIMAL_WEIGHT_3_ 0.3

Optimal value for the third fifth-order WENO weight

Definition at line 232 of file interpolation.h.

◆ _CRWENO_OPTIMAL_WEIGHT_1_

#define _CRWENO_OPTIMAL_WEIGHT_1_ 0.2

Optimal value for the first fifth-order CRWENO weight

Definition at line 234 of file interpolation.h.

◆ _CRWENO_OPTIMAL_WEIGHT_2_

#define _CRWENO_OPTIMAL_WEIGHT_2_ 0.5

Optimal value for the second fifth-order CRWENO weight

Definition at line 236 of file interpolation.h.

◆ _CRWENO_OPTIMAL_WEIGHT_3_

#define _CRWENO_OPTIMAL_WEIGHT_3_ 0.3

Optimal value for the third fifth-order CRWENO weight

Definition at line 238 of file interpolation.h.

◆ _WENOWeights_v_JS_

#define _WENOWeights_v_JS_	(	w1,
		w2,
		w3,
		c1,
		c2,
		c3,
		m3,
		m2,
		m1,
		p1,
		p2,
		eps,
		N
	)

Value:

{ \
    int idx; \
    /* calculate smoothness indicators and the WENO weights */\
    double b1, b2, b3, a1, a2, a3, a_sum_inv; \
    for (idx=0; idx<N; idx++) { \
      b1 = thirteen_by_twelve*(m3[idx]-2*m2[idx]+m1[idx])*(m3[idx]-2*m2[idx]+m1[idx]) \
           + one_fourth*(m3[idx]-4*m2[idx]+3*m1[idx])*(m3[idx]-4*m2[idx]+3*m1[idx]);  \
      a1 = c1 / ( (b1+eps) * (b1+eps) );  \
      b2 = thirteen_by_twelve*(m2[idx]-2*m1[idx]+p1[idx])*(m2[idx]-2*m1[idx]+p1[idx]) \
           + one_fourth*(m2[idx]-p1[idx])*(m2[idx]-p1[idx]);                \
      a2 = c2 / ( (b2+eps) * (b2+eps) );  \
      b3 = thirteen_by_twelve*(m1[idx]-2*p1[idx]+p2[idx])*(m1[idx]-2*p1[idx]+p2[idx]) \
           + one_fourth*(3*m1[idx]-4*p1[idx]+p2[idx])*(3*m1[idx]-4*p1[idx]+p2[idx]);  \
      a3 = c3 / ( (b3+eps) * (b3+eps) );  \
      a_sum_inv = 1.0 / (a1 + a2 + a3); \
      w1[idx] = a1 * a_sum_inv; \
      w2[idx] = a2 * a_sum_inv; \
      w3[idx] = a3 * a_sum_inv; \
    } \
  }

Compute the WENO weights according the the Jiang & Shu formulation:

\begin{equation} \omega_k = \frac {a_k} {\sum_{j=1}^3 a_j },\ a_k = \frac {c_k} {\left(\beta_k+\epsilon\right)^p},\ k = 1,2,3, \end{equation}

where \(c_k\) are the optimal weights, \(p\) is hardcoded to \(2\), and \(\epsilon\) is an input parameter (WENOParameters::eps) (typically \(10^{-6}\)). The smoothness indicators \(\beta_k\) are given by:

\begin{eqnarray} \beta_1 &=& \frac{13}{12} \left(f_{j-2}-2f_{j-1}+f_j\right)^2 + \frac{1}{4}\left(f_{j-2}-4f_{j-1}+3f_j\right)^2 \\ \beta_2 &=& \frac{13}{12} \left(f_{j-1}-2f_j+f_{j+1}\right)^2 + \frac{1}{4}\left(f_{j-1}-f_{j+1}\right)^2 \\ \beta_3 &=& \frac{13}{12} \left(f_j-2f_{j+1}+f_{j+2}\right)^2 + \frac{1}{4}\left(3f_j-4f_{j+1}+f_{j+2}\right)^2 \end{eqnarray}

Notes:

This macro computes the weights for one variable along one grid line.

Arguments:

w1, w2,w3 are the nonlinear WENO weights.
c1, c2,c3 are optimal coefficients.
m3, m2,m1,p1,p2 are the function values at stencil points corresponding to the interface j+1/2: j-2,j-1,j,j+1,j+2
weno is an object of type WENOParameters containing parameters for the WENO method.

Reference:

Jiang, Shu, J. Comput. Phys., 1996. http://dx.doi.org/10.1006/jcph.1996.0130

Definition at line 265 of file interpolation.h.

◆ _WENOWeights_v_M_

#define _WENOWeights_v_M_	(	w1,
		w2,
		w3,
		c1,
		c2,
		c3,
		m3,
		m2,
		m1,
		p1,
		p2,
		eps,
		N
	)

Value:

{ \
    int idx; \
    /* calculate smoothness indicators and the WENO weights */\
    double b1, b2, b3, a1, a2, a3, a_sum_inv; \
    for (idx=0; idx<N; idx++) { \
      b1 = thirteen_by_twelve*(m3[idx]-2*m2[idx]+m1[idx])*(m3[idx]-2*m2[idx]+m1[idx]) \
           + one_fourth*(m3[idx]-4*m2[idx]+3*m1[idx])*(m3[idx]-4*m2[idx]+3*m1[idx]);  \
      a1 = c1 / ( (b1+eps) * (b1+eps) );  \
      b2 = thirteen_by_twelve*(m2[idx]-2*m1[idx]+p1[idx])*(m2[idx]-2*m1[idx]+p1[idx]) \
           + one_fourth*(m2[idx]-p1[idx])*(m2[idx]-p1[idx]);                \
      a2 = c2 / ( (b2+eps) * (b2+eps) );  \
      b3 = thirteen_by_twelve*(m1[idx]-2*p1[idx]+p2[idx])*(m1[idx]-2*p1[idx]+p2[idx]) \
           + one_fourth*(3*m1[idx]-4*p1[idx]+p2[idx])*(3*m1[idx]-4*p1[idx]+p2[idx]);  \
      a3 = c3 / ( (b3+eps) * (b3+eps) );  \
      a_sum_inv = 1.0 / (a1 + a2 + a3); \
      w1[idx] = a1 * a_sum_inv; \
      w2[idx] = a2 * a_sum_inv; \
      w3[idx] = a3 * a_sum_inv; \
      a1 = w1[idx] * (c1 + c1*c1 - 3*c1*w1[idx] + w1[idx]*w1[idx]) / (c1*c1 + w1[idx]*(1.0-2.0*c1)); \
      a2 = w2[idx] * (c2 + c2*c2 - 3*c2*w2[idx] + w2[idx]*w2[idx]) / (c2*c2 + w2[idx]*(1.0-2.0*c2)); \
      a3 = w3[idx] * (c3 + c3*c3 - 3*c3*w3[idx] + w3[idx]*w3[idx]) / (c3*c3 + w3[idx]*(1.0-2.0*c3)); \
      a_sum_inv = 1.0 / (a1 + a2 + a3); \
      w1[idx] = a1 * a_sum_inv; \
      w2[idx] = a2 * a_sum_inv; \
      w3[idx] = a3 * a_sum_inv; \
    } \
  }

Compute the WENO weights according the the Mapped-WENO formulation:

\begin{eqnarray} \omega_k &=& \frac {a_k} {\sum_{j=1}^3 a_j },\ a_k = \frac {\tilde{\omega}_k \left( c_k + c_k^2 - 3c_k\tilde{\omega}_k + \tilde{\omega}_k^2\right)} {c_k^2 + \tilde{\omega}_k\left(1-2c_k\right)}, \\ \tilde{\omega}_k &=& \frac {\tilde{a}_k} {\sum_{j=1}^3 \tilde{a}_j },\ \tilde{a}_k = \frac {c_k} {\left(\beta_k+\epsilon\right)^p},\ k = 1,2,3, \end{eqnarray}

where \(c_k\) are the optimal weights, \(p\) is hardcoded to \(2\), and \(\epsilon\) is an input parameter (WENOParameters::eps) (typically \(10^{-6}\)). The smoothness indicators \(\beta_k\) are given by:

\begin{eqnarray} \beta_1 &=& \frac{13}{12} \left(f_{j-2}-2f_{j-1}+f_j\right)^2 + \frac{1}{4}\left(f_{j-2}-4f_{j-1}+3f_j\right)^2 \\ \beta_2 &=& \frac{13}{12} \left(f_{j-1}-2f_j+f_{j+1}\right)^2 + \frac{1}{4}\left(f_{j-1}-f_{j+1}\right)^2 \\ \beta_3 &=& \frac{13}{12} \left(f_j-2f_{j+1}+f_{j+2}\right)^2 + \frac{1}{4}\left(3f_j-4f_{j+1}+f_{j+2}\right)^2 \end{eqnarray}

Notes:

This macro computes the weights for one variable along one grid line.

Arguments:-

w1, w2,w3 are the nonlinear WENO weights.
c1, c2,c3 are optimal coefficients.
m3, m2,m1,p1,p2 are the function values at stencil points corresponding to the interface j+1/2: j-2,j-1,j,j+1,j+2
weno is an object of type WENOParameters containing parameters for the WENO method.
N is the number of interfaces along the grid line on which this WENO-type reconstruction is happening.

Reference:

Henrick, Aslam, Powers, J. Comput. Phys., 2005. http://dx.doi.org/10.1016/j.jcp.2005.01.023

Definition at line 314 of file interpolation.h.

◆ _WENOWeights_v_M_Scalar_

#define _WENOWeights_v_M_Scalar_	(	w1,
		w2,
		w3,
		c1,
		c2,
		c3,
		m3,
		m2,
		m1,
		p1,
		p2,
		eps,
		idx
	)

Value:

{ \
    /* calculate smoothness indicators and the WENO weights */\
    double b1, b2, b3, a1, a2, a3, a_sum_inv; \
    b1 = thirteen_by_twelve*(m3[idx]-2*m2[idx]+m1[idx])*(m3[idx]-2*m2[idx]+m1[idx]) \
          + one_fourth*(m3[idx]-4*m2[idx]+3*m1[idx])*(m3[idx]-4*m2[idx]+3*m1[idx]);  \
    a1 = c1 / ( (b1+eps) * (b1+eps) );  \
    b2 = thirteen_by_twelve*(m2[idx]-2*m1[idx]+p1[idx])*(m2[idx]-2*m1[idx]+p1[idx]) \
          + one_fourth*(m2[idx]-p1[idx])*(m2[idx]-p1[idx]);                \
    a2 = c2 / ( (b2+eps) * (b2+eps) );  \
    b3 = thirteen_by_twelve*(m1[idx]-2*p1[idx]+p2[idx])*(m1[idx]-2*p1[idx]+p2[idx]) \
          + one_fourth*(3*m1[idx]-4*p1[idx]+p2[idx])*(3*m1[idx]-4*p1[idx]+p2[idx]);  \
    a3 = c3 / ( (b3+eps) * (b3+eps) );  \
    a_sum_inv = 1.0 / (a1 + a2 + a3); \
    w1[idx] = a1 * a_sum_inv; \
    w2[idx] = a2 * a_sum_inv; \
    w3[idx] = a3 * a_sum_inv; \
    a1 = w1[idx] * (c1 + c1*c1 - 3*c1*w1[idx] + w1[idx]*w1[idx]) / (c1*c1 + w1[idx]*(1.0-2.0*c1)); \
    a2 = w2[idx] * (c2 + c2*c2 - 3*c2*w2[idx] + w2[idx]*w2[idx]) / (c2*c2 + w2[idx]*(1.0-2.0*c2)); \
    a3 = w3[idx] * (c3 + c3*c3 - 3*c3*w3[idx] + w3[idx]*w3[idx]) / (c3*c3 + w3[idx]*(1.0-2.0*c3)); \
    a_sum_inv = 1.0 / (a1 + a2 + a3); \
    w1[idx] = a1 * a_sum_inv; \
    w2[idx] = a2 * a_sum_inv; \
    w3[idx] = a3 * a_sum_inv; \
  }

Compute the WENO weights according the the Mapped-WENO formulation:

\begin{eqnarray} \omega_k &=& \frac {a_k} {\sum_{j=1}^3 a_j },\ a_k = \frac {\tilde{\omega}_k \left( c_k + c_k^2 - 3c_k\tilde{\omega}_k + \tilde{\omega}_k^2\right)} {c_k^2 + \tilde{\omega}_k\left(1-2c_k\right)}, \\ \tilde{\omega}_k &=& \frac {\tilde{a}_k} {\sum_{j=1}^3 \tilde{a}_j },\ \tilde{a}_k = \frac {c_k} {\left(\beta_k+\epsilon\right)^p},\ k = 1,2,3, \end{eqnarray}

where \(c_k\) are the optimal weights, \(p\) is hardcoded to \(2\), and \(\epsilon\) is an input parameter (WENOParameters::eps) (typically \(10^{-6}\)). The smoothness indicators \(\beta_k\) are given by:

\begin{eqnarray} \beta_1 &=& \frac{13}{12} \left(f_{j-2}-2f_{j-1}+f_j\right)^2 + \frac{1}{4}\left(f_{j-2}-4f_{j-1}+3f_j\right)^2 \\ \beta_2 &=& \frac{13}{12} \left(f_{j-1}-2f_j+f_{j+1}\right)^2 + \frac{1}{4}\left(f_{j-1}-f_{j+1}\right)^2 \\ \beta_3 &=& \frac{13}{12} \left(f_j-2f_{j+1}+f_{j+2}\right)^2 + \frac{1}{4}\left(3f_j-4f_{j+1}+f_{j+2}\right)^2 \end{eqnarray}

Notes:

This macro computes the weights for one variable along one grid line.

Arguments:-

w1, w2,w3 are the nonlinear WENO weights.
c1, c2,c3 are optimal coefficients.
m3, m2,m1,p1,p2 are the function values at stencil points corresponding to the interface j+1/2: j-2,j-1,j,j+1,j+2
weno is an object of type WENOParameters containing parameters for the WENO method.
N is the number of interfaces along the grid line on which this WENO-type reconstruction is happening.

Reference:

Henrick, Aslam, Powers, J. Comput. Phys., 2005. http://dx.doi.org/10.1016/j.jcp.2005.01.023

Definition at line 370 of file interpolation.h.

◆ _WENOWeights_v_Z_

#define _WENOWeights_v_Z_	(	w1,
		w2,
		w3,
		c1,
		c2,
		c3,
		m3,
		m2,
		m1,
		p1,
		p2,
		eps,
		N
	)

Value:

{ \
    int idx; \
    /* calculate smoothness indicators and the WENO weights */\
    double b1, b2, b3, a1, a2, a3, a_sum_inv, tau; \
    for (idx=0; idx<N; idx++) { \
      b1 = thirteen_by_twelve*(m3[idx]-2*m2[idx]+m1[idx])*(m3[idx]-2*m2[idx]+m1[idx]) \
           + one_fourth*(m3[idx]-4*m2[idx]+3*m1[idx])*(m3[idx]-4*m2[idx]+3*m1[idx]);  \
      b2 = thirteen_by_twelve*(m2[idx]-2*m1[idx]+p1[idx])*(m2[idx]-2*m1[idx]+p1[idx]) \
           + one_fourth*(m2[idx]-p1[idx])*(m2[idx]-p1[idx]);                \
      b3 = thirteen_by_twelve*(m1[idx]-2*p1[idx]+p2[idx])*(m1[idx]-2*p1[idx]+p2[idx]) \
           + one_fourth*(3*m1[idx]-4*p1[idx]+p2[idx])*(3*m1[idx]-4*p1[idx]+p2[idx]);  \
      tau = absolute(b3 - b1);  \
      a1 = c1 * (1.0 + (tau/(b1+eps)) * (tau/(b1+eps)) );  \
      a2 = c2 * (1.0 + (tau/(b2+eps)) * (tau/(b2+eps)) );  \
      a3 = c3 * (1.0 + (tau/(b3+eps)) * (tau/(b3+eps)) );  \
      a_sum_inv = 1.0 / (a1 + a2 + a3); \
      w1[idx] = a1 * a_sum_inv; \
      w2[idx] = a2 * a_sum_inv; \
      w3[idx] = a3 * a_sum_inv; \
    } \
  }

Compute the WENO weights according the the WENO-Z formulation:

\begin{equation} \omega_k = \frac {a_k} {\sum_{j=1}^3 a_j },\ a_k = c_k \left( 1 + \frac{\tau_5}{\beta_k+\epsilon} \right)^p,\ k = 1,2,3, \end{equation}

where \(c_k\) are the optimal weights, \(p\) is hardcoded to \(2\), and \(\epsilon\) is an input parameter (WENOParameters::eps) (typically \(10^{-6}\)). The smoothness indicators \(\beta_k\) are given by:

\begin{eqnarray} \beta_1 &=& \frac{13}{12} \left(f_{j-2}-2f_{j-1}+f_j\right)^2 + \frac{1}{4}\left(f_{j-2}-4f_{j-1}+3f_j\right)^2 \\ \beta_2 &=& \frac{13}{12} \left(f_{j-1}-2f_j+f_{j+1}\right)^2 + \frac{1}{4}\left(f_{j-1}-f_{j+1}\right)^2 \\ \beta_3 &=& \frac{13}{12} \left(f_j-2f_{j+1}+f_{j+2}\right)^2 + \frac{1}{4}\left(3f_j-4f_{j+1}+f_{j+2}\right)^2, \end{eqnarray}

and \(\tau_5 = \left|\beta_1 - \beta_3 \right|\).

Notes:

This macro computes the weights for one variable along one grid line.

Arguments:

w1, w2,w3 are the nonlinear WENO weights.
c1, c2,c3 are optimal coefficients.
m3, m2,m1,p1,p2 are the function values at stencil points corresponding to the interface j+1/2: j-2,j-1,j,j+1,j+2
weno is an object of type WENOParameters containing parameters for the WENO method.
N is the number of interfaces along the grid line on which this WENO-type reconstruction is happening.

Reference:

Borges, et. al., An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys., 2008. http://dx.doi.org/10.1016/j.jcp.2007.11.038
Castro, M., Costa, B., Don, W. S., High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys., 2011. http://dx.doi.org/10.1016/j.jcp.2010.11.028

Definition at line 426 of file interpolation.h.

◆ _WENOWeights_v_YC_

#define _WENOWeights_v_YC_	(	w1,
		w2,
		w3,
		c1,
		c2,
		c3,
		m3,
		m2,
		m1,
		p1,
		p2,
		eps,
		N
	)

Value:

{ \
    int idx; \
    /* calculate smoothness indicators and the WENO weights */\
    double b1, b2, b3, a1, a2, a3, a_sum_inv, tau; \
    for (idx=0; idx<N; idx++) { \
      b1 = thirteen_by_twelve*(m3[idx]-2*m2[idx]+m1[idx])*(m3[idx]-2*m2[idx]+m1[idx]) \
           + one_fourth*(m3[idx]-4*m2[idx]+3*m1[idx])*(m3[idx]-4*m2[idx]+3*m1[idx]);  \
      b2 = thirteen_by_twelve*(m2[idx]-2*m1[idx]+p1[idx])*(m2[idx]-2*m1[idx]+p1[idx]) \
           + one_fourth*(m2[idx]-p1[idx])*(m2[idx]-p1[idx]);                \
      b3 = thirteen_by_twelve*(m1[idx]-2*p1[idx]+p2[idx])*(m1[idx]-2*p1[idx]+p2[idx]) \
           + one_fourth*(3*m1[idx]-4*p1[idx]+p2[idx])*(3*m1[idx]-4*p1[idx]+p2[idx]);  \
      tau = (m3[idx]-4*m2[idx]+6*m1[idx]-4*p1[idx]+p2[idx])*(m3[idx]-4*m2[idx]+6*m1[idx]-4*p1[idx]+p2[idx]);  \
      a1 = c1 * (1.0 + (tau/(b1+eps)) * (tau/(b1+eps)) );  \
      a2 = c2 * (1.0 + (tau/(b2+eps)) * (tau/(b2+eps)) );  \
      a3 = c3 * (1.0 + (tau/(b3+eps)) * (tau/(b3+eps)) );  \
      a_sum_inv = 1.0 / (a1 + a2 + a3); \
      w1[idx] = a1 * a_sum_inv; \
      w2[idx] = a2 * a_sum_inv; \
      w3[idx] = a3 * a_sum_inv; \
    } \
  }

Compute the WENO weights according the the ESWENO formulation of Yamaleev & Carpenter. Note that only the formulation for the nonlinear weights is adopted and implemented here, not the ESWENO scheme as a whole.

\begin{equation} \omega_k = \frac {a_k} {\sum_{j=1}^3 a_j },\ a_k = c_k \left( 1 + \frac{\tau_5}{\beta_k+\epsilon} \right)^p,\ k = 1,2,3, \end{equation}

where \(c_k\) are the optimal weights, \(p\) is hardcoded to \(2\), and \(\epsilon\) is an input parameter (WENOParameters::eps) (typically \(10^{-6}\)). The smoothness indicators \(\beta_k\) are given by:

\begin{eqnarray} \beta_1 &=& \frac{13}{12} \left(f_{j-2}-2f_{j-1}+f_j\right)^2 + \frac{1}{4}\left(f_{j-2}-4f_{j-1}+3f_j\right)^2 \\ \beta_2 &=& \frac{13}{12} \left(f_{j-1}-2f_j+f_{j+1}\right)^2 + \frac{1}{4}\left(f_{j-1}-f_{j+1}\right)^2 \\ \beta_3 &=& \frac{13}{12} \left(f_j-2f_{j+1}+f_{j+2}\right)^2 + \frac{1}{4}\left(3f_j-4f_{j+1}+f_{j+2}\right)^2, \end{eqnarray}

and \(\tau_5 = \left( f_{j-2}-4f_{j-1}+6f_j-4f_{j+1}+f_{j+2} \right)^2\).

Notes:

This macro computes the weights for one variable along one grid line.

Arguments:

w1, w2,w3 are the nonlinear WENO weights.
c1, c2,c3 are optimal coefficients.
m3, m2,m1,p1,p2 are the function values at stencil points corresponding to the interface j+1/2: j-2,j-1,j,j+1,j+2
weno is an object of type WENOParameters containing parameters for the WENO method.
N is the number of interfaces along the grid line on which this WENO-type reconstruction is happening.
Reference:
- Yamaleev, Carpenter, A systematic methodology for constructing high-order energy stable WENO schemes, J. Comput. Phys., 2009. http://dx.doi.org/10.1016/j.jcp.2009.03.002

Definition at line 479 of file interpolation.h.

◆ _WENOWeights_v_YC_Scalar_

#define _WENOWeights_v_YC_Scalar_	(	w1,
		w2,
		w3,
		c1,
		c2,
		c3,
		m3,
		m2,
		m1,
		p1,
		p2,
		eps,
		idx
	)

Value:

{ \
    /* calculate smoothness indicators and the WENO weights */\
    double b1, b2, b3, a1, a2, a3, a_sum_inv, tau; \
    b1 = thirteen_by_twelve*(m3[idx]-2*m2[idx]+m1[idx])*(m3[idx]-2*m2[idx]+m1[idx]) \
          + one_fourth*(m3[idx]-4*m2[idx]+3*m1[idx])*(m3[idx]-4*m2[idx]+3*m1[idx]);  \
    b2 = thirteen_by_twelve*(m2[idx]-2*m1[idx]+p1[idx])*(m2[idx]-2*m1[idx]+p1[idx]) \
          + one_fourth*(m2[idx]-p1[idx])*(m2[idx]-p1[idx]);                \
    b3 = thirteen_by_twelve*(m1[idx]-2*p1[idx]+p2[idx])*(m1[idx]-2*p1[idx]+p2[idx]) \
          + one_fourth*(3*m1[idx]-4*p1[idx]+p2[idx])*(3*m1[idx]-4*p1[idx]+p2[idx]);  \
    tau = (m3[idx]-4*m2[idx]+6*m1[idx]-4*p1[idx]+p2[idx])*(m3[idx]-4*m2[idx]+6*m1[idx]-4*p1[idx]+p2[idx]);  \
    a1 = c1 * (1.0 + (tau/(b1+eps)) * (tau/(b1+eps)) );  \
    a2 = c2 * (1.0 + (tau/(b2+eps)) * (tau/(b2+eps)) );  \
    a3 = c3 * (1.0 + (tau/(b3+eps)) * (tau/(b3+eps)) );  \
    a_sum_inv = 1.0 / (a1 + a2 + a3); \
    w1[idx] = a1 * a_sum_inv; \
    w2[idx] = a2 * a_sum_inv; \
    w3[idx] = a3 * a_sum_inv; \
  }

Compute the WENO weights according the the ESWENO formulation of Yamaleev & Carpenter. Note that only the formulation for the nonlinear weights is adopted and implemented here, not the ESWENO scheme as a whole.

\begin{equation} \omega_k = \frac {a_k} {\sum_{j=1}^3 a_j },\ a_k = c_k \left( 1 + \frac{\tau_5}{\beta_k+\epsilon} \right)^p,\ k = 1,2,3, \end{equation}

where \(c_k\) are the optimal weights, \(p\) is hardcoded to \(2\), and \(\epsilon\) is an input parameter (WENOParameters::eps) (typically \(10^{-6}\)). The smoothness indicators \(\beta_k\) are given by:

\begin{eqnarray} \beta_1 &=& \frac{13}{12} \left(f_{j-2}-2f_{j-1}+f_j\right)^2 + \frac{1}{4}\left(f_{j-2}-4f_{j-1}+3f_j\right)^2 \\ \beta_2 &=& \frac{13}{12} \left(f_{j-1}-2f_j+f_{j+1}\right)^2 + \frac{1}{4}\left(f_{j-1}-f_{j+1}\right)^2 \\ \beta_3 &=& \frac{13}{12} \left(f_j-2f_{j+1}+f_{j+2}\right)^2 + \frac{1}{4}\left(3f_j-4f_{j+1}+f_{j+2}\right)^2, \end{eqnarray}

and \(\tau_5 = \left( f_{j-2}-4f_{j-1}+6f_j-4f_{j+1}+f_{j+2} \right)^2\).

Notes:

This macro computes the weights for one variable along one grid line.

Arguments:

w1, w2,w3 are the nonlinear WENO weights.
c1, c2,c3 are optimal coefficients.
m3, m2,m1,p1,p2 are the function values at stencil points corresponding to the interface j+1/2: j-2,j-1,j,j+1,j+2
weno is an object of type WENOParameters containing parameters for the WENO method.
N is the number of interfaces along the grid line on which this WENO-type reconstruction is happening.
Reference:
- Yamaleev, Carpenter, A systematic methodology for constructing high-order energy stable WENO schemes, J. Comput. Phys., 2009. http://dx.doi.org/10.1016/j.jcp.2009.03.002

Definition at line 532 of file interpolation.h.

Function Documentation

◆ Interp1PrimFirstOrderUpwind()

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

1st order upwind reconstruction (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the first-order upwind scheme

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

\begin{equation} \hat{\bf f}_{j+1/2} = \left\{\begin{array}{cc} {\bf f}_{j} & {\rm upw} > 0 \\ {\bf f}_{j+1} & {\rm upw} \le 0 \end{array}\right.. \end{equation}

Implementation Notes:

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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 61 of file Interp1PrimFirstOrderUpwind.c.

 {
   HyPar         *solver = (HyPar*)        s;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   int i;
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       indexC[dir] = (upw > 0 ? indexI[dir]-1 : indexI[dir]);
       int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0     ,p);
       int q; _ArrayIndex1D_(ndims,dim         ,indexC,ghosts,q);
       int v; for (v=0; v<nvars; v++)  fI[p*nvars+v] = fC[q*nvars+v];
     }
   }
 
   return(0);
 }

◆ Interp1PrimSecondOrderCentral()

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

2nd order central reconstruction (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the second-order central scheme

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

\begin{equation} \hat{\bf f}_{j+1/2} = \frac{1}{2}\left( {\bf f}_{j} + {\bf f}_{j+1} \right). \end{equation}

Implementation Notes:

The scalar interpolation method is applied to the vector function in a component-wise manner.
Since this is a central scheme, the input argument upw has no effect.
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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 62 of file Interp1PrimSecondOrderCentral.c.

 {
   HyPar         *solver = (HyPar*)        s;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexL[ndims], indexR[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   int i;
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexL,indexR,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexL,ndims);
     _ArrayCopy1D_(index_outer,indexR,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       indexL[dir] = indexI[dir]-1;
       indexR[dir] = indexI[dir];
       int p;  _ArrayIndex1D_(ndims,bounds_inter,indexI,0     ,p);
       int qL; _ArrayIndex1D_(ndims,dim         ,indexL,ghosts,qL);
       int qR; _ArrayIndex1D_(ndims,dim         ,indexR,ghosts,qR);
       int v; for (v=0; v<nvars; v++)  fI[p*nvars+v] = 0.5*(fC[qL*nvars+v]+fC[qR*nvars+v]);
     }
   }
 
   return(0);
 }

◆ Interp1PrimSecondOrderMUSCL()

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

Component-wise interpolation of the first primitive at the cell interfaces using the second-order MUSCL scheme

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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 76 of file Interp1PrimSecondOrderMUSCL.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MUSCLParameters *muscl  = (MUSCLParameters*) solver->interp;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   double half = 1.0/2.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   int i;
   if (upw > 0) {
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
         int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
         int qm2,qm1,qp1,v;
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         for (v=0; v<nvars; v++)  {
           /* Defining stencil points */
           double m2, m1, p1;
           m2 = fC[qm2*nvars+v];
           m1 = fC[qm1*nvars+v];
           p1 = fC[qp1*nvars+v];
 
           double slope_ratio = (m1 - m2) / ((p1 - m1) + 1e-40);
           double phi = muscl->LimiterFunction(slope_ratio);
           fI[p*nvars+v] =  m1 + 0.5 * phi * (p1-m1);
         }
       }
     }
   } else {
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
         int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
         int qm1,qp1,qp2,v;
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
         for (v=0; v<nvars; v++)  {
           /* Defining stencil points */
           double m1, p1, p2;
           m1 = fC[qm1*nvars+v];
           p1 = fC[qp1*nvars+v];
           p2 = fC[qp2*nvars+v];
 
           double slope_ratio = (p1 - m1) / ((p2 - p1) + 1e-40);
           double phi = muscl->LimiterFunction(slope_ratio);
           fI[p*nvars+v] = p1 + 0.5 * phi * (p1-p2);
         }
       }
     }
   }
   
   return(0);
 }

◆ Interp1PrimThirdOrderMUSCL()

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

3rd order MUSCL scheme with Koren's limiter (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the third-order MUSCL scheme

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 3rd order MUSCL scheme with Koren's limiter 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 3rd order MUSCL scheme with Koren's limiter as follows:

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

where

\begin{equation} \phi = \frac {3\left({\bf f}_j-{\bf f}_{j-1}\right)\left({\bf f}_{j-1}-{\bf f}_{j-2}\right) + \epsilon} {2\left[\left({\bf f}_j-{\bf f}_{j-1}\right)-\left({\bf f}_{j-1}-{\bf f}_{j-2}\right)\right]^2 + 3\left({\bf f}_j-{\bf f}_{j-1}\right)\left({\bf f}_{j-1}-{\bf f}_{j-2}\right) + \epsilon}. \end{equation}

and \(\epsilon\) is a small constant (typically \(10^{-3}\)).

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:

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 (nvarsdim[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., Towards the Ultimate Conservative Difference Scheme. 2: Monotonicity and Conservation Combined in a Second-Order Scheme, J. of Comput. Phys., 14 (4), 1974, pp.361-370, http://dx.doi.org/10.1016/0021-9991(74)90019-9
Koren, B., A Robust Upwind Discretization Method for Advection, Diffusion and Source Terms, Centrum voor Wiskunde en Informatica, Amsterdam, 1993

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 75 of file Interp1PrimThirdOrderMUSCL.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MUSCLParameters *muscl  = (MUSCLParameters*) solver->interp;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   double one_third = 1.0/3.0;
   double one_sixth = 1.0/6.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   int i;
   if (upw > 0) {
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
         int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
         int qm1,qm2,qp1,v;
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         for (v=0; v<nvars; v++)  {
           /* Defining stencil points */
           double m2, m1, p1;
           m2 = fC[qm2*nvars+v];
           m1 = fC[qm1*nvars+v];
           p1 = fC[qp1*nvars+v];
 
           double fdiff = p1 - m1;
           double bdiff = m1 - m2;
           double limit =  (3*fdiff*bdiff + muscl->eps) 
                         / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
 
           fI[p*nvars+v] = m1 +  limit * (one_third*fdiff + one_sixth*bdiff);
         }
       }
     }
   } else {
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
         int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
         int qm1,qp1,qp2,v;
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
         for (v=0; v<nvars; v++)  {
           /* Defining stencil points */
           double m1, p1, p2;
           m1 = fC[qm1*nvars+v];
           p1 = fC[qp1*nvars+v];
           p2 = fC[qp2*nvars+v];
 
           double fdiff = p2 - p1;
           double bdiff = p1 - m1;
           double limit =  (3*fdiff*bdiff + muscl->eps) 
                         / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
 
           fI[p*nvars+v] = p1 -  limit * (one_third*fdiff + one_sixth*bdiff);
         }
       }
     }
   }
   
   return(0);
 }

◆ Interp1PrimFourthOrderCentral()

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

4th order central reconstruction (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the fourth-order central scheme

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{\bf f}_{j+1/2} = -\frac{1}{12} {\bf f}_{j-1} + \frac{7}{12} {\bf f}_{j} + \frac{7}{12} {\bf f}_{j+1} - \frac{1}{12} {\bf f}_{j+2}. \end{equation}

Implementation Notes:

The scalar interpolation method is applied to the vector function in a component-wise manner.
Since this is a central scheme, the input argument upw has no effect.
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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 62 of file Interp1PrimFourthOrderCentral.c.

 {
   HyPar         *solver = (HyPar*)        s;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexL[ndims], 
       indexR[ndims], 
       indexLL[ndims], 
       indexRR[ndims], 
       indexI[ndims], 
       index_outer[ndims], 
       bounds_outer[ndims], 
       bounds_inter[ndims];
 
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   static const double c1 = 7.0 / 12.0;
   static const double c2 = -1.0 / 12.0;
 
   int i;
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexL,indexR,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexL,ndims);
     _ArrayCopy1D_(index_outer,indexLL,ndims);
     _ArrayCopy1D_(index_outer,indexR,ndims);
     _ArrayCopy1D_(index_outer,indexRR,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       indexLL[dir] = indexI[dir]-2;
       indexL[dir]  = indexI[dir]-1;
       indexR[dir]  = indexI[dir];
       indexRR[dir] = indexI[dir]+1;
 
       int p;    _ArrayIndex1D_(ndims,bounds_inter,indexI ,0     ,p);
       int qLL;  _ArrayIndex1D_(ndims,dim         ,indexLL,ghosts,qLL);
       int qL;   _ArrayIndex1D_(ndims,dim         ,indexL ,ghosts,qL );
       int qR;   _ArrayIndex1D_(ndims,dim         ,indexR ,ghosts,qR );
       int qRR;  _ArrayIndex1D_(ndims,dim         ,indexRR,ghosts,qRR);
 
       int v; for (v=0; v<nvars; v++)  {
         fI[p*nvars+v] =   c2 * fC[qLL*nvars+v]
                         + c1 * fC[qL *nvars+v]
                         + c1 * fC[qR *nvars+v]
                         + c2 * fC[qRR*nvars+v];
       }
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderUpwind()

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

5th order upwind reconstruction (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the fifth-order upwind scheme

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

\begin{align} \hat{\bf f}_{j+1/2} = \frac{1}{30} {\bf f}_{j-2} - \frac{13}{60}{\bf f}_{j-1} + \frac{47}{60}{\bf f}_j + \frac{27}{60}{\bf f}_{j+1} - \frac{1}{20}{\bf f}_{j+2}. \end{align}

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:

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 (nvarsdim[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:

Jiang, G.-S., Shu, C.-W., Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126 (1), 1996, pp. 202-228, http://dx.doi.org/10.1006/jcph.1996.0130

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 68 of file Interp1PrimFifthOrderUpwind.c.

 {
   HyPar           *solver = (HyPar*)          s;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
   int *stride= solver->stride_with_ghosts;
 
   /* define some constants */
   static const double one_by_thirty        = 1.0/30.0,
                       thirteen_by_sixty    = 13.0/60.0,
                       fortyseven_by_sixty  = 47.0/60.0,
                       twentyseven_by_sixty = 27.0/60.0,
                       one_by_twenty        = 1.0/20.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   int i;
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2,p;
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       if (upw > 0) {
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 - 2*stride[dir];
         qm2 = qm1 -   stride[dir];
         qp1 = qm1 +   stride[dir];
         qp2 = qm1 + 2*stride[dir];
       } else {
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 + 2*stride[dir];
         qm2 = qm1 +   stride[dir];
         qp1 = qm1 -   stride[dir];
         qp2 = qm1 - 2*stride[dir];
       }
 
       /* Defining stencil points */
       double *fm3, *fm2, *fm1, *fp1, *fp2;
       fm3 = (fC+qm3*nvars);
       fm2 = (fC+qm2*nvars);
       fm1 = (fC+qm1*nvars);
       fp1 = (fC+qp1*nvars);
       fp2 = (fC+qp2*nvars);
 
       int v;
       for (v=0; v<nvars; v++) {
         (fI+p*nvars)[v] =   one_by_thirty         * fm3[v]
                           - thirteen_by_sixty     * fm2[v]
                           + fortyseven_by_sixty   * fm1[v]
                           + twentyseven_by_sixty  * fp1[v]
                           - one_by_twenty         * fp2[v];
       }
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderCompactUpwind()

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

5th order compact upwind reconstruction (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the fifth-order compact upwind scheme

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

\begin{align} \frac{3}{10}\hat{\bf f}_{j-1/2} + \frac{6}{10}\hat{\bf f}_{j+1/2} + \frac{1}{10}\hat{\bf f}_{j+3/2} = \frac{1}{30}{\bf f}_{j-1} + \frac{19}{30}{\bf f}_j + \frac{1}{3}{\bf f}_{j+1}. \end{align}

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:

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 (nvarsdim[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

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 72 of file Interp1PrimFifthOrderCompactUpwind.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   CompactScheme   *compact= (CompactScheme*)  solver->compact;
   TridiagLU       *lu     = (TridiagLU*)      solver->lusolver;
   int             sys,Nsys,d,v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
   int *stride= solver->stride_with_ghosts;
 
   /* define some constants */
   static const double one_third            = 1.0/3.0,
                       thirteen_by_sixty    = 13.0/60.0,
                       fortyseven_by_sixty  = 47.0/60.0,
                       twentyseven_by_sixty = 27.0/60.0,
                       one_by_twenty        = 1.0/20.0,
                       one_by_thirty        = 1.0/30.0,
                       nineteen_by_thirty   = 19.0/30.0,
                       three_by_ten         = 3.0/10.0,
                       six_by_ten           = 6.0/10.0,
                       one_by_ten           = 1.0/10.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* calculate total number of tridiagonal systems to solve */
   _ArrayProduct1D_(bounds_outer,ndims,Nsys); Nsys *= nvars;
 
   /* Allocate arrays for tridiagonal system */
   double *A = compact->A;
   double *B = compact->B;
   double *C = compact->C;
   double *R = compact->R;
 
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
   for (sys=0; sys < N_outer; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2,p;
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       if (upw > 0) {
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 - 2*stride[dir];
         qm2 = qm1 -   stride[dir];
         qp1 = qm1 +   stride[dir];
         qp2 = qm1 + 2*stride[dir];
       } else {
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 + 2*stride[dir];
         qm2 = qm1 +   stride[dir];
         qp1 = qm1 -   stride[dir];
         qp2 = qm1 - 2*stride[dir];
       }
 
       /* Defining stencil points */
       double *fm3, *fm2, *fm1, *fp1, *fp2;
       fm3 = fC+qm3*nvars;
       fm2 = fC+qm2*nvars;
       fm1 = fC+qm1*nvars;
       fp1 = fC+qp1*nvars;
       fp2 = fC+qp2*nvars;
 
       if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
           || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
 
         /* Use 5th order upwind at the physical boundaries */
         _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
         _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,1.0)
         _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
         for (v=0; v<nvars; v++) {
           (R+Nsys*indexI[dir]+sys*nvars)[v] =   one_by_thirty         * fm3[v]
                                               - thirteen_by_sixty     * fm2[v]
                                               + fortyseven_by_sixty   * fm1[v]
                                               + twentyseven_by_sixty  * fp1[v]
                                               - one_by_twenty         * fp2[v];
         }
 
       } else {
 
         /* Use 5th order upwind at the physical boundaries */
         if (upw > 0) {
           _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,three_by_ten);
           _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,six_by_ten  );
           _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,one_by_ten  );
         } else {
           _ArraySetValue_((C+Nsys*indexI[dir]+sys*nvars),nvars,three_by_ten);
           _ArraySetValue_((B+Nsys*indexI[dir]+sys*nvars),nvars,six_by_ten  );
           _ArraySetValue_((A+Nsys*indexI[dir]+sys*nvars),nvars,one_by_ten  );
         }
         for (v=0; v<nvars; v++) {
           (R+Nsys*indexI[dir]+sys*nvars)[v] =   one_by_thirty      * fm2[v]
                                               + nineteen_by_thirty * fm1[v]
                                               + one_third          * fp1[v];
         }
       }
     }
   }
 
 #ifdef serial
 
   /* Solve the tridiagonal system */
   IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,NULL); CHECKERR(ierr);
 
 #else
 
   /* Solve the tridiagonal system */
   /* all processes except the last will solve without the last interface to avoid overlap */
   if (mpi->ip[dir] != mpi->iproc[dir]-1)  { IERR tridiagLU(A,B,C,R,dim[dir]  ,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
   else                                    { IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
 
   /* Now get the solution to the last interface from the next proc */
   double *sendbuf = compact->sendbuf;
   double *recvbuf = compact->recvbuf;
   MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
   if (mpi->ip[dir]) for (d=0; d<Nsys; d++) sendbuf[d] = R[d];
   if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
   if (mpi->ip[dir])                      MPI_Isend(sendbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
   MPI_Status status_arr[2];
   MPI_Waitall(2,&req[0],status_arr);
   if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys; d++) R[d+Nsys*dim[dir]] = recvbuf[d];
 
 #endif
 
   /* save the solution to fI */
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
   for (sys=0; sys < N_outer; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       _ArrayCopy1D_((R+sys*nvars+Nsys*indexI[dir]),(fI+nvars*p),nvars);
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderWENO()

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

5th order WENO reconstruction (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the fifth-order WENO scheme

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 WENO 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 WENO 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[ \hat{\bf f}_{j+1/2}^1 = \frac{1}{3} {\bf f}_{j-2} - \frac{7}{6} {\bf f}_{j-1} + \frac{11}{6} {\bf f}_j \right]\\ + &\ \omega_2\ \times\ \left[ \hat{\bf f}_{j+1/2}^2 = -\frac{1}{6} {\bf f}_{j-1} + \frac{5}{6} {\bf f}_j + \frac{1}{3} {\bf f}_{j+1} \right]\\ + &\ \omega_3\ \times\ \left[ \hat{\bf f}_{j+1/2}^3 = \frac{1}{3} {\bf f}_j + \frac{5}{6} {\bf f}_{j+1} - \frac{1}{6} {\bf f}_{j+2} \right]\\ \Rightarrow &\ \hat{\bf f}_{j+1/2} = \frac{\omega_1}{3} {\bf f}_{j-2} - \frac{1}{6}(7\omega_1+\omega_2){\bf f}_{j-1} + \frac{1}{6}(11\omega_1+5\omega_2+2\omega_3){\bf f}_j + \frac{1}{6}(2\omega_2+5\omega_3){\bf f}_{j+1} - \frac{\omega_3}{6}{\bf f}_{j+2}, \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}\)).

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:

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 (nvarsdim[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:

Jiang, G.-S., Shu, C.-W., Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126 (1), 1996, pp. 202-228, http://dx.doi.org/10.1006/jcph.1996.0130

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 74 of file Interp1PrimFifthOrderWENO.c.

 {
   HyPar           *solver = (HyPar*)          s;
   WENOParameters  *weno   = (WENOParameters*) solver->interp;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
   int *stride= solver->stride_with_ghosts;
 
   /* define some constants */
   static const double one_sixth          = 1.0/6.0;
 
   double *ww1, *ww2, *ww3;
   ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
 #if defined(CPU_STAT)
   clock_t cpu_start, cpu_end;
   cpu_start = clock();
 #endif
 
   int i;
 #pragma omp parallel for schedule(auto) default(shared) private(i,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2,p;
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       if (upw > 0) {
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 - 2*stride[dir];
         qm2 = qm1 -   stride[dir];
         qp1 = qm1 +   stride[dir];
         qp2 = qm1 + 2*stride[dir];
       } else {
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 + 2*stride[dir];
         qm2 = qm1 +   stride[dir];
         qp1 = qm1 -   stride[dir];
         qp2 = qm1 - 2*stride[dir];
       }
 
       /* Defining stencil points */
       double *fm3, *fm2, *fm1, *fp1, *fp2;
       fm3 = (fC+qm3*nvars);
       fm2 = (fC+qm2*nvars);
       fm1 = (fC+qm1*nvars);
       fp1 = (fC+qp1*nvars);
       fp2 = (fC+qp2*nvars);
 
       /* Candidate stencils and their optimal weights*/
       double f1[nvars], f2[nvars], f3[nvars];
       _ArrayAXBYCZ_(f1,(2*one_sixth),fm3,(-7*one_sixth) ,fm2,(11*one_sixth) ,fm1,nvars);
       _ArrayAXBYCZ_(f2,(-one_sixth) ,fm2,(5*one_sixth)  ,fm1,(2*one_sixth)  ,fp1,nvars);
       _ArrayAXBYCZ_(f3,(2*one_sixth),fm1,(5*one_sixth)  ,fp1,(-one_sixth)   ,fp2,nvars);
 
       /* calculate WENO weights */
       double *w1,*w2,*w3;
       w1 = (ww1+p*nvars);
       w2 = (ww2+p*nvars);
       w3 = (ww3+p*nvars);
 
       _ArrayMultiply3Add1D_((fI+p*nvars),w1,f1,w2,f2,w3,f3,nvars);
     }
   }
 
 #if defined(CPU_STAT)
   cpu_end = clock();
   printf("Interp1PrimFifthOrderWENO CPU time = %8.6lf dir = %d\n", (double)(cpu_end - cpu_start) / CLOCKS_PER_SEC, dir);
 #endif
 
   return(0);
 }

◆ 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

Component-wise interpolation of the first primitive at the cell interfaces using the fifth-order CRWENO scheme

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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 78 of file Interp1PrimFifthOrderCRWENO.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   CompactScheme   *compact= (CompactScheme*)  solver->compact;
   WENOParameters  *weno   = (WENOParameters*) solver->interp;
   TridiagLU       *lu     = (TridiagLU*)      solver->lusolver;
   int             sys,Nsys,d;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
   int *stride= solver->stride_with_ghosts;
 
   /* define some constants */
   static const double one_third          = 1.0/3.0;
   static const double one_sixth          = 1.0/6.0;
 
   double *ww1, *ww2, *ww3;
   ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* calculate total number of tridiagonal systems to solve */
   _ArrayProduct1D_(bounds_outer,ndims,Nsys); Nsys *= nvars;
 
   /* Allocate arrays for tridiagonal system */
   double *A = compact->A;
   double *B = compact->B;
   double *C = compact->C;
   double *R = compact->R;
 
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
   for (sys=0; sys < N_outer; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2,p;
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       if (upw > 0) {
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 - 2*stride[dir];
         qm2 = qm1 -   stride[dir];
         qp1 = qm1 +   stride[dir];
         qp2 = qm1 + 2*stride[dir];
       } else {
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         qm3 = qm1 + 2*stride[dir];
         qm2 = qm1 +   stride[dir];
         qp1 = qm1 -   stride[dir];
         qp2 = qm1 - 2*stride[dir];
       }
 
       /* Defining stencil points */
       double *fm3, *fm2, *fm1, *fp1, *fp2;
       fm3 = fC+qm3*nvars;
       fm2 = fC+qm2*nvars;
       fm1 = fC+qm1*nvars;
       fp1 = fC+qp1*nvars;
       fp2 = fC+qp2*nvars;
 
       /* Candidate stencils and their optimal weights*/
       double f1[nvars], f2[nvars], f3[nvars];
       if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
           || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
         /* Use WENO5 at the physical boundaries */
         _ArrayAXBYCZ_(f1,(2*one_sixth),fm3,(-7*one_sixth) ,fm2,(11*one_sixth) ,fm1,nvars);
         _ArrayAXBYCZ_(f2,(-one_sixth) ,fm2,(5*one_sixth)  ,fm1,(2*one_sixth)  ,fp1,nvars);
         _ArrayAXBYCZ_(f3,(2*one_sixth),fm1,(5*one_sixth)  ,fp1,(-one_sixth)   ,fp2,nvars);
       } else {
         /* CRWENO5 at the interior points */
         _ArrayAXBY_(f1,(one_sixth)  ,fm2,(5*one_sixth),fm1,nvars);
         _ArrayAXBY_(f2,(5*one_sixth),fm1,(one_sixth)  ,fp1,nvars);
         _ArrayAXBY_(f3,(one_sixth)  ,fm1,(5*one_sixth),fp1,nvars);
       }
 
       /* retrieve the WENO weights */
       double *w1, *w2, *w3;
       w1 = (ww1+p*nvars);
       w2 = (ww2+p*nvars);
       w3 = (ww3+p*nvars);
 
       if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
           || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
         _ArraySetValue_     ((A+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
         _ArraySetValue_     ((B+Nsys*indexI[dir]+sys*nvars),nvars,1.0)
         _ArraySetValue_     ((C+Nsys*indexI[dir]+sys*nvars),nvars,0.0)
       } else {
         if (upw > 0) {
           _ArrayAXBY_       ((A+Nsys*indexI[dir]+sys*nvars),(2*one_third) ,w1,(one_third)  ,w2,nvars);
           _ArrayAXBYCZ_     ((B+Nsys*indexI[dir]+sys*nvars),(one_third)   ,w1,(2*one_third),w2,(2*one_third),w3,nvars);
           _ArrayScaleCopy1D_(w3,(one_third),(C+Nsys*indexI[dir]+sys*nvars),nvars);
         } else {
           _ArrayAXBY_       ((C+Nsys*indexI[dir]+sys*nvars),(2*one_third) ,w1,(one_third)  ,w2,nvars);
           _ArrayAXBYCZ_     ((B+Nsys*indexI[dir]+sys*nvars),(one_third)   ,w1,(2*one_third),w2,(2*one_third),w3,nvars);
           _ArrayScaleCopy1D_(w3,(one_third),(A+Nsys*indexI[dir]+sys*nvars),nvars);
         }
       }
       _ArrayMultiply3Add1D_ ((R+Nsys*indexI[dir]+sys*nvars),w1,f1,w2,f2,w3,f3,nvars);
     }
   }
 
 #ifdef serial
 
   /* Solve the tridiagonal system */
   IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,NULL); CHECKERR(ierr);
 
 #else
 
   /* Solve the tridiagonal system */
   /* all processes except the last will solve without the last interface to avoid overlap */
   if (mpi->ip[dir] != mpi->iproc[dir]-1)  { IERR tridiagLU(A,B,C,R,dim[dir]  ,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
   else                                    { IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
 
   /* Now get the solution to the last interface from the next proc */
   double *sendbuf = compact->sendbuf;
   double *recvbuf = compact->recvbuf;
   MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
   if (mpi->ip[dir]) for (d=0; d<Nsys; d++) sendbuf[d] = R[d];
   if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
   if (mpi->ip[dir])                      MPI_Isend(sendbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
   MPI_Status status_arr[2];
   MPI_Waitall(2,&req[0],status_arr);
   if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys; d++) R[d+Nsys*dim[dir]] = recvbuf[d];
 
 #endif
 
   /* save the solution to fI */
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
   for (sys=0; sys < N_outer; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       _ArrayCopy1D_((R+sys*nvars+Nsys*indexI[dir]),(fI+nvars*p),nvars);
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderHCWENO()

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

5th order hybrid compact-WENO reconstruction (component-wise) on a uniform grid

Component-wise interpolation of the first primitive at the cell interfaces using the fifth-order hybrid-compact WENO scheme

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 hybrid compact-WENO scheme on a uniform grid. The tridiagonal system is solved using tridiagLU() (see also TridiagLU, tridiagLU.h). See references below for a complete description of the method implemented here.

Implementation Notes:

This method assumes a uniform grid in the spatial dimension corresponding to the interpolation.
The scalar interpolation method is applied to the vector function in a component-wise manner.
The WENO weights are computed in WENOFifthOrderCalculateWeights().
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:

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 (nvarsdim[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:

Pirozzoli, S., Conservative Hybrid Compact-WENO Schemes for Shock-Turbulence Interaction, J. Comput. Phys., 178 (1), 2002, pp. 81-117, http://dx.doi.org/10.1006/jcph.2002.7021
Ren, Y.-X., Liu, M., Zhang, H., A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws, J. Comput. Phys., 192 (2), 2003, pp. 365-386, http://dx.doi.org/10.1016/j.jcp.2003.07.006

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 63 of file Interp1PrimFifthOrderHCWENO.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   CompactScheme   *compact= (CompactScheme*)  solver->compact;
   WENOParameters  *weno   = (WENOParameters*) solver->interp;
   TridiagLU       *lu     = (TridiagLU*)      solver->lusolver;
   int             sys,Nsys,d;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   static const double one_half           = 1.0/2.0;
   static const double one_third          = 1.0/3.0;
   static const double one_sixth          = 1.0/6.0;
 
   double *ww1, *ww2, *ww3;
   ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* calculate total number of tridiagonal systems to solve */
   _ArrayProduct1D_(bounds_outer,ndims,Nsys); Nsys *= nvars;
 
   /* Allocate arrays for tridiagonal system */
   double *A = compact->A;
   double *B = compact->B;
   double *C = compact->C;
   double *R = compact->R;
 
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
   for (sys=0; sys < N_outer; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2,p;
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       if (upw > 0) {
         indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       } else {
         indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       }
       int v; 
       for (v=0; v<nvars; v++)  {
         /* Defining stencil points */
         double fm3, fm2, fm1, fp1, fp2;
         fm3 = fC[qm3*nvars+v];
         fm2 = fC[qm2*nvars+v];
         fm1 = fC[qm1*nvars+v];
         fp1 = fC[qp1*nvars+v];
         fp2 = fC[qp2*nvars+v];
 
         /* Candidate stencils and their optimal weights*/
         double f1, f2, f3;
         f1 = (2*one_sixth)*fm3 - (7.0*one_sixth)*fm2 + (11.0*one_sixth)*fm1;
         f2 = (-one_sixth)*fm2 + (5.0*one_sixth)*fm1 + (2*one_sixth)*fp1;
         f3 = (2*one_sixth)*fm1 + (5*one_sixth)*fp1 - (one_sixth)*fp2;
 
         /* calculate WENO weights */
         double w1,w2,w3;
         w1 = *(ww1+p*nvars+v);
         w2 = *(ww2+p*nvars+v);
         w3 = *(ww3+p*nvars+v);
 
         /* calculate the hybridization parameter */
         double sigma;
         if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
             || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
           /* Standard WENO at physical boundaries */
           sigma = 0.0;
         } else {
           double cuckoo, df_jm12, df_jp12, df_jp32, r_j, r_jp1, r_int;
           cuckoo = (0.9*weno->rc / (1.0-0.9*weno->rc)) * weno->xi * weno->xi;
           df_jm12 = fm1 - fm2;
           df_jp12 = fp1 - fm1;
           df_jp32 = fp2 - fp1;
           r_j   = (absolute(2*df_jp12*df_jm12)+cuckoo)/(df_jp12*df_jp12+df_jm12*df_jm12+cuckoo);
           r_jp1 = (absolute(2*df_jp32*df_jp12)+cuckoo)/(df_jp32*df_jp32+df_jp12*df_jp12+cuckoo);
           r_int = min(r_j, r_jp1);
           sigma = min((r_int/weno->rc), 1.0); 
         }
 
         if (upw > 0) {
           A[sys*nvars+v+Nsys*indexI[dir]] = one_half  * sigma;
           B[sys*nvars+v+Nsys*indexI[dir]] = 1.0;
           C[sys*nvars+v+Nsys*indexI[dir]] = one_sixth * sigma;
         } else {
           C[sys*nvars+v+Nsys*indexI[dir]] = one_half  * sigma;
           B[sys*nvars+v+Nsys*indexI[dir]] = 1.0;
           A[sys*nvars+v+Nsys*indexI[dir]] = one_sixth * sigma;
         }
 
         double fWENO, fCompact;
         fWENO    = w1*f1 + w2*f2 + w3*f3;
         fCompact = one_sixth * (one_third*fm2 + 19.0*one_third*fm1 + 10.0*one_third*fp1);
         R[sys*nvars+v+Nsys*indexI[dir]] =   sigma*fCompact + (1.0-sigma)*fWENO;
       }
     }
   }
 
 #ifdef serial
 
   /* Solve the tridiagonal system */
   IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,NULL); CHECKERR(ierr);
 
 #else
 
   /* Solve the tridiagonal system */
   /* all processes except the last will solve without the last interface to avoid overlap */
   if (mpi->ip[dir] != mpi->iproc[dir]-1)  { IERR tridiagLU(A,B,C,R,dim[dir]  ,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
   else                                    { IERR tridiagLU(A,B,C,R,dim[dir]+1,Nsys,lu,&mpi->comm[dir]); CHECKERR(ierr); }
 
   /* Now get the solution to the last interface from the next proc */
   double *sendbuf = compact->sendbuf;
   double *recvbuf = compact->recvbuf;
   MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
   if (mpi->ip[dir]) for (d=0; d<Nsys; d++) sendbuf[d] = R[d];
   if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
   if (mpi->ip[dir])                      MPI_Isend(sendbuf,Nsys,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
   MPI_Status status_arr[2];
   MPI_Waitall(2,&req[0],status_arr);
   if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys; d++) R[d+Nsys*dim[dir]] = recvbuf[d];
 
 #endif
 
   /* save the solution to fI */
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,index_outer,indexC,indexI)
   for (sys=0; sys < N_outer; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       _ArrayCopy1D_((R+sys*nvars+Nsys*indexI[dir]),(fI+nvars*p),nvars);
     }
   }
 
   return(0);
 }

◆ gpuInterp1PrimFifthOrderWENO()

int gpuInterp1PrimFifthOrderWENO	(	double *	,
		double *	,
		double *	,
		double *	,
		int	,
		int	,
		void *	,
		void *	,
		int
	)

GPU implementation of component-wise interpolation of the first primitive at the cell interfaces using the fifth-order WENO scheme

◆ Interp1PrimFirstOrderUpwindChar()

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

1st order upwind reconstruction (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the first-order upwind scheme

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

\begin{equation} \hat{\alpha}^k_{j+1/2} = \left\{\begin{array}{cc} {\alpha}^k_{j} & {\rm upw} > 0 \\ {\alpha}^k_{j+1} & {\rm upw} \le 0 \end{array}\right.. \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:

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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 77 of file Interp1PrimFirstOrderUpwindChar.c.

 {
   HyPar         *solver = (HyPar*) s;
   int           i, k, v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
 
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
 
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
       int pL, pR; /* 1D index of the left and right cells */
       indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pL);
       indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pR);
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*pL],&u[nvars*pR],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       /* For each characteristic field */
       for (v = 0; v < nvars; v++) {
 
         /* calculate the characteristic flux components along this characteristic */
         double fcL = 0, fcR = 0;
         for (k = 0; k < nvars; k++) {
           fcL += L[v*nvars+k] * fC[pL*nvars+k];
           fcR += L[v*nvars+k] * fC[pR*nvars+k];
         }
 
         /* first order upwind approximation of the characteristic flux */
         fchar[v] = (upw > 0 ? fcL : fcR);
 
       }
 
       /* calculate the interface u from the characteristic u */
       IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
 
     }
   }
 
   return(0);
 }

◆ Interp1PrimSecondOrderCentralChar()

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

2nd order central reconstruction (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the second-order central scheme

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 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 2nd 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}{2} \left( {\alpha}^k_{j+1} + {\alpha}^k_j \right), \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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 79 of file Interp1PrimSecondOrderCentralChar.c.

 {
   HyPar         *solver = (HyPar*) s;
   int           i, k, v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
 
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
 
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
       int pL, pR; /* 1D index of the left and right cells */
       indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pL);
       indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pR);
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*pL],&u[nvars*pR],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       /* For each characteristic field */
       for (v = 0; v < nvars; v++) {
 
         /* calculate the characteristic flux components along this characteristic */
         double fcL = 0, fcR = 0;
         for (k = 0; k < nvars; k++) {
           fcL += L[v*nvars+k] * fC[pL*nvars+k];
           fcR += L[v*nvars+k] * fC[pR*nvars+k];
         }
 
         /* first order upwind approximation of the characteristic flux */
         fchar[v] = 0.5 * (fcL + fcR);
 
       }
 
       /* calculate the interface u from the characteristic u */
       IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
 
     }
   }
 
   return(0);
 }

◆ 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

Characteristic-based interpolation of the first primitive at the cell interfaces using the second-order MUSCL scheme

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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 92 of file Interp1PrimSecondOrderMUSCLChar.c.

 {
   HyPar             *solver = (HyPar*) s;
   MUSCLParameters   *muscl   = (MUSCLParameters*) solver->interp;
   int               i, k, v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   double one_third = 1.0/3.0;
   double one_sixth = 1.0/6.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
 
   if (upw > 0) {
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
 
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
         /* 1D indices of the stencil grid points */
         int qm1,qm2,qp1;
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
 
         int p; /* 1D index of the interface */
         _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
         /* find averaged state at this interface */
         IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
         CHECKERR(ierr);
 
         /* Get the left and right eigenvectors */
         IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
         IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
         /* For each characteristic field */
         for (v = 0; v < nvars; v++) {
           /* calculate the characteristic flux components along this characteristic */
           double m2, m1, p1;
           m2 = m1 = p1 = 0;
           for (k = 0; k < nvars; k++) {
             m2 += L[v*nvars+k] * fC[qm2*nvars+k];
             m1 += L[v*nvars+k] * fC[qm1*nvars+k];
             p1 += L[v*nvars+k] * fC[qp1*nvars+k];
           }
 
           double slope_ratio = (m1 - m2) / ((p1 - m1) + 1e-40);
           double phi = muscl->LimiterFunction(slope_ratio);
           fchar[v] = m1 + 0.5 * phi * (p1-m1);
         }
 
         /* calculate the interface u from the characteristic u */
         IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
       }
     }
   } else {
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
 
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
         /* 1D indices of the stencil grid points */
         int qm1,qp1,qp2;
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
 
         int p; /* 1D index of the interface */
         _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
         /* find averaged state at this interface */
         IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
         CHECKERR(ierr);
 
         /* Get the left and right eigenvectors */
         IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
         IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
         /* For each characteristic field */
         for (v = 0; v < nvars; v++) {
           /* calculate the characteristic flux components along this characteristic */
           double m1, p1, p2;
           m1 = p1 = p2 = 0;
           for (k = 0; k < nvars; k++) {
             m1 += L[v*nvars+k] * fC[qm1*nvars+k];
             p1 += L[v*nvars+k] * fC[qp1*nvars+k];
             p2 += L[v*nvars+k] * fC[qp2*nvars+k];
           }
 
           double slope_ratio = (p1 - m1) / ((p2 - p1) + 1e-40);
           double phi = muscl->LimiterFunction(slope_ratio);
           fchar[v] = p1 + 0.5 * phi * (p1-p2);
         }
 
         /* calculate the interface u from the characteristic u */
         IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
       }
     }
   }
 
   return(0);
 }

◆ Interp1PrimThirdOrderMUSCLChar()

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

3rd order MUSCL scheme with Koren's limiter (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the third-order MUSCL scheme

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 3rd order MUSCL scheme with Koren's limiter 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 3rd order MUSCL scheme with Koren's limiter as follows:

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

where

\begin{equation} \phi = \frac {3\left({\alpha}^k_j-{\alpha}^k_{j-1}\right)\left({\alpha}^k_{j-1}-{\alpha}^k_{j-2}\right) + \epsilon} {2\left[\left({\alpha}^k_j-{\alpha}^k_{j-1}\right)-\left({\alpha}^k_{j-1}-{\alpha}^k_{j-2}\right)\right]^2 + 3\left({\alpha}^k_j-{\alpha}^k_{j-1}\right)\left({\alpha}^k_{j-1}-{\alpha}^k_{j-2}\right) + \epsilon}, \end{equation}

\(\epsilon\) is a small constant (typically \(10^{-3}\)), 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:

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 (nvarsdim[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., Towards the Ultimate Conservative Difference Scheme. 2: Monotonicity and Conservation Combined in a Second-Order Scheme, J. of Comput. Phys., 14 (4), 1974, pp.361-370, http://dx.doi.org/10.1016/0021-9991(74)90019-9
Koren, B., A Robust Upwind Discretization Method for Advection, Diffusion and Source Terms, Centrum voor Wiskunde en Informatica, Amsterdam, 1993

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 91 of file Interp1PrimThirdOrderMUSCLChar.c.

 {
   HyPar             *solver = (HyPar*) s;
   MUSCLParameters   *muscl   = (MUSCLParameters*) solver->interp;
   int               i, k, v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   double one_third = 1.0/3.0;
   double one_sixth = 1.0/6.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
 
   if (upw > 0) {
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
 
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
         /* 1D indices of the stencil grid points */
         int qm1,qm2,qp1;
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
 
         int p; /* 1D index of the interface */
         _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
         /* find averaged state at this interface */
         IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
         CHECKERR(ierr);
 
         /* Get the left and right eigenvectors */
         IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
         IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
         /* For each characteristic field */
         for (v = 0; v < nvars; v++) {
           /* calculate the characteristic flux components along this characteristic */
           double m2, m1, p1;
           m2 = m1 = p1 = 0;
           for (k = 0; k < nvars; k++) {
             m2 += L[v*nvars+k] * fC[qm2*nvars+k];
             m1 += L[v*nvars+k] * fC[qm1*nvars+k];
             p1 += L[v*nvars+k] * fC[qp1*nvars+k];
           }
           double fdiff = p1 - m1;
           double bdiff = m1 - m2;
           double limit =  (3*fdiff*bdiff + muscl->eps) 
                         / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
           fchar[v] = m1 +  limit * (one_third*fdiff + one_sixth*bdiff);
         }
 
         /* calculate the interface u from the characteristic u */
         IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
       }
     }
   } else {
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
     for (i=0; i<N_outer; i++) {
       _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
       _ArrayCopy1D_(index_outer,indexC,ndims);
       _ArrayCopy1D_(index_outer,indexI,ndims);
 
       for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
         /* 1D indices of the stencil grid points */
         int qm1,qp1,qp2;
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
 
         int p; /* 1D index of the interface */
         _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
         /* find averaged state at this interface */
         IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics);
         CHECKERR(ierr);
 
         /* Get the left and right eigenvectors */
         IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
         IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
         /* For each characteristic field */
         for (v = 0; v < nvars; v++) {
           /* calculate the characteristic flux components along this characteristic */
           double m1, p1, p2;
           m1 = p1 = p2 = 0;
           for (k = 0; k < nvars; k++) {
             m1 += L[v*nvars+k] * fC[qm1*nvars+k];
             p1 += L[v*nvars+k] * fC[qp1*nvars+k];
             p2 += L[v*nvars+k] * fC[qp2*nvars+k];
           }
           double fdiff = p2 - p1;
           double bdiff = p1 - m1;
           double limit =  (3*fdiff*bdiff + muscl->eps) 
                         / (2*(fdiff-bdiff)*(fdiff-bdiff) + 3*fdiff*bdiff + muscl->eps);
           fchar[v] = p1 -  limit * (one_third*fdiff + one_sixth*bdiff);
         }
 
         /* calculate the interface u from the characteristic u */
         IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
       }
     }
   }
 
   return(0);
 }

◆ Interp1PrimFourthOrderCentralChar()

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

Characteristic-based interpolation of the first primitive at the cell interfaces using the fourth-order central scheme

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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 79 of file Interp1PrimFourthOrderCentralChar.c.

 {
   HyPar         *solver = (HyPar*) s;
   int           i, k, v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
 
   static const double c1 = 7.0 / 12.0;
   static const double c2 = -1.0 / 12.0;
 
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
 
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
       int pLL, pL, pR, pRR; /* 1D index of the left and right cells */
       indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pLL);
       indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pL );
       indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pR );
       indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,pRR);
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*pL],&u[nvars*pR],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       /* For each characteristic field */
       for (v = 0; v < nvars; v++) {
 
         /* calculate the characteristic flux components along this characteristic */
         double fcLL = 0, fcL = 0, fcR = 0, fcRR = 0;
         for (k = 0; k < nvars; k++) {
           fcLL += L[v*nvars+k] * fC[pLL*nvars+k];
           fcL  += L[v*nvars+k] * fC[pL *nvars+k];
           fcR  += L[v*nvars+k] * fC[pR *nvars+k];
           fcRR += L[v*nvars+k] * fC[pRR*nvars+k];
         }
 
         /* first order upwind approximation of the characteristic flux */
         fchar[v] = c2*fcLL + c1*fcL + c1*fcR + c2*fcRR;
 
       }
 
       /* calculate the interface u from the characteristic u */
       IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
 
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderUpwindChar()

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

5th order upwind reconstruction (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the fifth-order upwind scheme

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

\begin{align} \hat{\alpha}^k_{j+1/2} = \frac{1}{30} {\alpha}^k_{j-2} - \frac{13}{60}{\alpha}^k_{j-1} + \frac{47}{60}{\alpha}^k_j + \frac{27}{60}{\alpha}^k_{j+1} - \frac{1}{20}{\alpha}^k_{j+2}, \end{align}

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

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 (nvarsdim[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:

Jiang, G.-S., Shu, C.-W., Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126 (1), 1996, pp. 202-228, http://dx.doi.org/10.1006/jcph.1996.0130

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 84 of file Interp1PrimFifthOrderUpwindChar.c.

 {
   HyPar           *solver = (HyPar*) s;
   int             i, k, v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   static const double one_by_thirty        = 1.0/30.0,
                       thirteen_by_sixty    = 13.0/60.0,
                       fortyseven_by_sixty  = 47.0/60.0,
                       twentyseven_by_sixty = 27.0/60.0,
                       one_by_twenty        = 1.0/20.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
 
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
 
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
       /* 1D indices of the stencil grid points */
       int qm1,qm2,qm3,qp1,qp2;
       if (upw > 0) {
         indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       } else {
         indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       }
 
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       /* For each characteristic field */
       for (v = 0; v < nvars; v++) {
 
         /* calculate the characteristic flux components along this characteristic */
         double fm3, fm2, fm1, fp1, fp2;
         fm3 = fm2 = fm1 = fp1 = fp2 = 0;
         for (k = 0; k < nvars; k++) {
           fm3 += L[v*nvars+k] * fC[qm3*nvars+k];
           fm2 += L[v*nvars+k] * fC[qm2*nvars+k];
           fm1 += L[v*nvars+k] * fC[qm1*nvars+k];
           fp1 += L[v*nvars+k] * fC[qp1*nvars+k];
           fp2 += L[v*nvars+k] * fC[qp2*nvars+k];
         }
 
         /* fifth order upwind approximation of the characteristic flux */
         fchar[v] =    one_by_thirty         * fm3
                    -  thirteen_by_sixty     * fm2
                    +  fortyseven_by_sixty   * fm1
                    +  twentyseven_by_sixty  * fp1
                    -  one_by_twenty         * fp2;
       }
 
       /* calculate the interface u from the characteristic u */
       IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
 
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderCompactUpwindChar()

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

5th order compact upwind reconstruction (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the fifth-order compact upwind scheme

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

\begin{align} \frac{3}{10}\hat{\alpha}^k_{j-1/2} + \frac{6}{10}\hat{\alpha}^k_{j+1/2} + \frac{1}{10}\hat{\alpha}^k_{j+3/2} = \frac{}{30}{\alpha}^k_{j-1} + \frac{19}{30}{\alpha}^k_j + \frac{1}{3}{\alpha}^k_{j+1}, \end{align}

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 resulting block tridiagonal system is solved using blocktridiagLU() (see also TridiagLU, tridiagLU.h). 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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 90 of file Interp1PrimFifthOrderCompactUpwindChar.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   CompactScheme   *compact= (CompactScheme*)  solver->compact;
   TridiagLU       *lu     = (TridiagLU*)      solver->lusolver;
   int             sys,Nsys,d,v,k;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   static const double three_by_ten          = 3.0/10.0,
                       six_by_ten            = 6.0/10.0,
                       one_by_ten            = 1.0/10.0,
                       one_by_thirty         = 1.0/30.0,
                       nineteen_by_thirty    = 19.0/30.0,
                       one_third             = 1.0/3.0,
                       thirteen_by_sixty     = 13.0/60.0,
                       fortyseven_by_sixty   = 47.0/60.0,
                       twentyseven_by_sixty  = 27.0/60.0,
                       one_by_twenty         = 1.0/20.0;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
 
   /* calculate total number of block tridiagonal systems to solve */
   _ArrayProduct1D_(bounds_outer,ndims,Nsys);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars];
 
   /* Allocate arrays for tridiagonal system */
   double *A = compact->A;
   double *B = compact->B;
   double *C = compact->C;
   double *F = compact->R;
 
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,v,k,R,L,uavg,index_outer,indexC,indexI)
   for (sys=0; sys<Nsys; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2;
       if (upw > 0) {
         indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       } else {
         indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       }
 
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       for (v=0; v<nvars; v++)  {
 
         /* calculate the characteristic flux components along this characteristic */
         double fm3, fm2, fm1, fp1, fp2;
         fm3 = fm2 = fm1 = fp1 = fp2 = 0;
         for (k = 0; k < nvars; k++) {
           fm3 += L[v*nvars+k] * fC[qm3*nvars+k];
           fm2 += L[v*nvars+k] * fC[qm2*nvars+k];
           fm1 += L[v*nvars+k] * fC[qm1*nvars+k];
           fp1 += L[v*nvars+k] * fC[qp1*nvars+k];
           fp2 += L[v*nvars+k] * fC[qp2*nvars+k];
         }
 
         if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
             || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
           /* Use 5th order upwind at the physical boundaries */
           for (k=0; k<nvars; k++) {
             A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = 0.0;
             C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = 0.0;
             B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = L[v*nvars+k];
           }
           F[(Nsys*indexI[dir]+sys)*nvars+v] =   one_by_thirty         * fm3
                                               - thirteen_by_sixty     * fm2
                                               + fortyseven_by_sixty   * fm1
                                               + twentyseven_by_sixty  * fp1
                                               - one_by_twenty         * fp2;
         } else {
           /* 5th order compact upwind at the interior points */
           if (upw > 0) {
             for (k=0; k<nvars; k++) {
               A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = three_by_ten * L[v*nvars+k];
               B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = six_by_ten   * L[v*nvars+k];
               C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = one_by_ten   * L[v*nvars+k];
             }
           } else {
             for (k=0; k<nvars; k++) {
               C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = three_by_ten * L[v*nvars+k];
               B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = six_by_ten   * L[v*nvars+k];
               A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = one_by_ten   * L[v*nvars+k];
             }
           }
           F[(Nsys*indexI[dir]+sys)*nvars+v] =   one_by_thirty      * fm2
                                               + nineteen_by_thirty * fm1
                                               + one_third          * fp1;
         }
       }
     }
   }
 
 #ifdef serial
 
   /* Solve the tridiagonal system */
   IERR blocktridiagLU(A,B,C,F,dim[dir]+1,Nsys,nvars,lu,NULL); CHECKERR(ierr);
 
 #else
 
   /* Solve the tridiagonal system */
   /* all processes except the last will solve without the last interface to avoid overlap */
   if (mpi->ip[dir] != mpi->iproc[dir]-1)  { 
     IERR blocktridiagLU(A,B,C,F,dim[dir]  ,Nsys,nvars,lu,&mpi->comm[dir]); CHECKERR(ierr); 
   } else { 
     IERR blocktridiagLU(A,B,C,F,dim[dir]+1,Nsys,nvars,lu,&mpi->comm[dir]); CHECKERR(ierr);
   }
 
   /* Now get the solution to the last interface from the next proc */
   double *sendbuf = compact->sendbuf;
   double *recvbuf = compact->recvbuf;
   MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
   if (mpi->ip[dir]) for (d=0; d<Nsys*nvars; d++) sendbuf[d] = F[d];
   if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys*nvars,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
   if (mpi->ip[dir])                      MPI_Isend(sendbuf,Nsys*nvars,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
   MPI_Status status_arr[2];
   MPI_Waitall(2,&req[0],status_arr);
   if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys*nvars; d++) F[d+Nsys*nvars*dim[dir]] = recvbuf[d];
 
 #endif
 
   /* save the solution to fI */
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,v,k,R,L,uavg,index_outer,indexC,indexI)
   for (sys=0; sys<Nsys; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       _ArrayCopy1D_((F+sys*nvars+Nsys*nvars*indexI[dir]),(fI+nvars*p),nvars);
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderWENOChar()

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

5th order WENO reconstruction (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the fifth-order WENO scheme

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 WENO 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 WENO 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[ \hat{\alpha}^k_{j+1/2} = \frac{1}{3} {\alpha}^k_{j-2} - \frac{7}{6} {\alpha}^k_{j-1} + \frac{11}{6} {\alpha}^k_j \right]\\ + &\ \omega_2\ \times\ \left[ \hat{\alpha}^k_{j+1/2} = -\frac{1}{6} {\alpha}^k_{j-1} + \frac{5}{6} {\alpha}^k_j + \frac{1}{3} {\alpha}^k_{j+1} \right]\\ + &\ \omega_3\ \times\ \left[ \hat{\alpha}^k_{j+1/2} = \frac{1}{3} {\alpha}^k_j + \frac{5}{6} {\alpha}^k_{j+1} - \frac{1}{6} {\alpha}^k_{j+2} \right]\\ \Rightarrow &\ \hat{\alpha}^k_{j+1/2} = \frac{\omega_1}{3} {\alpha}^k_{j-2} - \frac{1}{6}(7\omega_1+\omega_2){\alpha}^k_{j-1} + \frac{1}{6}(11\omega_1+5\omega_2+2\omega_3){\alpha}^k_j + \frac{1}{6}(2\omega_2+5\omega_3){\alpha}^k_{j+1} - \frac{\omega_3}{6}{\alpha}^k_{j+2}, \end{align}

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 nonlinear weights \(\omega_k; k=1,2,3\) are the WENO weights computed in WENOFifthOrderCalculateWeightsChar(). 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:

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 (nvarsdim[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:

Jiang, G.-S., Shu, C.-W., Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126 (1), 1996, pp. 202-228, http://dx.doi.org/10.1006/jcph.1996.0130

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 86 of file Interp1PrimFifthOrderWENOChar.c.

 {
   HyPar           *solver = (HyPar*) s;
   WENOParameters  *weno   = (WENOParameters*) solver->interp;
   int             i, k, v;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   static const double one_sixth          = 1.0/6.0;
 
   double *ww1, *ww2, *ww3;
   ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
   int N_outer; _ArrayProduct1D_(bounds_outer,ndims,N_outer);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars], fchar[nvars];
 
 #pragma omp parallel for schedule(auto) default(shared) private(i,k,v,R,L,uavg,fchar,index_outer,indexC,indexI)
   for (i=0; i<N_outer; i++) {
     _ArrayIndexnD_(ndims,i,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
 
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
 
       /* 1D indices of the stencil grid points */
       int qm1,qm2,qm3,qp1,qp2;
       if (upw > 0) {
         indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       } else {
         indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       }
 
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       /* For each characteristic field */
       for (v = 0; v < nvars; v++) {
 
         /* calculate the characteristic flux components along this characteristic */
         double fm3, fm2, fm1, fp1, fp2;
         fm3 = fm2 = fm1 = fp1 = fp2 = 0;
         for (k = 0; k < nvars; k++) {
           fm3 += L[v*nvars+k] * fC[qm3*nvars+k];
           fm2 += L[v*nvars+k] * fC[qm2*nvars+k];
           fm1 += L[v*nvars+k] * fC[qm1*nvars+k];
           fp1 += L[v*nvars+k] * fC[qp1*nvars+k];
           fp2 += L[v*nvars+k] * fC[qp2*nvars+k];
         }
 
         /* Candidate stencils and their optimal weights*/
         double f1, f2, f3;
         f1 = (2*one_sixth)*fm3 - (7.0*one_sixth)*fm2 + (11.0*one_sixth)*fm1;
         f2 = (-one_sixth)*fm2 + (5.0*one_sixth)*fm1 + (2*one_sixth)*fp1;
         f3 = (2*one_sixth)*fm1 + (5*one_sixth)*fp1 - (one_sixth)*fp2;
 
         /* calculate WENO weights */
         double w1,w2,w3;
         w1 = *(ww1+p*nvars+v);
         w2 = *(ww2+p*nvars+v);
         w3 = *(ww3+p*nvars+v);
 
         /* fifth order WENO approximation of the characteristic flux */
         fchar[v] = w1*f1 + w2*f2 + w3*f3;
 
       }
 
       /* calculate the interface u from the characteristic u */
       IERR MatVecMult(nvars,(fI+nvars*p),R,fchar); CHECKERR(ierr);
 
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderCRWENOChar()

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

5th order CRWENO reconstruction (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the fifth-order CRWENO scheme

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{\alpha}^k_{j-1/2} + \frac{1}{3}\hat{\alpha}^k_{j+1/2} = \frac{1}{6} \left( f_{j-1} + 5f_j \right) \right]\\ + &\ \omega_2\ \times\ \left[ \frac{1}{3}\hat{\alpha}^k_{j-1/2}+\frac{2}{3}\hat{\alpha}^k_{j+1/2} = \frac{1}{6} \left( 5f_j + f_{j+1} \right) \right] \\ + &\ \omega_3\ \times\ \left[ \frac{2}{3}\hat{\alpha}^k_{j+1/2} + \frac{1}{3}\hat{\alpha}^k_{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{\alpha}^k_{j-1/2} + \left[\frac{1}{3}\omega_1+\frac{2}{3}(\omega_2+\omega_3)\right]\hat{\alpha}^k_{j+1/2} + \frac{1}{3}\omega_3\hat{\alpha}^k_{j+3/2} = \frac{\omega_1}{6}{\alpha}^k_{j-1} + \frac{5(\omega_1+\omega_2)+\omega_3}{6}{\alpha}^k_j + \frac{\omega_2+5\omega_3}{6}{\alpha}^k_{j+1}, \end{align}

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 nonlinear weights \(\omega_k; k=1,2,3\) are the WENO weights computed in WENOFifthOrderCalculateWeightsChar(). The resulting block tridiagonal system is solved using blocktridiagLU() (see also TridiagLU, tridiagLU.h). 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:

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 (nvarsdim[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

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 93 of file Interp1PrimFifthOrderCRWENOChar.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   CompactScheme   *compact= (CompactScheme*)  solver->compact;
   WENOParameters  *weno   = (WENOParameters*) solver->interp;
   TridiagLU       *lu     = (TridiagLU*)      solver->lusolver;
   int             sys,Nsys,d,v,k;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   static const double one_third          = 1.0/3.0;
   static const double one_sixth          = 1.0/6.0;
 
   double *ww1, *ww2, *ww3;
   ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
 
   /* calculate total number of block tridiagonal systems to solve */
   _ArrayProduct1D_(bounds_outer,ndims,Nsys);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars];
 
   /* Allocate arrays for tridiagonal system */
   double *A = compact->A;
   double *B = compact->B;
   double *C = compact->C;
   double *F = compact->R;
 
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,v,k,R,L,uavg,index_outer,indexC,indexI)
   for (sys=0; sys<Nsys; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2;
       if (upw > 0) {
         indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       } else {
         indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       }
 
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       for (v=0; v<nvars; v++)  {
 
         /* calculate the characteristic flux components along this characteristic */
         double fm3, fm2, fm1, fp1, fp2;
         fm3 = fm2 = fm1 = fp1 = fp2 = 0;
         for (k = 0; k < nvars; k++) {
           fm3 += L[v*nvars+k] * fC[qm3*nvars+k];
           fm2 += L[v*nvars+k] * fC[qm2*nvars+k];
           fm1 += L[v*nvars+k] * fC[qm1*nvars+k];
           fp1 += L[v*nvars+k] * fC[qp1*nvars+k];
           fp2 += L[v*nvars+k] * fC[qp2*nvars+k];
         }
 
         /* Candidate stencils and their optimal weights*/
         double f1, f2, f3;
         if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
             || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
           /* Use WENO5 at the physical boundaries */
           f1 = (2*one_sixth)*fm3 - (7.0*one_sixth)*fm2 + (11.0*one_sixth)*fm1;
           f2 = (-one_sixth)*fm2 + (5.0*one_sixth)*fm1 + (2*one_sixth)*fp1;
           f3 = (2*one_sixth)*fm1 + (5*one_sixth)*fp1 - (one_sixth)*fp2;
         } else {
           /* CRWENO5 at the interior points */
           f1 = (one_sixth) * (fm2 + 5*fm1);
           f2 = (one_sixth) * (5*fm1 + fp1);
           f3 = (one_sixth) * (fm1 + 5*fp1);
         }
 
         /* calculate WENO weights */
         double w1,w2,w3;
         w1 = *(ww1+p*nvars+v);
         w2 = *(ww2+p*nvars+v);
         w3 = *(ww3+p*nvars+v);
 
         if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
             || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
           for (k=0; k<nvars; k++) {
             A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = 0.0;
             C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = 0.0;
             B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = L[v*nvars+k];
           }
         } else {
           if (upw > 0) {
             for (k=0; k<nvars; k++) {
               A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = ((2*one_third)*w1 + (one_third)*w2)      * L[v*nvars+k];
               B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = ((one_third)*w1 + (2*one_third)*(w2+w3)) * L[v*nvars+k];
               C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = ((one_third)*w3)                         * L[v*nvars+k];
             }
           } else {
             for (k=0; k<nvars; k++) {
               C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = ((2*one_third)*w1 + (one_third)*w2)      * L[v*nvars+k];
               B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = ((one_third)*w1 + (2*one_third)*(w2+w3)) * L[v*nvars+k];
               A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = ((one_third)*w3)                         * L[v*nvars+k];
             }
           }
         }
         F[(Nsys*indexI[dir]+sys)*nvars+v] = w1*f1 + w2*f2 + w3*f3;
       }
     }
   }
 
 #ifdef serial
 
   /* Solve the tridiagonal system */
   IERR blocktridiagLU(A,B,C,F,dim[dir]+1,Nsys,nvars,lu,NULL); CHECKERR(ierr);
 
 #else
 
   /* Solve the tridiagonal system */
   /* all processes except the last will solve without the last interface to avoid overlap */
   if (mpi->ip[dir] != mpi->iproc[dir]-1)  { 
     IERR blocktridiagLU(A,B,C,F,dim[dir]  ,Nsys,nvars,lu,&mpi->comm[dir]); CHECKERR(ierr); 
   } else { 
     IERR blocktridiagLU(A,B,C,F,dim[dir]+1,Nsys,nvars,lu,&mpi->comm[dir]); CHECKERR(ierr);
   }
 
   /* Now get the solution to the last interface from the next proc */
   double *sendbuf = compact->sendbuf;
   double *recvbuf = compact->recvbuf;
   MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
   if (mpi->ip[dir]) for (d=0; d<Nsys*nvars; d++) sendbuf[d] = F[d];
   if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys*nvars,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
   if (mpi->ip[dir])                      MPI_Isend(sendbuf,Nsys*nvars,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
   MPI_Status status_arr[2];
   MPI_Waitall(2,&req[0],status_arr);
   if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys*nvars; d++) F[d+Nsys*nvars*dim[dir]] = recvbuf[d];
 
 #endif
 
   /* save the solution to fI */
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,v,k,R,L,uavg,index_outer,indexC,indexI)
   for (sys=0; sys<Nsys; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       _ArrayCopy1D_((F+sys*nvars+Nsys*nvars*indexI[dir]),(fI+nvars*p),nvars);
     }
   }
 
   return(0);
 }

◆ Interp1PrimFifthOrderHCWENOChar()

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

5th order hybrid compact-WENO reconstruction (characteristic-based) on a uniform grid

Characteristic-based interpolation of the first primitive at the cell interfaces using the fifth-order hybrid-compact WENO scheme

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 hybrid compact-WENO order scheme on a uniform grid. The reconstruction is carried out in terms of

\begin{equation} \alpha^k = {\bf l}_k \cdot {\bf f},\ k=1,\cdots,n, \end{equation}

which 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 block tridiagonal system is solved using blocktridiagLU() (see also TridiagLU, tridiagLU.h). 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}

See the references below for details of the hybrid compact-WENO scheme.

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 WENO weights are computed in WENOFifthOrderCalculateWeightsChar().
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:

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 (nvarsdim[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:

Pirozzoli, S., Conservative Hybrid Compact-WENO Schemes for Shock-Turbulence Interaction, J. Comput. Phys., 178 (1), 2002, pp. 81-117, http://dx.doi.org/10.1006/jcph.2002.7021
Ren, Y.-X., Liu, M., Zhang, H., A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws, J. Comput. Phys., 192 (2), 2003, pp. 365-386,

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
u	Array of cell-centered values of the solution \({\bf u}\)
x	Grid coordinates
upw	Upwind direction (left or right biased)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables
uflag	Flag to indicate if \(f(u) \equiv u\), i.e, if the solution is being reconstructed

Definition at line 79 of file Interp1PrimFifthOrderHCWENOChar.c.

 {
   HyPar           *solver = (HyPar*)          s;
   MPIVariables    *mpi    = (MPIVariables*)   m;
   CompactScheme   *compact= (CompactScheme*)  solver->compact;
   WENOParameters  *weno   = (WENOParameters*) solver->interp;
   TridiagLU       *lu     = (TridiagLU*)      solver->lusolver;
   int             sys,Nsys,d,v,k;
   _DECLARE_IERR_;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* define some constants */
   static const double one_half           = 1.0/2.0;
   static const double one_third          = 1.0/3.0;
   static const double one_sixth          = 1.0/6.0;
 
   double *ww1, *ww2, *ww3;
   ww1 = weno->w1 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww2 = weno->w2 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
   ww3 = weno->w3 + (upw < 0 ? 2*weno->size : 0) + (uflag ? weno->size : 0) + weno->offset[dir];
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
 
   /* calculate total number of block tridiagonal systems to solve */
   _ArrayProduct1D_(bounds_outer,ndims,Nsys);
 
   /* allocate arrays for the averaged state, eigenvectors and characteristic interpolated f */
   double R[nvars*nvars], L[nvars*nvars], uavg[nvars];
 
   /* Allocate arrays for tridiagonal system */
   double *A = compact->A;
   double *B = compact->B;
   double *C = compact->C;
   double *F = compact->R;
 
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,v,k,R,L,uavg,index_outer,indexC,indexI)
   for (sys=0; sys<Nsys; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int qm1,qm2,qm3,qp1,qp2;
       if (upw > 0) {
         indexC[dir] = indexI[dir]-3; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       } else {
         indexC[dir] = indexI[dir]+2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm3);
         indexC[dir] = indexI[dir]+1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm2);
         indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qm1);
         indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp1);
         indexC[dir] = indexI[dir]-2; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qp2);
       }
 
       int p; /* 1D index of the interface */
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
 
       /* find averaged state at this interface */
       IERR solver->AveragingFunction(uavg,&u[nvars*qm1],&u[nvars*qp1],solver->physics); CHECKERR(ierr);
 
       /* Get the left and right eigenvectors */
       IERR solver->GetLeftEigenvectors  (uavg,L,solver->physics,dir); CHECKERR(ierr);
       IERR solver->GetRightEigenvectors (uavg,R,solver->physics,dir); CHECKERR(ierr);
 
       for (v=0; v<nvars; v++)  {
 
         /* calculate the characteristic flux components along this characteristic */
         double fm3, fm2, fm1, fp1, fp2;
         fm3 = fm2 = fm1 = fp1 = fp2 = 0;
         for (k = 0; k < nvars; k++) {
           fm3 += L[v*nvars+k] * fC[qm3*nvars+k];
           fm2 += L[v*nvars+k] * fC[qm2*nvars+k];
           fm1 += L[v*nvars+k] * fC[qm1*nvars+k];
           fp1 += L[v*nvars+k] * fC[qp1*nvars+k];
           fp2 += L[v*nvars+k] * fC[qp2*nvars+k];
         }
 
         /* Candidate stencils and their optimal weights*/
         double f1, f2, f3;
         if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
             || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
           /* Use WENO5 at the physical boundaries */
           f1 = (2*one_sixth)*fm3 - (7.0*one_sixth)*fm2 + (11.0*one_sixth)*fm1;
           f2 = (-one_sixth)*fm2 + (5.0*one_sixth)*fm1 + (2*one_sixth)*fp1;
           f3 = (2*one_sixth)*fm1 + (5*one_sixth)*fp1 - (one_sixth)*fp2;
         } else {
           /* HCWENO5 at the interior points */
           f1 = (one_sixth) * (fm2 + 5*fm1);
           f2 = (one_sixth) * (5*fm1 + fp1);
           f3 = (one_sixth) * (fm1 + 5*fp1);
         }
 
         /* calculate WENO weights */
         double w1,w2,w3;
         w1 = *(ww1+p*nvars+v);
         w2 = *(ww2+p*nvars+v);
         w3 = *(ww3+p*nvars+v);
 
         /* calculate the hybridization parameter */
         double sigma;
         if (   ((mpi->ip[dir] == 0                ) && (indexI[dir] == 0       ))
             || ((mpi->ip[dir] == mpi->iproc[dir]-1) && (indexI[dir] == dim[dir])) ) {
           /* Standard WENO at physical boundaries */
           sigma = 0.0;
         } else {
           double cuckoo, df_jm12, df_jp12, df_jp32, r_j, r_jp1, r_int;
           cuckoo = (0.9*weno->rc / (1.0-0.9*weno->rc)) * weno->xi * weno->xi;
           df_jm12 = fm1 - fm2;
           df_jp12 = fp1 - fm1;
           df_jp32 = fp2 - fp1;
           r_j   = (absolute(2*df_jp12*df_jm12)+cuckoo)/(df_jp12*df_jp12+df_jm12*df_jm12+cuckoo);
           r_jp1 = (absolute(2*df_jp32*df_jp12)+cuckoo)/(df_jp32*df_jp32+df_jp12*df_jp12+cuckoo);
           r_int = min(r_j, r_jp1);
           sigma = min((r_int/weno->rc), 1.0); 
         }
 
         if (upw > 0) {
           for (k=0; k<nvars; k++) {
             A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = (one_half*sigma)  * L[v*nvars+k];
             B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = (1.0)             * L[v*nvars+k];
             C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = (one_sixth*sigma) * L[v*nvars+k];
           }
         } else {
           for (k=0; k<nvars; k++) {
             C[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = (one_half*sigma)  * L[v*nvars+k];
             B[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = (1.0)             * L[v*nvars+k];
             A[(Nsys*indexI[dir]+sys)*nvars*nvars+v*nvars+k] = (one_sixth*sigma) * L[v*nvars+k];
           }
         }
         double fWENO, fCompact;
         fWENO    = w1*f1 + w2*f2 + w3*f3;
         fCompact = one_sixth * (one_third*fm2 + 19.0*one_third*fm1 + 10.0*one_third*fp1);
         F[(Nsys*indexI[dir]+sys)*nvars+v] = sigma*fCompact + (1.0-sigma)*fWENO;
       }
     }
   }
 
 #ifdef serial
 
   /* Solve the tridiagonal system */
   IERR blocktridiagLU(A,B,C,F,dim[dir]+1,Nsys,nvars,lu,NULL); CHECKERR(ierr);
 
 #else
 
   /* Solve the tridiagonal system */
   /* all processes except the last will solve without the last interface to avoid overlap */
   if (mpi->ip[dir] != mpi->iproc[dir]-1)  { 
     IERR blocktridiagLU(A,B,C,F,dim[dir]  ,Nsys,nvars,lu,&mpi->comm[dir]); CHECKERR(ierr); 
   } else { 
     IERR blocktridiagLU(A,B,C,F,dim[dir]+1,Nsys,nvars,lu,&mpi->comm[dir]); CHECKERR(ierr);
   }
 
   /* Now get the solution to the last interface from the next proc */
   double *sendbuf = compact->sendbuf;
   double *recvbuf = compact->recvbuf;
   MPI_Request req[2] = {MPI_REQUEST_NULL,MPI_REQUEST_NULL};
   if (mpi->ip[dir]) for (d=0; d<Nsys*nvars; d++) sendbuf[d] = F[d];
   if (mpi->ip[dir] != mpi->iproc[dir]-1) MPI_Irecv(recvbuf,Nsys*nvars,MPI_DOUBLE,mpi->ip[dir]+1,214,mpi->comm[dir],&req[0]);
   if (mpi->ip[dir])                      MPI_Isend(sendbuf,Nsys*nvars,MPI_DOUBLE,mpi->ip[dir]-1,214,mpi->comm[dir],&req[1]);
   MPI_Status status_arr[2];
   MPI_Waitall(2,&req[0],status_arr);
   if (mpi->ip[dir] != mpi->iproc[dir]-1) for (d=0; d<Nsys*nvars; d++) F[d+Nsys*nvars*dim[dir]] = recvbuf[d];
 
 #endif
 
   /* save the solution to fI */
 #pragma omp parallel for schedule(auto) default(shared) private(sys,d,v,k,R,L,uavg,index_outer,indexC,indexI)
   for (sys=0; sys<Nsys; sys++) {
     _ArrayIndexnD_(ndims,sys,bounds_outer,index_outer,0);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int p; _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       _ArrayCopy1D_((F+sys*nvars+Nsys*nvars*indexI[dir]),(fI+nvars*p),nvars);
     }
   }
 
   return(0);
 }

◆ Interp2PrimSecondOrder()

int Interp2PrimSecondOrder	(	double *	fI,
		double *	fC,
		int	dir,
		void *	s,
		void *	m
	)

2nd order component-wise interpolation of the 2nd primitive on a uniform grid

Interpolation of the second primitive at the cell interfaces using the second-order central scheme

Computes the interpolated values of the second primitive of a function \({\bf f}\left({\bf u}\right)\) at the interfaces from the cell-centered values of the function using the second order central method on a uniform grid. The second 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^2} \int_{x-\Delta x/2}^{x+\Delta x/2} \left( \int_{\eta-\Delta x/2}^{\eta+\Delta x/2} {\bf h}\left({\bf u}\left(\zeta\right)\right) d\zeta \right) d\eta, \end{equation}

where \(x\) is the spatial coordinate along the dimension of the interpolation. This function the 2nd order central numerical approximation \(\hat{\bf f}_{j+1/2} \approx {\bf h}_{j+1/2}\) as \(\hat{\bf f}_{j+1/2} \approx {\bf h}_{j+1/2}\) as:

\begin{equation} \hat{\bf f}_{j+1/2} = \frac{1}{2}\left( {\bf f}_{j+1} - {\bf f}_j \right). \end{equation}

Implementation Notes:

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:

Reference:

Liu, Y., Shu, C.-W., Zhang, M., High Order Finite Difference WENO Schemes for Nonlinear Degenerate Parabolic Equations, SIAM J. Sci. Comput., 33 (2), 2011, pp. 939-965, http://dx.doi.org/10.1137/100791002.

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.
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.

Parameters

fI	Array of interpolated function values at the interfaces
fC	Array of cell-centered values of the function \({\bf f}\left({\bf u}\right)\)
dir	Spatial dimension along which to interpolation
s	Object of type HyPar containing solver-related variables
m	Object of type MPIVariables containing MPI-related variables

Definition at line 56 of file Interp2PrimSecondOrder.c.

 {
   HyPar         *solver = (HyPar*)        s;
 
   int ghosts = solver->ghosts;
   int ndims  = solver->ndims;
   int nvars  = solver->nvars;
   int *dim   = solver->dim_local;
 
   /* create index and bounds for the outer loop, i.e., to loop over all 1D lines along
      dimension "dir"                                                                    */
   int indexC[ndims], indexI[ndims], index_outer[ndims], bounds_outer[ndims], bounds_inter[ndims];
   _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] =  1;
   _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
 
   int done = 0; _ArraySetValue_(index_outer,ndims,0);
   while (!done) {
     _ArrayCopy1D_(index_outer,indexC,ndims);
     _ArrayCopy1D_(index_outer,indexI,ndims);
     for (indexI[dir] = 0; indexI[dir] < dim[dir]+1; indexI[dir]++) {
       int p, qL, qR;
       _ArrayIndex1D_(ndims,bounds_inter,indexI,0,p);
       indexC[dir] = indexI[dir]-1; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qL);
       indexC[dir] = indexI[dir]  ; _ArrayIndex1D_(ndims,dim,indexC,ghosts,qR);
       int v; for (v=0; v<nvars; v++) fI[p*nvars+v] = (fC[qR*nvars+v] - fC[qL*nvars+v]);
     }
     _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
   }
   
   return(0);
 }

◆ InterpSetLimiterVar()

int InterpSetLimiterVar	(	double *	,
		double *	,
		double *	,
		int	,
		void *	,
		void *
	)

Function to calculate and save the nonlinear interpolation coefficients for a non-linear (solution-dependent) interpolation scheme (eg. _FIFTH_ORDER_WENO_, _FIFTH_ORDER_CRWENO_, _FIFTH_ORDER_HCWENO_, etc).

◆ MUSCLInitialize()

int MUSCLInitialize	(	void *	s,
		void *	m
	)

Initialize the 2nd or 3rd order MUSCL scheme.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables

Definition at line 17 of file MUSCLInitialize.c.

 {
   HyPar           *solver = (HyPar*) s;
   MPIVariables    *mpi    = (MPIVariables*) m;
   MUSCLParameters *muscl  = (MUSCLParameters*) solver->interp;
 
   /* default values */
   muscl->eps = 1e-3;
   strcpy(muscl->limiter_type, _LIM_GM_);
 
   if (!mpi->rank) {
     FILE *in;
     int ferr;
     in = fopen("muscl.inp","r");
     if (!in) printf("Warning: File muscl.inp not found. Using default parameters for muscl2/muscl3 scheme.\n");
     else {
       printf("Reading MUSCL parameters from muscl.inp.\n");
       char word[_MAX_STRING_SIZE_];
       ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
       if (!strcmp(word, "begin")){
           while (strcmp(word, "end")){
               ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
           if      (!strcmp(word,"epsilon")) { ferr = fscanf(in,"%lf",&muscl->eps);         if (ferr != 1) return(1); }
           else if (!strcmp(word,"limiter")) { ferr = fscanf(in,"%s" ,muscl->limiter_type); if (ferr != 1) return(1); }
           else if (strcmp(word,"end")) {
             char useless[_MAX_STRING_SIZE_];
             ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
             printf("Warning: keyword %s in file \"muscl.inp\" with value %s not ",word,useless);
             printf("recognized or extraneous. Ignoring.\n");
           }
         }
         } else {
           fprintf(stderr,"Error: Illegal format in file \"muscl.inp\".\n");
         return(1);
         }
       fclose(in);
     }
   }
 
   MPIBroadcast_double(&muscl->eps, 1, 0, &mpi->world);
   MPIBroadcast_character(muscl->limiter_type, _MAX_STRING_SIZE_, 0, &mpi->world);
 
   if (!strcmp(muscl->limiter_type, _LIM_GM_)) {
     muscl->LimiterFunction = LimiterGeneralizedMinMod;
   } else if (!strcmp(muscl->limiter_type, _LIM_MM_)) {
     muscl->LimiterFunction = LimiterMinMod;
   } else if (!strcmp(muscl->limiter_type, _LIM_VANLEER_)) {
     muscl->LimiterFunction = LimiterVanLeer;
   } else if (!strcmp(muscl->limiter_type, _LIM_SUPERBEE_)) {
     muscl->LimiterFunction = LimiterSuperBee;
   } else {
     if (!mpi->rank) {
       fprintf(stderr, "Warning: %s is an invalid limiter type. Using default (Generalized MinMod).\n", 
               muscl->limiter_type);
     }
     muscl->LimiterFunction = LimiterGeneralizedMinMod;
   }
   
   return(0);
 }

◆ WENOInitialize()

int WENOInitialize	(	void *	s,
		void *	m,
		char *	scheme,
		char *	type
	)

Initialize the structure containing variables and parameters for WENO-type schemes

This function initializes the WENO-type methods.

Sets the parameters to default values.
Reads in the parameters from optional input file "weno.inp", if available.
Allocates memory for and initializes the nonlinear weights used by WENO-type schemes.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
scheme	Name of scheme
type	Type of interpolation

Definition at line 34 of file WENOInitialize.c.

 {
   HyPar           *solver = (HyPar*) s;
   MPIVariables    *mpi    = (MPIVariables*) m;
   WENOParameters  *weno   = (WENOParameters*) solver->interp;
 
   static int count = 0;
 
   int nvars = solver->nvars;
   int ndims = solver->ndims;
 
   /* default parameters */
   weno->mapped      = 0;
   weno->borges      = 0;
   weno->yc          = 0;
   weno->no_limiting = 0;
   weno->eps         = 1e-6;
   weno->p           = 2.0;
 
   weno->rc          = 0.3;
   weno->xi          = 0.001;
   weno->tol         = 1e-16;
 
   if (!mpi->rank) {
     FILE *in;
     int ferr;
     in = fopen("weno.inp","r");
     if (!in) printf("Warning: File weno.inp not found. Using default parameters for WENO5/CRWENO5/HCWENO5 scheme.\n");
     else {
       if (!count) printf("Reading WENO parameters from weno.inp.\n");
       char word[_MAX_STRING_SIZE_];
       ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
       if (!strcmp(word, "begin")){
           while (strcmp(word, "end")){
               ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
           if      (!strcmp(word,"mapped"     )) { ferr = fscanf(in,"%d" ,&weno->mapped     ); if (ferr != 1) return(1); }
           else if (!strcmp(word,"borges"     )) { ferr = fscanf(in,"%d" ,&weno->borges     ); if (ferr != 1) return(1); }
           else if (!strcmp(word,"yc"         )) { ferr = fscanf(in,"%d" ,&weno->yc         ); if (ferr != 1) return(1); }
           else if (!strcmp(word,"no_limiting")) { ferr = fscanf(in,"%d" ,&weno->no_limiting); if (ferr != 1) return(1); }
           else if (!strcmp(word,"epsilon"    )) { ferr = fscanf(in,"%lf",&weno->eps        ); if (ferr != 1) return(1); }
           else if (!strcmp(word,"p"          )) { ferr = fscanf(in,"%lf",&weno->p          ); if (ferr != 1) return(1); }
           else if (!strcmp(word,"rc"         )) { ferr = fscanf(in,"%lf",&weno->rc         ); if (ferr != 1) return(1); }
           else if (!strcmp(word,"xi"         )) { ferr = fscanf(in,"%lf",&weno->xi         ); if (ferr != 1) return(1); }
           else if (!strcmp(word,"tol"        )) { ferr = fscanf(in,"%lf",&weno->tol        ); if (ferr != 1) return(1); }
           else if (strcmp(word,"end")) {
             char useless[_MAX_STRING_SIZE_];
             ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
             printf("Warning: keyword %s in file \"weno.inp\" with value %s not ",word,useless);
             printf("recognized or extraneous. Ignoring.\n");
           }
         }
         } else {
           fprintf(stderr,"Error: Illegal format in file \"weno.inp\".\n");
         return(1);
         }
       fclose(in);
     }
   }
 
   int     integer_data[4];
   double  real_data[5];
   if (!mpi->rank) {
     integer_data[0] = weno->mapped;
     integer_data[1] = weno->borges;
     integer_data[2] = weno->yc;
     integer_data[3] = weno->no_limiting;
     real_data[0]    = weno->eps;
     real_data[1]    = weno->p;
     real_data[2]    = weno->rc;
     real_data[3]    = weno->xi;
     real_data[4]    = weno->tol;
   }
   MPIBroadcast_integer(integer_data,4,0,&mpi->world);
   MPIBroadcast_double (real_data   ,5,0,&mpi->world);
 
   weno->mapped      = integer_data[0];
   weno->borges      = integer_data[1];
   weno->yc          = integer_data[2];
   weno->no_limiting = integer_data[3];
   weno->eps         = real_data   [0];
   weno->p           = real_data   [1];
   weno->rc          = real_data   [2];
   weno->xi          = real_data   [3];
   weno->tol         = real_data   [4];
 
   /* WENO weight calculation is hard-coded for p=2, so return error if p != 2 in
    * user input file, so that there's no confusion */
   if (weno->p != 2.0) {
     if (!mpi->rank && !count) printf("Warning from WENOInitialize(): \"p\" parameter is 2.0. Any other value will be ignored!\n");
   }
 
   weno->offset = NULL;
   weno->w1 = NULL;
   weno->w2 = NULL;
   weno->w3 = NULL;
 
   weno->offset = (int*) calloc (ndims,sizeof(int));
   int dir,d;
   for (dir=0; dir<ndims; dir++) {
     weno->offset[dir] = 0;
     for (d=0; d<dir; d++) {
       int size = nvars, i;
       for (i=0; i<ndims; i++)
         size *= ( i==d ? solver->dim_local[i]+1 : solver->dim_local[i] );
       weno->offset[dir] += size;
     }
   }
 
   int total_size = 0;
   for (d=0; d<ndims; d++) {
     int size = nvars, i;
     for (i=0; i<ndims; i++)
       size *= ( i==d ? solver->dim_local[i]+1 : solver->dim_local[i] );
     total_size += size;
   }
   weno->size = total_size;
 
   if ((!strcmp(type,_CHARACTERISTIC_)) && (nvars > 1))
     solver->SetInterpLimiterVar = WENOFifthOrderCalculateWeightsChar;
   else {
 #if defined(HAVE_CUDA)
     if (solver->use_gpu) {
       solver->SetInterpLimiterVar = gpuWENOFifthOrderCalculateWeights;
     } else {
 #endif
       solver->SetInterpLimiterVar = WENOFifthOrderCalculateWeights;
 #if defined(HAVE_CUDA)
     }
 #endif
   }
 
   /* initialize WENO weights to their optimal values */
   double* tmp_w1 = (double*) calloc (4*total_size,sizeof(double));
   double* tmp_w2 = (double*) calloc (4*total_size,sizeof(double));
   double* tmp_w3 = (double*) calloc (4*total_size,sizeof(double));
   for (d=0; d<ndims; d++) WENOFifthOrderInitializeWeights(  tmp_w1, 
                                                             tmp_w2, 
                                                             tmp_w3, 
                                                             weno->offset, 
                                                             d,
                                                             solver,
                                                             mpi);
   count++;
 
 #if defined(HAVE_CUDA)
   if (solver->use_gpu) {
     //gpuMalloc((void**)&weno->gpu_offset, solver->ndims*sizeof(int));
     gpuMalloc((void**)&weno->w1, 4*total_size*sizeof(double));
     gpuMalloc((void**)&weno->w2, 4*total_size*sizeof(double));
     gpuMalloc((void**)&weno->w3, 4*total_size*sizeof(double));
 
     //gpuMemcpy(weno->gpu_offset, weno->offset, solver->ndims*sizeof(int), gpuMemcpyHostToDevice);
     gpuMemcpy(weno->w1, tmp_w1, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
     gpuMemcpy(weno->w2, tmp_w2, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
     gpuMemcpy(weno->w3, tmp_w3, 4*total_size*sizeof(double), gpuMemcpyHostToDevice);
   } else {
 #endif
     weno->w1 = (double*) calloc (4*total_size,sizeof(double));
     weno->w2 = (double*) calloc (4*total_size,sizeof(double));
     weno->w3 = (double*) calloc (4*total_size,sizeof(double));
 
     _ArrayCopy1D_(tmp_w1, weno->w1, 4*total_size);
     _ArrayCopy1D_(tmp_w2, weno->w2, 4*total_size);
     _ArrayCopy1D_(tmp_w3, weno->w3, 4*total_size);
 #if defined(HAVE_CUDA)
   }
 #endif
 
   free(tmp_w1);
   free(tmp_w2);
   free(tmp_w3);
 
   return 0;
 }

◆ WENOCleanup()

int WENOCleanup	(	void *	s,
		int	flag_gpu
	)

Clean up the structure containing variables and parameters for WENO-type schemes

Cleans up all allocations related to the WENO-type methods.

Parameters

s	WENO object of type WENOParameters
flag_gpu	flag to indicate if on GPU

Definition at line 17 of file WENOCleanup.c.

 {
   WENOParameters  *weno   = (WENOParameters*) s;
 
   if (weno->offset) free(weno->offset);
 #if defined(HAVE_CUDA)
   if (flag_gpu) {
     if (weno->w1) gpuFree(weno->w1);
     if (weno->w2) gpuFree(weno->w2);
     if (weno->w3) gpuFree(weno->w3);
   } else {
 #endif
     if (weno->w1) free(weno->w1);
     if (weno->w2) free(weno->w2);
     if (weno->w3) free(weno->w3);
 #if defined(HAVE_CUDA)
   }
 #endif
 
   return(0);
 }

◆ CompactSchemeInitialize()

int CompactSchemeInitialize	(	void *	s,
		void *	m,
		char *	type
	)

Initialize the structure containing variables and parameters for compact schemes

This function initializes the compact finite-difference methods: allocates the arrays to store the tridiagonal system.

Parameters

s	Solver object of type HyPar
m	MPI object of type MPIVariables
type	Type of interpolation

Definition at line 17 of file CompactSchemeInitialize.c.

 {
   HyPar         *solver = (HyPar*) s;
   MPIVariables  *mpi    = (MPIVariables*) m;
   CompactScheme *compact   = (CompactScheme*) solver->compact;
 
   int nvars = solver->nvars;
   int ndims = solver->ndims;
 
   int size = 1, d;
   for (d=0; d<solver->ndims; d++) size *= (solver->dim_local[d]+1);
   size *= solver->nvars;
   if (!strcmp(solver->interp_type,_CHARACTERISTIC_)) size *= solver->nvars;
 
   compact->A = (double*) calloc (size, sizeof(double));
   compact->B = (double*) calloc (size, sizeof(double));
   compact->C = (double*) calloc (size, sizeof(double));
   compact->R = (double*) calloc (size, sizeof(double));
 
   compact->sendbuf = (double*) calloc (size, sizeof(double));
   compact->recvbuf = (double*) calloc (size, sizeof(double));
 
   return(0);
 }

◆ CompactSchemeCleanup()

int CompactSchemeCleanup ( void * s )

Clean up the structure containing variables and parameters for compact schemes

Cleans up all allocations related to the compact finite difference methods.

Parameters

s	CompactScheme object of type CompactScheme

Definition at line 14 of file CompactSchemeCleanup.c.

 {
   CompactScheme *compact   = (CompactScheme*) s;
 
   if (compact->A) free(compact->A);
   if (compact->B) free(compact->B);
   if (compact->C) free(compact->C);
   if (compact->R) free(compact->R);
 
   if (compact->sendbuf) free(compact->sendbuf);
   if (compact->recvbuf) free(compact->recvbuf);
 
   return(0);
 }

Data Structures

Macros

Functions

Detailed Description

Data Structure Documentation

◆ WENOParameters

◆ CompactScheme

Macro Definition Documentation

◆ _FIRST_ORDER_UPWIND_

◆ _SECOND_ORDER_CENTRAL_

◆ _SECOND_ORDER_MUSCL_

◆ _THIRD_ORDER_MUSCL_

◆ _FOURTH_ORDER_CENTRAL_

◆ _FIFTH_ORDER_UPWIND_

◆ _FIFTH_ORDER_COMPACT_UPWIND_

◆ _FIFTH_ORDER_WENO_

◆ _FIFTH_ORDER_CRWENO_

◆ _FIFTH_ORDER_HCWENO_

◆ _CHARACTERISTIC_

◆ _COMPONENTS_

◆ _WENO_OPTIMAL_WEIGHT_1_

◆ _WENO_OPTIMAL_WEIGHT_2_

◆ _WENO_OPTIMAL_WEIGHT_3_

◆ _CRWENO_OPTIMAL_WEIGHT_1_

◆ _CRWENO_OPTIMAL_WEIGHT_2_

◆ _CRWENO_OPTIMAL_WEIGHT_3_

◆ _WENOWeights_v_JS_

◆ _WENOWeights_v_M_

◆ _WENOWeights_v_M_Scalar_

◆ _WENOWeights_v_Z_

◆ _WENOWeights_v_YC_

◆ _WENOWeights_v_YC_Scalar_

Function Documentation

◆ Interp1PrimFirstOrderUpwind()

◆ Interp1PrimSecondOrderCentral()

◆ Interp1PrimSecondOrderMUSCL()

◆ Interp1PrimThirdOrderMUSCL()

◆ Interp1PrimFourthOrderCentral()

◆ Interp1PrimFifthOrderUpwind()

◆ Interp1PrimFifthOrderCompactUpwind()

◆ Interp1PrimFifthOrderWENO()

◆ Interp1PrimFifthOrderCRWENO()

◆ Interp1PrimFifthOrderHCWENO()

◆ gpuInterp1PrimFifthOrderWENO()

◆ Interp1PrimFirstOrderUpwindChar()

◆ Interp1PrimSecondOrderCentralChar()

◆ Interp1PrimSecondOrderMUSCLChar()

◆ Interp1PrimThirdOrderMUSCLChar()

◆ Interp1PrimFourthOrderCentralChar()

◆ Interp1PrimFifthOrderUpwindChar()

◆ Interp1PrimFifthOrderCompactUpwindChar()

◆ Interp1PrimFifthOrderWENOChar()

◆ Interp1PrimFifthOrderCRWENOChar()

◆ Interp1PrimFifthOrderHCWENOChar()

◆ Interp2PrimSecondOrder()

◆ InterpSetLimiterVar()

◆ MUSCLInitialize()

◆ WENOInitialize()

◆ WENOCleanup()

◆ CompactSchemeInitialize()

◆ CompactSchemeCleanup()