numMethods/finiteDifferences.h

Go to the documentation of this file.

/**
 *  @file      finiteDifferences.h
 *  @brief     Macros to emit various finite differences (e.g., for approximating
               derivatives, or for the Keisser-Oliger term).
 *  @author    Mikica Kocic
 *  @copyright GNU General Public License (GPLv3).
 */

/////////////////////////////////////////////////////////////////////////////////////////
/** @defgroup g11 Numerical Methods                                                    */
/////////////////////////////////////////////////////////////////////////////////////////
                                                                                   /*@{*/
#ifndef CFDS_ORDER
    /** The order of the centered final difference scheme.
     */
    #define CFDS_ORDER 4
#endif

/// @todo A run-time CFDS_ORDER?

/////////////////////////////////////////////////////////////////////////////////////////
// Macros to emit the spatial derivatives (centered in space)
//
// Convention: if f() denotes a function, then f_r() denotes a spatial derivative,
// f_t() denotes a time derivative, f_rr() denotes a Real spatial derivative.
//
// @note The execution of the code where the spatial derivatives are
// converted to differences in Mathematica (and then simplified) is around two times
// faster than the code where the derivatives use the f_r and f_rr inline functions.

#if CFDS_ORDER == 2 // Second order central finite difference scheme

    /** `emitDerivative_r` is a macro that emits a method for the first derivative
     *  in `r` using the 2nd order central finite difference scheme.
     */
    #define emitDerivative_r( f ) inline Real f##_r( Int m, Int n ) \
        { return 0.5 * ( f(m,n+1) - f(m,n-1) ) * inv_delta_r; }

    /** `emitDerivative_rr` is a macro that emits a method for the second derivative
     *  in `r` using the 2nd order central finite difference scheme.
     */
    #define emitDerivative_rr( f ) inline Real f##_rr( Int m, Int n ) \
        { return ( f(m,n+1) - 2 * f(m,n) + f(m,n-1) ) * inv_delta_rr; }

#elif CFDS_ORDER == 4 // Fourth order central finite difference scheme

    /** `emitDerivative_r` is a macro that emits a method for the first derivative
     *  in `r` using the 4th order central finite difference scheme.
     */
    #define emitDerivative_r( f ) inline Real f##_r( Int m, Int n ) \
        { return ( - 1./12. * f(m,n+2) + 2./3.  * f(m,n+1) \
                   - 2./3.  * f(m,n-1) + 1./12. * f(m,n-2) ) * inv_delta_r; }

    /** `emitDerivative_rr` is a macro that emits a method for the second derivative
     *  in `r` using the 4th order central finite difference scheme.
     */
    #define emitDerivative_rr( f ) inline Real f##_rr( Int m, Int n ) \
        { return ( - 1./12.* f(m,n+2) + 4./3.  * f(m,n+1) - 5./2. * f(m,n) \
                   + 4./3. * f(m,n-1) - 1./12. * f(m,n-2) ) * inv_delta_rr; }

#elif CFDS_ORDER == 6 // Sixth order central finite difference scheme

    /** `emitDerivative_r` is a macro that emits a method for the first derivative
     *  in `r` using the 6th order central finite difference scheme.
     */
    #define emitDerivative_r( f ) inline Real f##_r( Int m, Int n ) \
        { return ( 1./60. * f(m,n+3) - 9./60. * f(m,n+2) + 45./60. * f(m,n+1) \
           - 45./60. * f(m,n-1) + 9./60. * f(m,n-2) - 1./60. * f(m,n-3) ) * inv_delta_r; }

    /** `emitDerivative_rr` is a macro that emits a method for the second derivative
     *  in `r` using the 6th order central finite difference scheme.
     */
    #define emitDerivative_rr( f ) inline Real f##_rr( Int m, Int n ) \
        { return ( 2./180. * f(m,n+3) - 27./180. * f(m,n+2) + 270./180. * f(m,n+1) \
                  - 490./180. * f(m,n) + 270./180. * f(m,n-1) - 27./180. * f(m,n-2) \
                  + 2./180. * f(m,n-3) ) * inv_delta_rr; }

#else

    #error "CFDS_ORDER must be either 2, 4, or 6"

#endif

/////////////////////////////////////////////////////////////////////////////////////////

/** extrapolate_R is an optimized version of 4th order in accuracy extrapolation
 *  using the 4th order Taylor expansion.
 */
#define extrapolate_R(f,m,n)  f(m,n) = \
    -   7./48. * f(m,n-8) +  209./144. * f(m,n-7) - 103./16. * f(m,n-6) \
    + 793./48. * f(m,n-5) - 3835./144. * f(m,n-4) + 439./16. * f(m,n-3) \
    - 281./16. * f(m,n-2) +  917./144. * f(m,n-1)

#define extrapolate_RGF(f,m,n)  GF(f,m,n) = \
    -   7./48. * GF(f,m,n-8) +  209./144. * GF(f,m,n-7) - 103./16. * GF(f,m,n-6) \
    + 793./48. * GF(f,m,n-5) - 3835./144. * GF(f,m,n-4) + 439./16. * GF(f,m,n-3) \
    - 281./16. * GF(f,m,n-2) +  917./144. * GF(f,m,n-1)

/** extrapolate_D is a linear extrapolation of the 4th order in accuracy
 *  (used for derivatives).
 */
#define extrapolate_D(f,m,n)  f(m,n) = \
    1./4. * f(m,n-5) - 4./3. * f(m,n-4) + 3 * f(m,n-3) - 4 * f(m,n-2) + 37./12. * f(m,n-1)

#define extrapolate_DGF(f,m,n)  GF(f,m,n) = \
    1./4. * GF(f,m,n-5) - 4./3. * GF(f,m,n-4) + 3 * GF(f,m,n-3) - 4 * GF(f,m,n-2) + 37./12. * GF(f,m,n-1)

/////////////////////////////////////////////////////////////////////////////////////////
    /** `KreissOligerTerm` is a macro that gives a Kreiss-Oliger dissipation term.
     *
     *  Coefficients d_k at gridpoints for the centred finite-difference discretisation
     *  of the derivative `partial_x^{2N}` (of the dissipative Kreiss-Oliger operator).
     *  <pre>
     *      ---------------------------------------------
     *       N  -4  -3   -2   -1    0   +1   +2  +3  +4
     *      ---------------------------------------------
     *       1                +1   -2   +1
     *       2           -1   +4   -6   +4   -1
     *       3      +1   -6  +15  -20  +15   -6  +1
     *       4  -1  +8  -28  +56  -70  +56  -28  +8  -1
     *      ---------------------------------------------  </pre>
     *
     *  Here, `N >= 1` is an integer and `D_h^{2N}` is a centered difference operator
     *  approximating a partial spatial derivative of order `2N`, e.g.,
     *  a second- (`N = 1`) or fourth-order (`N = 2`) spatial derivative.
     *
     *  To be subtracted from the evolution equation:
     *
     *      `epsilon (-1)^N  delta_t {delta_x}^{2N-1}  D_h^{2N}`
     *
     *  where:
     *
     *      `D_h^{2N} = (2 delta_x)^{-2N} * ( sum_k d_k f[n+k] )`
     *
     *  For `N = 2`, subtract from the evolution equation `u_t`:
     *
     *      `dissip(f,dt) = epsilon 2^{-2N} * dt/delta_x * ( sum_k d_k f[n+k] )`
     */
#if CFDS_ORDER == 2
        #define KreissOligerTerm(f,dt) \
            ( - ( f(m,n-1) - 2* f(m,n) + f(m,n+1) ) * dissip_delta_r * (dt) )
        #define KreissOligerTermGF(f,dt) \
            ( - ( GF(f,m,n-1) - 2* GF(f,m,n) + GF(f,m,n+1) ) * dissip_delta_r * (dt) )
#else
    #if CFDS_ORDER == 4
        #define KreissOligerTerm(f,dt) \
            ( ( f(m,n-2) - 4* f(m,n-1) + 6* f(m,n) - 4* f(m,n+1) + f(m,n+2) ) \
               * dissip_delta_r * (dt) )
        #define KreissOligerTermGF(f,dt) \
            ( ( GF(f,m,n-2) - 4* GF(f,m,n-1) + 6* GF(f,m,n) - 4* GF(f,m,n+1) + GF(f,m,n+2) ) \
               * dissip_delta_r * (dt) )
    #else
        #define KreissOligerTerm(f,dt) \
            ( - ( f(m,n-3) - 6* f(m,n-2) + 15* f(m,n-1) - 20* f(m,n) + 15* f(m,n+1) \
               - 6* f(m,n+2) + f(m,n+3) ) * dissip_delta_r * (dt) )
        #define KreissOligerTermGF(f,dt) \
            ( - ( GF(f,m,n-3) - 6* GF(f,m,n-2) + 15* GF(f,m,n-1) - 20* GF(f,m,n) + 15* GF(f,m,n+1) \
               - 6* GF(f,m,n+2) + GF(f,m,n+3) ) * dissip_delta_r * (dt) )
    #endif
#endif
                                                                                   /*@}*/
/////////////////////////////////////////////////////////////////////////////////////////