kernel flux and volume written

kernels.cl

@@ -145,23 +145,65 @@ int Varindex(int ipg, int iv) {
 
 void NumFlux(float *wL, float *wR, float *vnorm, float *flux)
 {
-  const float transport_v2d[] = {sqrt(0.5), sqrt(0.5), 0};
-  //const real transport_v2d[] = {1,0, 0};
-  float vn
-    = transport_v2d[0] * vnorm[0]
-    + transport_v2d[1] * vnorm[1]
-    + transport_v2d[2] * vnorm[2];
-  float vnp = vn > 0 ? vn : 0;
+  float vn = sqrt(0.5) * vnorm[0] + sqrt(0.5) * vnorm[1];
+  float vnp = vn > 0.0 ? vn : 0.0;
   float vnm = vn - vnp;
   flux[0] = vnp * wL[0] + vnm * wR[0];
-  /* if (fabs(vnorm[2])>1e-6) { */
-  /*   printf("vnds %lf %lf %lf \n", vnorm[0], vnorm[1], vnorm[2]); */
-  /* } */
-  // verify that 2d computations are actually
-  // activated
-  //assert(fabs(vnorm[2]) < 1e-8);
+  vn = -vn;
+  vnp = vn > 0.0 ? vn : 0.0;
+  vnm = vn - vnp;
+  flux[1] = vnp * wL[1] + vnm * wR[1];
 }
+void ImposedData(const float *x, const float t, float* w)
+{
+  float vx = sqrt(0.5) * x[0] + sqrt(0.5) * x[1];
+  float xx = vx - t;
+  w[0] = xx * xx;
+  w[0] = sin(xx);
+  xx = vx + t;
+  w[1] = xx * xx;
+  w[1] = sin(xx);
+}
+
+
+void BoundaryFlux(float *x, float t,
+			    float *wL, float *vnorm,
+			    float *flux)
+{
+  float wR[2];
+  ImposedData(x, t, wR);
+  NumFlux(wL, wR, vnorm, flux);
+}
+
+void get_tau(float x, float y, float z,
+	     __constant float *p, float xphy[]){
 
+  float phi[_NHEXA];
+  
+  float t1 = -1 + z;
+  float t2 = -1 + y;
+  float t4 = -1 + x;
+  float t8 = x * y;
+  phi[0] = -t1 * t2 * t4;
+  phi[1] = x * t1 * t2;
+  phi[2] = -t8 * t1;
+  phi[3] = y * t1 * t4;
+  phi[4] = z * t2 * t4;
+  phi[5] = -x * z * t2;
+  phi[6] = t8 * z;
+  phi[7] = -y * z * t4;
+
+  xphy[0] = 0;
+  xphy[1] = 0;
+  xphy[2] = 0;
+ 
+ for(int ii = 0; ii < _NHEXA; ii++){
+   xphy[0] += p[ii * 3 + 0] * phi[ii]; 
+   xphy[1] += p[ii * 3 + 1] * phi[ii]; 
+   xphy[2] += p[ii * 3 + 2] * phi[ii];
+ }
+ 
+}
 
 void get_dtau(float x, float y, float z,
 	      __constant float *p, float dtau[][3]){
@@ -296,6 +338,7 @@ void DGVolume(int ie,                    // 1: macrocel index
 __kernel
 void DGFlux(int ie,                      // 1: macrocell index
 	    __constant float *physnodes, // 3: macrocell nodes
+	    __constant float *tnow,
 	    __global   float *wn,       // 4: field values
 	    __global   float *dtwn     // 5: time derivative
 	    )
@@ -331,25 +374,6 @@ void DGFlux(int ie,                      // 1: macrocell index
   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;
@@ -383,13 +407,39 @@ void DGFlux(int ie,                      // 1: macrocell index
   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];
+
+  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);
+
+  float xphy[3];
+  if (ipgL < 0 || ipgR < 0){
+    get_tau(x, y, z, physnode, xphy);
   }
+
   
+  float wL[_M], wR[_M];
+  float tim = *tnow;
+  for(int iv = 0; iv < _M; iv++) {
+    if (ipgL >= 0){
+      int imemL =  Varindex(ipgL, iv) + woffset;
+      wL[iv] = wn[imemL];
+    } else {
+      ImposedData(xphy, tim, wL);
+    }
+    if (ipgR >= 0){
+      int imemR =  Varindex(ipgR, iv) + woffset;
+      wR[iv] = wn[imemR];
+    } else {
+      ImposedData(xphy, tim, wR);      
+    }
+  }
+
+
+
   // TODO: test ipgL and ipgR for Boundary terms !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
   float flux[_M];
   NumFlux(wL, wR, vnds, flux); //