Merge pull request #1866 from CEED/jeremy/fix-mpm

jeremylt · web-flow · commit 66c2c3811ed2 · 2025-07-18T07:57:45.000-06:00
gpu - fix gen AtPoints transpose
diff --git a/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h b/include/ceed/jit-source/cuda/cuda-shared-basis-tensor-at-points-templates.h
@@ -194,7 +194,7 @@ inline __device__ void InterpTransposeAtPoints2d(SharedData_Cuda &data, const Ce
       for (CeedInt j = 0; j < Q_1D; j++) {
         const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-        atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+        if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
       }
     }
     // Pull from shared to register
@@ -270,7 +270,7 @@ inline __device__ void GradTransposeAtPoints2d(SharedData_Cuda &data, const Ceed
         for (CeedInt j = 0; j < Q_1D; j++) {
           const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-          atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+          if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
         }
       }
     }
@@ -355,7 +355,7 @@ inline __device__ void InterpTransposeAtPoints3d(SharedData_Cuda &data, const Ce
         for (CeedInt j = 0; j < Q_1D; j++) {
           const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-          atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+          if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
         }
       }
       // Pull from shared to register
@@ -439,11 +439,12 @@ inline __device__ void GradTransposeAtPoints3d(SharedData_Cuda &data, const Ceed
         if (dim == 1) ChebyshevDerivativeAtPoint<Q_1D>(r_X[1], chebyshev_x);
         else ChebyshevPolynomialsAtPoint<Q_1D>(r_X[1], chebyshev_x);
         const CeedScalar zz  = dim == 2 ? dz : z;
-        const CeedScalar r_u = p < NUM_POINTS ? r_U[comp + dim * NUM_COMP] : 0.0;
+        const CeedScalar r_u = (p < NUM_POINTS) ? r_U[comp + dim * NUM_COMP] : 0.0;
 
         for (CeedInt i = 0; i < Q_1D; i++) {
           buffer[i] = chebyshev_x[i] * r_u * zz;
         }
+
         // Contract x direction
         if (dim == 0) ChebyshevDerivativeAtPoint<Q_1D>(r_X[0], chebyshev_x);
         else ChebyshevPolynomialsAtPoint<Q_1D>(r_X[0], chebyshev_x);
@@ -454,7 +455,7 @@ inline __device__ void GradTransposeAtPoints3d(SharedData_Cuda &data, const Ceed
           for (CeedInt j = 0; j < Q_1D; j++) {
             const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-            atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+            if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
           }
         }
       }
diff --git a/include/ceed/jit-source/hip/hip-shared-basis-tensor-at-points-templates.h b/include/ceed/jit-source/hip/hip-shared-basis-tensor-at-points-templates.h
@@ -195,7 +195,7 @@ inline __device__ void InterpTransposeAtPoints2d(SharedData_Hip &data, const Cee
       for (CeedInt j = 0; j < Q_1D; j++) {
         const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-        atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+        if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
       }
     }
     // Pull from shared to register
@@ -271,7 +271,7 @@ inline __device__ void GradTransposeAtPoints2d(SharedData_Hip &data, const CeedI
         for (CeedInt j = 0; j < Q_1D; j++) {
           const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-          atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+          if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
         }
       }
     }
@@ -356,7 +356,7 @@ inline __device__ void InterpTransposeAtPoints3d(SharedData_Hip &data, const Cee
         for (CeedInt j = 0; j < Q_1D; j++) {
           const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-          atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+          if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
         }
       }
       // Pull from shared to register
@@ -455,7 +455,7 @@ inline __device__ void GradTransposeAtPoints3d(SharedData_Hip &data, const CeedI
           for (CeedInt j = 0; j < Q_1D; j++) {
             const CeedInt jj = (j + data.t_id_x) % Q_1D;
 
-            atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
+            if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);
           }
         }
       }

Original file line number	Diff line number	Diff line change
`@@ -194,7 +194,7 @@ inline __device__ void InterpTransposeAtPoints2d(SharedData_Cuda &data, const Ce`
`194`	`194`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`195`	`195`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`196`	`196`
`197`		`- atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`197`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`198`	`198`	`}`
`199`	`199`	`}`
`200`	`200`	`// Pull from shared to register`
`@@ -270,7 +270,7 @@ inline __device__ void GradTransposeAtPoints2d(SharedData_Cuda &data, const Ceed`
`270`	`270`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`271`	`271`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`272`	`272`
`273`		`- atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`273`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`274`	`274`	`}`
`275`	`275`	`}`
`276`	`276`	`}`
`@@ -355,7 +355,7 @@ inline __device__ void InterpTransposeAtPoints3d(SharedData_Cuda &data, const Ce`
`355`	`355`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`356`	`356`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`357`	`357`
`358`		`- atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`358`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`359`	`359`	`}`
`360`	`360`	`}`
`361`	`361`	`// Pull from shared to register`
`@@ -439,11 +439,12 @@ inline __device__ void GradTransposeAtPoints3d(SharedData_Cuda &data, const Ceed`
`439`	`439`	`if (dim == 1) ChebyshevDerivativeAtPoint<Q_1D>(r_X[1], chebyshev_x);`
`440`	`440`	`else ChebyshevPolynomialsAtPoint<Q_1D>(r_X[1], chebyshev_x);`
`441`	`441`	`const CeedScalar zz = dim == 2 ? dz : z;`
`442`		`- const CeedScalar r_u = p < NUM_POINTS ? r_U[comp + dim * NUM_COMP] : 0.0;`
	`442`	`+ const CeedScalar r_u = (p < NUM_POINTS) ? r_U[comp + dim * NUM_COMP] : 0.0;`
`443`	`443`
`444`	`444`	`for (CeedInt i = 0; i < Q_1D; i++) {`
`445`	`445`	`buffer[i] = chebyshev_x[i] * r_u * zz;`
`446`	`446`	`}`
	`447`	`+`
`447`	`448`	`// Contract x direction`
`448`	`449`	`if (dim == 0) ChebyshevDerivativeAtPoint<Q_1D>(r_X[0], chebyshev_x);`
`449`	`450`	`else ChebyshevPolynomialsAtPoint<Q_1D>(r_X[0], chebyshev_x);`
`@@ -454,7 +455,7 @@ inline __device__ void GradTransposeAtPoints3d(SharedData_Cuda &data, const Ceed`
`454`	`455`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`455`	`456`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`456`	`457`
`457`		`- atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`458`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd_block(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`458`	`459`	`}`
`459`	`460`	`}`
`460`	`461`	`}`
Original file line number	Diff line number	Diff line change
`@@ -195,7 +195,7 @@ inline __device__ void InterpTransposeAtPoints2d(SharedData_Hip &data, const Cee`
`195`	`195`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`196`	`196`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`197`	`197`
`198`		`- atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`198`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`199`	`199`	`}`
`200`	`200`	`}`
`201`	`201`	`// Pull from shared to register`
`@@ -271,7 +271,7 @@ inline __device__ void GradTransposeAtPoints2d(SharedData_Hip &data, const CeedI`
`271`	`271`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`272`	`272`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`273`	`273`
`274`		`- atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`274`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`275`	`275`	`}`
`276`	`276`	`}`
`277`	`277`	`}`
`@@ -356,7 +356,7 @@ inline __device__ void InterpTransposeAtPoints3d(SharedData_Hip &data, const Cee`
`356`	`356`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`357`	`357`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`358`	`358`
`359`		`- atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`359`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`360`	`360`	`}`
`361`	`361`	`}`
`362`	`362`	`// Pull from shared to register`
`@@ -455,7 +455,7 @@ inline __device__ void GradTransposeAtPoints3d(SharedData_Hip &data, const CeedI`
`455`	`455`	`for (CeedInt j = 0; j < Q_1D; j++) {`
`456`	`456`	`const CeedInt jj = (j + data.t_id_x) % Q_1D;`
`457`	`457`
`458`		`- atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
	`458`	`+ if (data.t_id_x < Q_1D && data.t_id_y < Q_1D) atomicAdd(&data.slice[jj + ii * Q_1D], chebyshev_x[jj] * buffer[ii]);`
`459`	`459`	`}`
`460`	`460`	`}`
`461`	`461`	`}`