text: // NBN: 2020/11/16
  // NBN: 2020/11/20 - "do-nothing" outflow
  
  
  #define eulerInitialConditions2D(gamma, t, cx, cy, x, y, r, u, v, p)  \
    {                 \
      *(r) = p_rbar;  \
      *(u) = p_ubar;  \
      *(v) = p_vbar;  \
      *(p) = p_pbar;  \
    }
  
  // Body force
  #define eulerBodyForce2D(gamma, t, x, y, r, u, v, p, fx, fy)	\
    {               \
      *(fx) = 0.0;  \
      *(fy) = 0.0;  \
    }
  
  // Boundary conditions (actually run with periodic) 
  /* wall 1, inflow 2, outflow 3, x-slip 4, y-slip 5 */
  #define p_WALL 1
  #define p_INFLOW 2
  #define p_OUTFLOW 3
  
  #define eulerBoundaryConditions2D(bc, gamma,          \
                                    t, x, y, nx, ny,    \
                                    rM, uM, vM, pM,     \
                                    rB, uB, vB, pB)     \
    {                                                   \
      switch(bc){                          \
      case p_WALL:{ /* wall */             \
        *(rB) = rM;                        \
        *(uB) = uM-2.*(nx*uM + ny*vM)*nx;  \
        *(vB) = vM-2.*(nx*uM + ny*vM)*ny;  \
        *(pB) = pM;                        \
        break;                             \
      }                                    \
      case p_INFLOW:{ /* inflow */    \
        *(rB) = p_rbar;               \
        *(uB) = p_ubar;               \
        *(vB) = p_vbar;               \
        *(pB) = p_pbar;               \
        break;                        \
      }                               \
      case p_OUTFLOW:{ /* outflow */  \
        *(rB) = rM;                   \
        *(uB) = uM;                   \
        *(vB) = vM;                   \
        *(pB) = pM;                   \
        break;                        \
      }                               \
      default:{   \
        ;         \
      }           \
      }           \
    }
  
  #define eulerGradientBoundaryConditions2D(bc, gamma,            \
                                            t, x, y, nx, ny,      \
                                            v1M, v2M, v3M, v4M,   \
                                            v1B, v2B, v3B, v4B)   \
    { /* only apply bcs for fields 1:3 */   \
      switch(bc){                     \
      case p_WALL:{ /* wall */        \
        *(v2B) =  -v2M;               \
        *(v3B) =  -v3M;               \
        *(v4B) =  v4M;                \
        break;                        \
      }                               \
      case p_INFLOW:{ /* inflow */    \
        *(v2B) =  v2M;                \
        *(v3B) =  v3M;                \
        *(v4B) =  v4M;                \
        break;                        \
      }                               \
      case p_OUTFLOW:{ /* outflow */  \
        *(v2B) =  v2M;      \
        *(v3B) =  v3M;      \
        *(v4B) =  v4M;      \
        break;      \
      }             \
      default:{     \
        ;           \
      }             \
      }             \
    }
  
  
  #define eulerDivergenceBoundaryConditions2D(bc, gamma, \
      t, x, y, nx, ny,                  \
      dv1dnM, dv2dnM, dv3dnM, dv4dnM,   \
      dv1dnB, dv2dnB, dv3dnB, dv4dnB)   \
    {                                   \
      switch(bc){                       \
      case p_WALL:{ /* wall */  \
        *(dv2dnB) =  dv2dnM;    \
        *(dv3dnB) =  dv3dnM;    \
        *(dv4dnB) = -dv4dnM;    \
        break;                      \
      }                             \
      case p_INFLOW:{ /* inflow */  \
        *(dv2dnB) = -dv2dnM;        \
        *(dv3dnB) = -dv3dnM;        \
        *(dv4dnB) = -dv4dnM;        \
        break;                      \
      }                               \
      case p_OUTFLOW:{ /* outflow */  \
        *(dv2dnB) = dv2dnM;           \
        *(dv3dnB) = dv3dnM;           \
        *(dv4dnB) = dv4dnM;           \
        break;      \
      }             \
      default:{			\
        ;           \
      }             \
      }             \
    }
  