kernel volume computation ok

kernels.cl

@@ -205,7 +205,7 @@ void DGVolume(int ie,                    // 1: macrocel index
   __constant float *physnode = physnodes + ie * 3 * _NHEXA;
 
   int ix = get_global_id(0);
-  int iy = get_global_id(1);
+  int iy = get_global_id(1); 
   int iz = get_global_id(2);
  
   int woffset = ie * _M * _NPG_MACROCELL;
@@ -292,3 +292,123 @@ void DGVolume(int ie,                    // 1: macrocel index
 
 }
 
+// Compute the surface terms inside one macrocell
+__kernel
+void DGFlux(int ie,                      // 1: macrocell index
+	    __constant float *physnodes, // 3: macrocell nodes
+	    __global   float *wn,       // 4: field values
+	    __global   float *dtwn     // 5: time derivative
+	    )
+{
+
+  __constant float *physnode = physnodes + ie * 3 * _NHEXA;
+
+  int ifx = get_global_id(0); // face id in x direction: 0.._RAF
+  int iy = get_global_id(1);// pg id in other directions: 0.._NPG_ROW
+  int iz = get_global_id(2);
+ 
+  int woffset = ie * _M * _NPG_MACROCELL;
+
+  // subcell id
+  int icL[3];
+  icL[0] = ifx - 1;
+  icL[1] = iy / _RAF;
+  icL[2] = iz / _RAF;
+
+  int pL[3];
+  pL[0] = _RAF - 1;
+  pL[1] = iy % _RAF;
+  pL[2] = iz % _RAF;
+
+  int icR[3];
+
+  icR[0] = ifx;
+  icR[1] = icL[1];
+  icR[2] = icL[2];
+
+  int pR[3];
+  pR[0] = 0;
+  pR[1] = iy % _RAF;
+  pR[2] = iz % _RAF;
+  
+  float wL[_M], wR[_M];
+  int ipgL = -1;
+  if (ifx > 0)
+    ipgL = pL[0] + _RAF * icL[0] + _NPG_ROW * (iy + _NPG_ROW * iz);
+  int ipgR = -1;
+  if (ifx < _RAF)
+    ipgR = pR[0] + _RAF * icR[0] + _NPG_ROW * (iy + _NPG_ROW * iz);
+  
+  for(int iv = 0; iv < _M; iv++) {
+    if (ipgL >= 0){
+      int imemL =  Varindex(ipgL, iv) + woffset;
+      wL[iv] = wn[imemL];
+    }
+    if (ipgR >= 0){
+      int imemR =  Varindex(ipgR, iv) + woffset;
+      wR[iv] = wn[imemR];
+    }
+  }
+
+  float hx = 1.0 / _RAF;
+  float hy = 1.0 / _RAF;
+  float hz = 1.0 / _RAF;
+
+  int offset[3] = {gauss_lob_offset[_DEG] + pL[0],
+  		   gauss_lob_offset[_DEG] + pL[1],
+  		   gauss_lob_offset[_DEG] + pL[2]};
+
+  float x = hx * (icL[0] + gauss_lob_point[offset[0]]);
+  float y = hy * (icL[1] + gauss_lob_point[offset[1]]);
+  float z = hz * (icL[2] + gauss_lob_point[offset[2]]);
+
+  float codtau[3][3];
+  {
+    float dtau[3][3];
+    get_dtau(x, y, z, physnode, dtau);
+    codtau[0][0] =  dtau[1][1] * dtau[2][2] - dtau[1][2] * dtau[2][1];
+    codtau[0][1] = -dtau[1][0] * dtau[2][2] + dtau[1][2] * dtau[2][0];
+    codtau[0][2] =  dtau[1][0] * dtau[2][1] - dtau[1][1] * dtau[2][0];
+    codtau[1][0] = -dtau[0][1] * dtau[2][2] + dtau[0][2] * dtau[2][1];
+    codtau[1][1] =  dtau[0][0] * dtau[2][2] - dtau[0][2] * dtau[2][0];
+    codtau[1][2] = -dtau[0][0] * dtau[2][1] + dtau[0][1] * dtau[2][0];
+    codtau[2][0] =  dtau[0][1] * dtau[1][2] - dtau[0][2] * dtau[1][1];
+    codtau[2][1] = -dtau[0][0] * dtau[1][2] + dtau[0][2] * dtau[1][0];
+    codtau[2][2] =  dtau[0][0] * dtau[1][1] - dtau[0][1] * dtau[1][0];
+  }
+
+  float h1h2 = 1.0 / _RAF / _RAF;
+  float vnds[3];
+  vnds[0] = codtau[0][0] * h1h2;
+  vnds[1] = codtau[1][0] * h1h2;
+  vnds[2] = codtau[2][0] * h1h2;
+
+  for(int iv = 0; iv < _M; iv++) {
+    int imemL =  Varindex(ipgL, iv) + woffset;
+    wL[iv] = wn[imemL];
+    int imemR =  Varindex(ipgR, iv) + woffset;
+    wR[iv] = wn[imemR];
+  }
+  
+  // TODO: test ipgL and ipgR for Boundary terms !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
+  float flux[_M];
+  NumFlux(wL, wR, vnds, flux); //
+  
+  float wpgs = wglop(_DEG, pL[1]) * wglop(_DEG, pL[2]);
+  
+  for(int iv = 0; iv < _M; ++iv) {
+    if (ipgL >= 0){
+      int imemL =  Varindex(ipgL, iv) + woffset;
+      dtwn[imemL] -= flux[iv] * wpgs;
+    }
+    
+    
+    if (ipgR >= 0){
+      int imemR =  Varindex(ipgR, iv) + woffset;
+      dtwn[imemR] += flux[iv] * wpgs;
+    }
+  }
+  
+}
+
+