#include "sp.h"

//---------------------------------------------------------------------
// compute the reciprocal of density, and the kinetic energy, 
// and the speed of sound. 
//---------------------------------------------------------------------
__kernel void compute_rhs1(__global const double *g_u,
                           __global double *g_us,
                           __global double *g_vs,
                           __global double *g_ws,
                           __global double *g_qs,
                           __global double *g_rho_i,
                           __global double *g_speed,
                           __global double *g_square,
                           int gp0,
                           int gp1,
                           int gp2)
{
  int i, j, k;
  double aux, rho_inv;
  double p_u[5];
  double p_square;

  k = get_global_id(2);
  j = get_global_id(1);
  i = get_global_id(0);
  if (k >= gp2 || j >= gp1 || i >= gp0) return;

  __global double (*u)[JMAXP+1][IMAXP+1][5] = 
    (__global double (*)[JMAXP+1][IMAXP+1][5])g_u;
  __global double (*us)[JMAXP+1][IMAXP+1] = 
    (__global double (*)[JMAXP+1][IMAXP+1])g_us;
  __global double (*vs)[JMAXP+1][IMAXP+1] = 
    (__global double (*)[JMAXP+1][IMAXP+1])g_vs;
  __global double (*ws)[JMAXP+1][IMAXP+1] = 
    (__global double (*)[JMAXP+1][IMAXP+1])g_ws;
  __global double (*qs)[JMAXP+1][IMAXP+1] = 
    (__global double (*)[JMAXP+1][IMAXP+1])g_qs;
  __global double (*rho_i)[JMAXP+1][IMAXP+1] = 
    (__global double (*)[JMAXP+1][IMAXP+1])g_rho_i;
  __global double (*speed)[JMAXP+1][IMAXP+1] = 
    (__global double (*)[JMAXP+1][IMAXP+1])g_speed;
  __global double (*square)[JMAXP+1][IMAXP+1] = 
    (__global double (*)[JMAXP+1][IMAXP+1])g_square;

  p_u[0] = u[k][j][i][0];
  p_u[1] = u[k][j][i][1];
  p_u[2] = u[k][j][i][2];
  p_u[3] = u[k][j][i][3];
  p_u[4] = u[k][j][i][4];

  rho_inv = 1.0/p_u[0];
  rho_i[k][j][i] = rho_inv;
  us[k][j][i] = p_u[1] * rho_inv;
  vs[k][j][i] = p_u[2] * rho_inv;
  ws[k][j][i] = p_u[3] * rho_inv;
  p_square = 0.5* (p_u[1]*p_u[1] + p_u[2]*p_u[2] + p_u[3]*p_u[3]) * rho_inv;
  square[k][j][i] = p_square;
  qs[k][j][i] = p_square * rho_inv;
  //-------------------------------------------------------------------
  // (don't need speed and ainx until the lhs computation)
  //-------------------------------------------------------------------
  aux = c1c2*rho_inv* (p_u[4] - p_square);
  speed[k][j][i] = sqrt(aux);
}

