interpolation.h

Go to the documentation of this file.

#ifndef _INTERP_H_
#define _INTERP_H_

#include "basic.h"

#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"

/* interpolation type definitions */
#define _CHARACTERISTIC_        "characteristic" 
#define _COMPONENTS_            "components"     
/*
  One-dimensional Interpolation Functions:-
    Functions to interpolate the primitive of a function at interfaces from its
    cell-centered values.

  Arguments:-

    fI        double*   array of size (N+1) in the interpolation direction; will contain
                        the interpolated interface function primitive
                        (**needs to be allocated by the calling function)
                        (**does not have ghost points))

    fC        double*   array of size (N+2*ghosts) in the interpolation direction of the
                        cell centered function values
                        (**does have ghost points))

    u         double*   used only by the characteristic-based interpolation schemes
                        array of size (N+2*ghosts) in the interpolation direction of the
                        cell centered conserved variable (needed to calculate the eigen-
                        decomposition at the interfaces)

    x         double*   grid point locations along the 1D line on which the interpolation
                        is being carried out
                        array of size (N+2*ghosts)
                        Used only by non-uniform-grid interpolation schemes

    upw       int       upwind direction for non-central schemes
                        (1 -> left biased, -1 -> right biased)

    dir       int       direction/dimension along which to carry out interpolation for
                        multi-dimensional domains
                        (eg: 0 for 1D; 0 or 1 for 2D; 0,1 or 2 for 3D)

    s         void*     object of type depending on the main solver, should have at least
                        the following data:
                        + ghosts    : number of ghosts points in the domain
                        + ndims     : number of dimensions
                        + nvars     : number of variables per grid point
                                      (i.e. size of vector variable being interpolated)
                        + dim_local : local (of this process) dimensions of the domain
                                      (integer array of size ndims, with the number of
                                      points in each dimension as elements)
    m         void*     object containing MPI domain decomposition related variables
                        (this information is used only by compact interpolation schemes)

    uflag     int       flag to indicate where the flux function or the solution function
                        is being interpolated
                        (1 -> u; 0 -> f(u) )

    Notes:
      + arrangement of points along the direction of interpolation is

              |<--ghosts-->|0 1 2 ..(interior).. N-1|<--ghosts-->|

      + number of interfaces are one more than number of interior points
      + interface i lies between grid points i-1+ghosts and i+ghosts

            0   1   2   .....        ....          N-3 N-2 N-1
          | x | x | x | x | x | x | x | x | x | x | x | x | x |
          0   1   2   ...             ...    ...     N-2 N-1  N

    Returns 0 on normal execution, non-zero on error.
*/

#ifdef __cplusplus
extern "C" {
#endif

/* functions to interpolate the first primitive in a component-wise way
   (for conservative discretization of the 1st derivative) on a uniform grid */
int Interp1PrimFirstOrderUpwind           (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimSecondOrderCentral         (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimSecondOrderMUSCL           (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimThirdOrderMUSCL            (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFourthOrderCentral         (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderUpwind           (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderCompactUpwind    (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderWENO             (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderCRWENO           (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderHCWENO           (double*,double*,double*,double*,int,int,void*,void*,int);

#if defined(HAVE_CUDA)

int gpuInterp1PrimFifthOrderWENO (double*,double*,double*,double*,int,int,void*,void*,int);
#endif

/* functions to interpolate the first primitive in a characteristic-based way
   (for conservative discretization of the 1st derivative) on a uniform grid */
int Interp1PrimFirstOrderUpwindChar       (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimSecondOrderCentralChar     (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimSecondOrderMUSCLChar       (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimThirdOrderMUSCLChar        (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFourthOrderCentralChar     (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderUpwindChar       (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderCompactUpwindChar(double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderWENOChar         (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderCRWENOChar       (double*,double*,double*,double*,int,int,void*,void*,int);
int Interp1PrimFifthOrderHCWENOChar       (double*,double*,double*,double*,int,int,void*,void*,int);

/* functions to interpolate the second primitive
   (for conservative discretization of the 2nd derivative) */
int Interp2PrimSecondOrder  (double*,double*,int,void*,void*);

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

#ifdef __cplusplus
}
#endif

typedef struct paramters_muscl {

  char limiter_type[_MAX_STRING_SIZE_]; 
  double eps; 
  double (*LimiterFunction) (double);

} MUSCLParameters;
int MUSCLInitialize(void*,void*);

typedef struct parameters_weno {
  int     mapped;           
  int     borges;           
  int     yc;               
  int     no_limiting;  
  double  eps;            
  double    p;                
  double  tol;          
  /* hybrid compact-WENO scheme related parameters
   * **References**:
   * + http://dx.doi.org/10.1006/jcph.2002.7021
   * + http://dx.doi.org/10.1016/j.jcp.2003.07.006
  */
  double  rc, 
          xi; 
  /* Arrays to save the WENO weights */
  double *w1, 
         *w2, 
         *w3;
  /* size and offset for the WENO weights arrays */
  int *offset ,
      size ;

} WENOParameters;

int WENOInitialize(void*,void*,char*,char*);
int WENOCleanup(void*, int);

/* define optimal weights */
#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) \
  { \
    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; \
    } \
  }

#define _WENOWeights_v_M_(w1,w2,w3,c1,c2,c3,m3,m2,m1,p1,p2,eps,N) \
  { \
    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; \
    } \
  }

#define _WENOWeights_v_M_Scalar_(w1,w2,w3,c1,c2,c3,m3,m2,m1,p1,p2,eps,idx) \
  { \
    /* 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; \
  }

#define _WENOWeights_v_Z_(w1,w2,w3,c1,c2,c3,m3,m2,m1,p1,p2,eps,N) \
  { \
    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; \
    } \
  }

#define _WENOWeights_v_YC_(w1,w2,w3,c1,c2,c3,m3,m2,m1,p1,p2,eps,N) \
  { \
    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; \
    } \
  }

#define _WENOWeights_v_YC_Scalar_(w1,w2,w3,c1,c2,c3,m3,m2,m1,p1,p2,eps,idx) \
  { \
    /* 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; \
  }

typedef struct compact_scheme {

  double *A, 
         *B, 
         *C, 
         *R;
  /* CRWENO scheme: buffer arrays for sending and receiving data */
  double *sendbuf, 
         *recvbuf; 
} CompactScheme;

int CompactSchemeInitialize(void*,void*,char*);
int CompactSchemeCleanup(void*);

#endif