1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #define CEED_DEBUG_COLOR 12 18 19 #include <ceed/ceed.h> 20 #include <ceed/backend.h> 21 #include <cuda_runtime.h> 22 #include <iostream> 23 #include <sstream> 24 #include "ceed-cuda-gen.h" 25 #include "../cuda-shared/ceed-cuda-shared.h" 26 27 static const char *atomicAdd = QUOTE( 28 //------------------------------------------------------------------------------ 29 // Atomic add, for older CUDA 30 //------------------------------------------------------------------------------ 31 __device__ double atomicAdd(double *address, double val) { 32 unsigned long long int *address_as_ull = (unsigned long long int *)address; 33 unsigned long long int old = *address_as_ull, assumed; 34 do { 35 assumed = old; 36 old = 37 atomicCAS(address_as_ull, assumed, 38 __double_as_longlong(val + 39 __longlong_as_double(assumed))); 40 // Note: uses integer comparison to avoid hang in case of NaN 41 // (since NaN != NaN) 42 } while (assumed != old); 43 return __longlong_as_double(old); 44 } 45 ); 46 47 static const char *deviceFunctions = QUOTE( 48 49 //------------------------------------------------------------------------------ 50 // Typedefs 51 //------------------------------------------------------------------------------ 52 typedef struct { const CeedScalar* in[16]; CeedScalar* out[16]; } CudaFields; 53 typedef struct { CeedInt* in[16]; CeedInt* out[16]; } CudaFieldsInt; 54 55 typedef struct { 56 CeedInt tidx; 57 CeedInt tidy; 58 CeedInt tidz; 59 CeedInt tid; 60 CeedScalar* slice; 61 } BackendData; 62 63 //------------------------------------------------------------------------------ 64 // Load matrices for basis actions 65 //------------------------------------------------------------------------------ 66 template <int P, int Q> 67 inline __device__ void loadMatrix(BackendData& data, const CeedScalar* d_B, CeedScalar* B) { 68 for (CeedInt i = data.tid; i < P*Q; i += blockDim.x*blockDim.y*blockDim.z) 69 B[i] = d_B[i]; 70 } 71 72 //------------------------------------------------------------------------------ 73 // 1D 74 //------------------------------------------------------------------------------ 75 76 //------------------------------------------------------------------------------ 77 // L-vector -> E-vector, offsets provided 78 //------------------------------------------------------------------------------ 79 template <int NCOMP, int COMPSTRIDE, int P1d> 80 inline __device__ void readDofsOffset1d(BackendData& data, const CeedInt nnodes, const CeedInt elem, const CeedInt* indices, const CeedScalar* d_u, CeedScalar* r_u) { 81 if (data.tidx < P1d) { 82 const CeedInt node = data.tidx; 83 const CeedInt ind = indices[node + elem * P1d]; 84 for (CeedInt comp = 0; comp < NCOMP; ++comp) 85 r_u[comp] = d_u[ind + COMPSTRIDE * comp]; 86 } 87 } 88 89 //------------------------------------------------------------------------------ 90 // L-vector -> E-vector, strided 91 //------------------------------------------------------------------------------ 92 template <int NCOMP, int P1d, int STRIDES_NODE, int STRIDES_COMP, int STRIDES_ELEM> 93 inline __device__ void readDofsStrided1d(BackendData& data, const CeedInt elem, const CeedScalar* d_u, CeedScalar* r_u) { 94 if (data.tidx < P1d) { 95 const CeedInt node = data.tidx; 96 const CeedInt ind = node * STRIDES_NODE + elem * STRIDES_ELEM; 97 for (CeedInt comp = 0; comp < NCOMP; ++comp) 98 r_u[comp] = d_u[ind + comp * STRIDES_COMP]; 99 } 100 } 101 102 //------------------------------------------------------------------------------ 103 // E-vector -> L-vector, offsets provided 104 //------------------------------------------------------------------------------ 105 template <int NCOMP, int COMPSTRIDE, int P1d> 106 inline __device__ void writeDofsOffset1d(BackendData& data, const CeedInt nnodes, const CeedInt elem, const CeedInt* indices, const CeedScalar* r_v, CeedScalar* d_v) { 107 if (data.tidx < P1d) { 108 const CeedInt node = data.tidx; 109 const CeedInt ind = indices[node + elem * P1d]; 110 for (CeedInt comp = 0; comp < NCOMP; ++comp) 111 atomicAdd(&d_v[ind + COMPSTRIDE * comp], r_v[comp]); 112 } 113 } 114 115 //------------------------------------------------------------------------------ 116 // E-vector -> L-vector, strided 117 //------------------------------------------------------------------------------ 118 template <int NCOMP, int P1d, int STRIDES_NODE, int STRIDES_COMP, int STRIDES_ELEM> 119 inline __device__ void writeDofsStrided1d(BackendData& data, const CeedInt elem, const CeedScalar* r_v, CeedScalar* d_v) { 120 if (data.tidx < P1d) { 121 const CeedInt node = data.tidx; 122 const CeedInt ind = node * STRIDES_NODE + elem * STRIDES_ELEM; 123 for (CeedInt comp = 0; comp < NCOMP; ++comp) 124 d_v[ind + comp * STRIDES_COMP] += r_v[comp]; 125 } 126 } 127 128 //------------------------------------------------------------------------------ 129 // 1D tensor contraction x 130 //------------------------------------------------------------------------------ 131 template <int NCOMP, int P1d, int Q1d> 132 inline __device__ void ContractX1d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 133 data.slice[data.tidx] = *U; 134 __syncthreads(); 135 *V = 0.0; 136 if (data.tidx < Q1d) 137 for (CeedInt i = 0; i < P1d; ++i) 138 *V += B[i + data.tidx*P1d] * data.slice[i]; // Contract x direction 139 __syncthreads(); 140 } 141 142 //------------------------------------------------------------------------------ 143 // 1D transpose tensor contraction x 144 //------------------------------------------------------------------------------ 145 template <int NCOMP, int P1d, int Q1d> 146 inline __device__ void ContractTransposeX1d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 147 data.slice[data.tidx] = *U; 148 __syncthreads(); 149 *V = 0.0; 150 if (data.tidx < P1d) 151 for (CeedInt i = 0; i < Q1d; ++i) 152 *V += B[data.tidx + i*P1d] * data.slice[i]; // Contract x direction 153 __syncthreads(); 154 } 155 156 //------------------------------------------------------------------------------ 157 // 1D interpolate to quadrature points 158 //------------------------------------------------------------------------------ 159 template <int NCOMP, int P1d, int Q1d> 160 inline __device__ void interp1d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, CeedScalar *__restrict__ r_V) { 161 for (CeedInt comp = 0; comp < NCOMP; comp++) 162 ContractX1d<NCOMP, P1d, Q1d>(data, r_U + comp, c_B, r_V + comp); 163 } 164 165 //------------------------------------------------------------------------------ 166 // 1D interpolate transpose 167 //------------------------------------------------------------------------------ 168 template <int NCOMP, int P1d, int Q1d> 169 inline __device__ void interpTranspose1d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, CeedScalar *__restrict__ r_V) { 170 for (CeedInt comp=0; comp<NCOMP; comp++) 171 ContractTransposeX1d<NCOMP, P1d, Q1d>(data, r_U + comp, c_B, r_V + comp); 172 } 173 174 //------------------------------------------------------------------------------ 175 // 1D derivatives at quadrature points 176 //------------------------------------------------------------------------------ 177 template <int NCOMP, int P1d, int Q1d> 178 inline __device__ void grad1d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 179 for (CeedInt comp = 0; comp < NCOMP; comp++) 180 ContractX1d<NCOMP, P1d, Q1d>(data, r_U + comp, c_G, r_V + comp); 181 } 182 183 //------------------------------------------------------------------------------ 184 // 1D derivatives transpose 185 //------------------------------------------------------------------------------ 186 template <int NCOMP, int P1d, int Q1d> 187 inline __device__ void gradTranspose1d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 188 for (CeedInt comp = 0; comp < NCOMP; comp++) 189 ContractTransposeX1d<NCOMP, P1d, Q1d>(data, r_U + comp, c_G, r_V + comp); 190 } 191 192 //------------------------------------------------------------------------------ 193 // 2D 194 //------------------------------------------------------------------------------ 195 196 //------------------------------------------------------------------------------ 197 // L-vector -> E-vector, offsets provided 198 //------------------------------------------------------------------------------ 199 template <int NCOMP, int COMPSTRIDE, int P1d> 200 inline __device__ void readDofsOffset2d(BackendData& data, const CeedInt nnodes, const CeedInt elem, const CeedInt* indices, const CeedScalar* d_u, CeedScalar* r_u) { 201 if (data.tidx < P1d && data.tidy < P1d) { 202 const CeedInt node = data.tidx + data.tidy*P1d; 203 const CeedInt ind = indices[node + elem * P1d*P1d]; 204 for (CeedInt comp = 0; comp < NCOMP; ++comp) 205 r_u[comp] = d_u[ind + COMPSTRIDE * comp]; 206 } 207 } 208 209 //------------------------------------------------------------------------------ 210 // L-vector -> E-vector, strided 211 //------------------------------------------------------------------------------ 212 template <int NCOMP, int P1d, int STRIDES_NODE, int STRIDES_COMP, int STRIDES_ELEM> 213 inline __device__ void readDofsStrided2d(BackendData& data, const CeedInt elem, const CeedScalar* d_u, CeedScalar* r_u) { 214 if (data.tidx < P1d && data.tidy < P1d) { 215 const CeedInt node = data.tidx + data.tidy*P1d; 216 const CeedInt ind = node * STRIDES_NODE + elem * STRIDES_ELEM; 217 for (CeedInt comp = 0; comp < NCOMP; ++comp) 218 r_u[comp] = d_u[ind + comp * STRIDES_COMP]; 219 } 220 } 221 222 //------------------------------------------------------------------------------ 223 // E-vector -> L-vector, offsets provided 224 //------------------------------------------------------------------------------ 225 template <int NCOMP, int COMPSTRIDE, int P1d> 226 inline __device__ void writeDofsOffset2d(BackendData& data, const CeedInt nnodes, const CeedInt elem, const CeedInt* indices, const CeedScalar* r_v, CeedScalar* d_v) { 227 if (data.tidx < P1d && data.tidy < P1d) { 228 const CeedInt node = data.tidx + data.tidy*P1d; 229 const CeedInt ind = indices[node + elem * P1d*P1d]; 230 for (CeedInt comp = 0; comp < NCOMP; ++comp) 231 atomicAdd(&d_v[ind + COMPSTRIDE * comp], r_v[comp]); 232 } 233 } 234 235 //------------------------------------------------------------------------------ 236 // E-vector -> L-vector, strided 237 //------------------------------------------------------------------------------ 238 template <int NCOMP, int P1d, int STRIDES_NODE, int STRIDES_COMP, int STRIDES_ELEM> 239 inline __device__ void writeDofsStrided2d(BackendData& data, const CeedInt elem, const CeedScalar* r_v, CeedScalar* d_v) { 240 if (data.tidx < P1d && data.tidy < P1d) { 241 const CeedInt node = data.tidx + data.tidy*P1d; 242 const CeedInt ind = node * STRIDES_NODE + elem * STRIDES_ELEM; 243 for (CeedInt comp = 0; comp < NCOMP; ++comp) 244 d_v[ind + comp * STRIDES_COMP] += r_v[comp]; 245 } 246 } 247 248 //------------------------------------------------------------------------------ 249 // 2D tensor contraction x 250 //------------------------------------------------------------------------------ 251 template <int NCOMP, int P1d, int Q1d> 252 inline __device__ void ContractX2d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 253 data.slice[data.tidx+data.tidy*T1d] = *U; 254 __syncthreads(); 255 *V = 0.0; 256 if (data.tidx < Q1d && data.tidy < P1d) 257 for (CeedInt i = 0; i < P1d; ++i) 258 *V += B[i + data.tidx*P1d] * data.slice[i + data.tidy*T1d]; // Contract x direction 259 __syncthreads(); 260 } 261 262 //------------------------------------------------------------------------------ 263 // 2D tensor contract y 264 //------------------------------------------------------------------------------ 265 template <int NCOMP, int P1d, int Q1d> 266 inline __device__ void ContractY2d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 267 data.slice[data.tidx+data.tidy*T1d] = *U; 268 __syncthreads(); 269 *V = 0.0; 270 if (data.tidx < Q1d && data.tidy < Q1d) 271 for (CeedInt i = 0; i < P1d; ++i) 272 *V += B[i + data.tidy*P1d] * data.slice[data.tidx + i*T1d]; // Contract y direction 273 __syncthreads(); 274 } 275 276 //------------------------------------------------------------------------------ 277 // 2D transpose tensor contract y 278 //------------------------------------------------------------------------------ 279 template <int NCOMP, int P1d, int Q1d> 280 inline __device__ void ContractYTranspose2d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 281 data.slice[data.tidx+data.tidy*T1d] = *U; 282 __syncthreads(); 283 *V = 0.0; 284 if (data.tidx < Q1d && data.tidy < P1d) 285 for (CeedInt i = 0; i < Q1d; ++i) 286 *V += B[data.tidy + i*P1d] * data.slice[data.tidx + i*T1d]; // Contract y direction 287 __syncthreads(); 288 } 289 290 //------------------------------------------------------------------------------ 291 // 2D transpose tensor contract x 292 //------------------------------------------------------------------------------ 293 template <int NCOMP, int P1d, int Q1d> 294 inline __device__ void ContractXTranspose2d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 295 data.slice[data.tidx+data.tidy*T1d] = *U; 296 __syncthreads(); 297 *V = 0.0; 298 if (data.tidx < P1d && data.tidy < P1d) 299 for (CeedInt i = 0; i < Q1d; ++i) 300 *V += B[data.tidx + i*P1d] * data.slice[i + data.tidy*T1d]; // Contract x direction 301 __syncthreads(); 302 } 303 304 //------------------------------------------------------------------------------ 305 // 2D transpose tensor contract and add x 306 //------------------------------------------------------------------------------ 307 template <int NCOMP, int P1d, int Q1d> 308 inline __device__ void ContractXTransposeAdd2d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 309 data.slice[data.tidx+data.tidy*T1d] = *U; 310 __syncthreads(); 311 if (data.tidx < P1d && data.tidy < P1d) 312 for (CeedInt i = 0; i < Q1d; ++i) 313 *V += B[data.tidx + i*P1d] * data.slice[i + data.tidy*T1d]; // Contract x direction 314 __syncthreads(); 315 } 316 317 //------------------------------------------------------------------------------ 318 // 2D interpolate to quadrature points 319 //------------------------------------------------------------------------------ 320 template <int NCOMP, int P1d, int Q1d> 321 inline __device__ void interp2d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, CeedScalar *__restrict__ r_V) { 322 CeedScalar r_t[1]; 323 for (CeedInt comp = 0; comp < NCOMP; comp++) { 324 ContractX2d<NCOMP, P1d, Q1d>(data, r_U + comp, c_B, r_t); 325 ContractY2d<NCOMP, P1d, Q1d>(data, r_t, c_B, r_V + comp); 326 } 327 } 328 329 //------------------------------------------------------------------------------ 330 // 2D interpolate transpose 331 //------------------------------------------------------------------------------ 332 template <int NCOMP, int P1d, int Q1d> 333 inline __device__ void interpTranspose2d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, CeedScalar *__restrict__ r_V) { 334 CeedScalar r_t[1]; 335 for (CeedInt comp = 0; comp < NCOMP; comp++) { 336 ContractYTranspose2d<NCOMP, P1d, Q1d>(data, r_U + comp, c_B, r_t); 337 ContractXTranspose2d<NCOMP, P1d, Q1d>(data, r_t, c_B, r_V + comp); 338 } 339 } 340 341 //------------------------------------------------------------------------------ 342 // 2D derivatives at quadrature points 343 //------------------------------------------------------------------------------ 344 template <int NCOMP, int P1d, int Q1d> 345 inline __device__ void grad2d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 346 CeedScalar r_t[1]; 347 for (CeedInt comp = 0; comp < NCOMP; comp++) { 348 ContractX2d<NCOMP, P1d, Q1d>(data, r_U + comp, c_G, r_t); 349 ContractY2d<NCOMP, P1d, Q1d>(data, r_t, c_B, r_V + comp + 0*NCOMP); 350 ContractX2d<NCOMP, P1d, Q1d>(data, r_U + comp, c_B, r_t); 351 ContractY2d<NCOMP, P1d, Q1d>(data, r_t, c_G, r_V + comp + 1*NCOMP); 352 } 353 } 354 355 //------------------------------------------------------------------------------ 356 // 2D derivatives transpose 357 //------------------------------------------------------------------------------ 358 template <int NCOMP, int P1d, int Q1d> 359 inline __device__ void gradTranspose2d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 360 CeedScalar r_t[1]; 361 for (CeedInt comp = 0; comp < NCOMP; comp++) { 362 ContractYTranspose2d<NCOMP, P1d, Q1d>(data, r_U + comp + 0*NCOMP, c_B, r_t); 363 ContractXTranspose2d<NCOMP, P1d, Q1d>(data, r_t, c_G, r_V + comp); 364 ContractYTranspose2d<NCOMP, P1d, Q1d>(data, r_U + comp + 1*NCOMP, c_G, r_t); 365 ContractXTransposeAdd2d<NCOMP, P1d, Q1d>(data, r_t, c_B, r_V + comp); 366 } 367 } 368 369 //------------------------------------------------------------------------------ 370 // 3D 371 //------------------------------------------------------------------------------ 372 373 //------------------------------------------------------------------------------ 374 // L-vector -> E-vector, offsets provided 375 //------------------------------------------------------------------------------ 376 template <int NCOMP, int COMPSTRIDE, int P1d> 377 inline __device__ void readDofsOffset3d(BackendData& data, const CeedInt nnodes, const CeedInt elem, const CeedInt* indices, const CeedScalar* d_u, CeedScalar* r_u) { 378 if (data.tidx < P1d && data.tidy < P1d) 379 for (CeedInt z = 0; z < P1d; ++z) { 380 const CeedInt node = data.tidx + data.tidy*P1d + z*P1d*P1d; 381 const CeedInt ind = indices[node + elem * P1d*P1d*P1d]; 382 for (CeedInt comp = 0; comp < NCOMP; ++comp) 383 r_u[z+comp*P1d] = d_u[ind + COMPSTRIDE * comp]; 384 } 385 } 386 387 //------------------------------------------------------------------------------ 388 // L-vector -> E-vector, strided 389 //------------------------------------------------------------------------------ 390 template <int NCOMP, int P1d, int STRIDES_NODE, int STRIDES_COMP, int STRIDES_ELEM> 391 inline __device__ void readDofsStrided3d(BackendData& data, const CeedInt elem, const CeedScalar* d_u, CeedScalar* r_u) { 392 if (data.tidx < P1d && data.tidy < P1d) 393 for (CeedInt z = 0; z < P1d; ++z) { 394 const CeedInt node = data.tidx + data.tidy*P1d + z*P1d*P1d; 395 const CeedInt ind = node * STRIDES_NODE + elem * STRIDES_ELEM; 396 for (CeedInt comp = 0; comp < NCOMP; ++comp) 397 r_u[z+comp*P1d] = d_u[ind + comp * STRIDES_COMP]; 398 } 399 } 400 401 //------------------------------------------------------------------------------ 402 // E-vector -> Q-vector, offests provided 403 //------------------------------------------------------------------------------ 404 template <int NCOMP, int COMPSTRIDE, int Q1d> 405 inline __device__ void readSliceQuadsOffset3d(BackendData& data, const CeedInt nquads, const CeedInt elem, const CeedInt q, const CeedInt* indices, const CeedScalar* d_u, CeedScalar* r_u) { 406 if (data.tidx < Q1d && data.tidy < Q1d) { 407 const CeedInt node = data.tidx + data.tidy*Q1d + q*Q1d*Q1d; 408 const CeedInt ind = indices[node + elem * Q1d*Q1d*Q1d];; 409 for (CeedInt comp = 0; comp < NCOMP; ++comp) 410 r_u[comp] = d_u[ind + COMPSTRIDE * comp]; 411 } 412 } 413 414 //------------------------------------------------------------------------------ 415 // E-vector -> Q-vector, strided 416 //------------------------------------------------------------------------------ 417 template <int NCOMP, int Q1d, int STRIDES_NODE, int STRIDES_COMP, int STRIDES_ELEM> 418 inline __device__ void readSliceQuadsStrided3d(BackendData& data, const CeedInt elem, const CeedInt q, const CeedScalar* d_u, CeedScalar* r_u) { 419 if (data.tidx < Q1d && data.tidy < Q1d) { 420 const CeedInt node = data.tidx + data.tidy*Q1d + q*Q1d*Q1d; 421 const CeedInt ind = node * STRIDES_NODE + elem * STRIDES_ELEM; 422 for (CeedInt comp = 0; comp < NCOMP; ++comp) 423 r_u[comp] = d_u[ind + comp * STRIDES_COMP]; 424 } 425 } 426 427 //------------------------------------------------------------------------------ 428 // E-vector -> L-vector, offsets provided 429 //------------------------------------------------------------------------------ 430 template <int NCOMP, int COMPSTRIDE, int P1d> 431 inline __device__ void writeDofsOffset3d(BackendData& data, const CeedInt nnodes, const CeedInt elem, const CeedInt* indices, const CeedScalar* r_v, CeedScalar* d_v) { 432 if (data.tidx < P1d && data.tidy < P1d) 433 for (CeedInt z = 0; z < P1d; ++z) { 434 const CeedInt node = data.tidx + data.tidy*P1d + z*P1d*P1d; 435 const CeedInt ind = indices[node + elem * P1d*P1d*P1d]; 436 for (CeedInt comp = 0; comp < NCOMP; ++comp) 437 atomicAdd(&d_v[ind + COMPSTRIDE * comp], r_v[z+comp*P1d]); 438 } 439 } 440 441 //------------------------------------------------------------------------------ 442 // E-vector -> L-vector, strided 443 //------------------------------------------------------------------------------ 444 template <int NCOMP, int P1d, int STRIDES_NODE, int STRIDES_COMP, int STRIDES_ELEM> 445 inline __device__ void writeDofsStrided3d(BackendData& data, const CeedInt elem, const CeedScalar* r_v, CeedScalar* d_v) { 446 if (data.tidx < P1d && data.tidy < P1d) 447 for (CeedInt z = 0; z < P1d; ++z) { 448 const CeedInt node = data.tidx + data.tidy*P1d + z*P1d*P1d; 449 const CeedInt ind = node * STRIDES_NODE + elem * STRIDES_ELEM; 450 for (CeedInt comp = 0; comp < NCOMP; ++comp) 451 d_v[ind + comp * STRIDES_COMP] += r_v[z+comp*P1d]; 452 } 453 } 454 455 //------------------------------------------------------------------------------ 456 // 3D tensor contract x 457 //------------------------------------------------------------------------------ 458 template <int NCOMP, int P1d, int Q1d> 459 inline __device__ void ContractX3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 460 CeedScalar r_B[P1d]; 461 for (CeedInt i = 0; i < P1d; ++i) 462 r_B[i] = B[i + data.tidx*P1d]; 463 464 for (CeedInt k = 0; k < P1d; ++k) { 465 data.slice[data.tidx+data.tidy*T1d] = U[k]; 466 __syncthreads(); 467 V[k] = 0.0; 468 if (data.tidx < Q1d && data.tidy < P1d) 469 for (CeedInt i = 0; i < P1d; ++i) 470 V[k] += r_B[i] * data.slice[i + data.tidy*T1d]; // Contract x direction 471 __syncthreads(); 472 } 473 } 474 475 //------------------------------------------------------------------------------ 476 // 3D tensor contract y 477 //------------------------------------------------------------------------------ 478 template <int NCOMP, int P1d, int Q1d> 479 inline __device__ void ContractY3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 480 CeedScalar r_B[P1d]; 481 for (CeedInt i = 0; i < P1d; ++i) 482 r_B[i] = B[i + data.tidy*P1d]; 483 484 for (CeedInt k = 0; k < P1d; ++k) { 485 data.slice[data.tidx+data.tidy*T1d] = U[k]; 486 __syncthreads(); 487 V[k] = 0.0; 488 if (data.tidx < Q1d && data.tidy < Q1d) 489 for (CeedInt i = 0; i < P1d; ++i) 490 V[k] += r_B[i] * data.slice[data.tidx + i*T1d]; // Contract y direction 491 __syncthreads(); 492 } 493 } 494 495 //------------------------------------------------------------------------------ 496 // 3D tensor contract z 497 //------------------------------------------------------------------------------ 498 template <int NCOMP, int P1d, int Q1d> 499 inline __device__ void ContractZ3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 500 for (CeedInt k = 0; k < Q1d; ++k) { 501 V[k] = 0.0; 502 if (data.tidx < Q1d && data.tidy < Q1d) 503 for (CeedInt i = 0; i < P1d; ++i) 504 V[k] += B[i + k*P1d] * U[i]; // Contract z direction 505 } 506 } 507 508 //------------------------------------------------------------------------------ 509 // 3D transpose tensor contract z 510 //------------------------------------------------------------------------------ 511 template <int NCOMP, int P1d, int Q1d> 512 inline __device__ void ContractTransposeZ3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 513 for (CeedInt k = 0; k < P1d; ++k) { 514 V[k] = 0.0; 515 if (data.tidx < Q1d && data.tidy < Q1d) 516 for (CeedInt i = 0; i < Q1d; ++i) 517 V[k] += B[k + i*P1d] * U[i]; // Contract z direction 518 } 519 } 520 521 //------------------------------------------------------------------------------ 522 // 3D transpose tensor contract y 523 //------------------------------------------------------------------------------ 524 template <int NCOMP, int P1d, int Q1d> 525 inline __device__ void ContractTransposeY3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 526 CeedScalar r_B[Q1d]; 527 for (CeedInt i = 0; i < Q1d; ++i) 528 r_B[i] = B[data.tidy + i*P1d]; 529 530 for (CeedInt k = 0; k < P1d; ++k) { 531 data.slice[data.tidx+data.tidy*T1d] = U[k]; 532 __syncthreads(); 533 V[k] = 0.0; 534 if (data.tidx < Q1d && data.tidy < P1d) 535 for (CeedInt i = 0; i < Q1d; ++i) 536 V[k] += r_B[i] * data.slice[data.tidx + i*T1d]; // Contract y direction 537 __syncthreads(); 538 } 539 } 540 541 //------------------------------------------------------------------------------ 542 // 3D transpose tensor contract add y 543 //------------------------------------------------------------------------------ 544 template <int NCOMP, int P1d, int Q1d> 545 inline __device__ void ContractTransposeAddY3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 546 CeedScalar r_B[Q1d]; 547 for (CeedInt i = 0; i < Q1d; ++i) 548 r_B[i] = B[data.tidy + i*P1d]; 549 550 for (CeedInt k = 0; k < P1d; ++k) { 551 data.slice[data.tidx+data.tidy*T1d] = U[k]; 552 __syncthreads(); 553 if (data.tidx < Q1d && data.tidy < P1d) 554 for (CeedInt i = 0; i < Q1d; ++i) 555 V[k] += r_B[i] * data.slice[data.tidx + i*T1d]; // Contract y direction 556 __syncthreads(); 557 } 558 } 559 560 //------------------------------------------------------------------------------ 561 // 3D transpose tensor contract x 562 //------------------------------------------------------------------------------ 563 template <int NCOMP, int P1d, int Q1d> 564 inline __device__ void ContractTransposeX3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 565 CeedScalar r_B[Q1d]; 566 for (CeedInt i = 0; i < Q1d; ++i) 567 r_B[i] = B[data.tidx + i*P1d]; 568 569 for (CeedInt k = 0; k < P1d; ++k) { 570 data.slice[data.tidx+data.tidy*T1d] = U[k]; 571 __syncthreads(); 572 V[k] = 0.0; 573 if (data.tidx < P1d && data.tidy < P1d) 574 for (CeedInt i = 0; i < Q1d; ++i) 575 V[k] += r_B[i] * data.slice[i + data.tidy*T1d]; // Contract x direction 576 __syncthreads(); 577 } 578 } 579 580 //------------------------------------------------------------------------------ 581 // 3D transpose tensor contract add x 582 //------------------------------------------------------------------------------ 583 template <int NCOMP, int P1d, int Q1d> 584 inline __device__ void ContractTransposeAddX3d(BackendData& data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) { 585 CeedScalar r_B[Q1d]; 586 for (CeedInt i = 0; i < Q1d; ++i) 587 r_B[i] = B[data.tidx + i*P1d]; 588 589 for (CeedInt k = 0; k < P1d; ++k) { 590 data.slice[data.tidx+data.tidy*T1d] = U[k]; 591 __syncthreads(); 592 if (data.tidx < P1d && data.tidy < P1d) 593 for (CeedInt i = 0; i < Q1d; ++i) 594 V[k] += r_B[i] * data.slice[i + data.tidy*T1d]; // Contract x direction 595 __syncthreads(); 596 } 597 } 598 599 //------------------------------------------------------------------------------ 600 // 3D interpolate to quadrature points 601 //------------------------------------------------------------------------------ 602 template <int NCOMP, int P1d, int Q1d> 603 inline __device__ void interp3d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, CeedScalar *__restrict__ r_V) { 604 CeedScalar r_t1[T1d]; 605 CeedScalar r_t2[T1d]; 606 for (CeedInt comp = 0; comp < NCOMP; comp++) { 607 ContractX3d<NCOMP, P1d, Q1d>(data, r_U + comp*P1d, c_B, r_t1); 608 ContractY3d<NCOMP, P1d, Q1d>(data, r_t1, c_B, r_t2); 609 ContractZ3d<NCOMP, P1d, Q1d>(data, r_t2, c_B, r_V + comp*Q1d); 610 } 611 } 612 613 //------------------------------------------------------------------------------ 614 // 3D interpolate transpose 615 //------------------------------------------------------------------------------ 616 template <int NCOMP, int P1d, int Q1d> 617 inline __device__ void interpTranspose3d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, CeedScalar *__restrict__ r_V) { 618 CeedScalar r_t1[T1d]; 619 CeedScalar r_t2[T1d]; 620 for (CeedInt comp = 0; comp < NCOMP; comp++) { 621 ContractTransposeZ3d<NCOMP, P1d, Q1d>(data, r_U + comp*Q1d, c_B, r_t1); 622 ContractTransposeY3d<NCOMP, P1d, Q1d>(data, r_t1, c_B, r_t2); 623 ContractTransposeX3d<NCOMP, P1d, Q1d>(data, r_t2, c_B, r_V + comp*P1d); 624 } 625 } 626 627 //------------------------------------------------------------------------------ 628 // 3D derivatives at quadrature points 629 //------------------------------------------------------------------------------ 630 template <int NCOMP, int P1d, int Q1d> 631 inline __device__ void grad3d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 632 CeedScalar r_t1[T1d]; 633 CeedScalar r_t2[T1d]; 634 for (CeedInt comp = 0; comp < NCOMP; comp++) { 635 ContractX3d<NCOMP, P1d, Q1d>(data, r_U + comp*P1d, c_G, r_t1); 636 ContractY3d<NCOMP, P1d, Q1d>(data, r_t1, c_B, r_t2); 637 ContractZ3d<NCOMP, P1d, Q1d>(data, r_t2, c_B, r_V + comp*Q1d + 0*NCOMP*Q1d); 638 ContractX3d<NCOMP, P1d, Q1d>(data, r_U + comp*P1d, c_B, r_t1); 639 ContractY3d<NCOMP, P1d, Q1d>(data, r_t1, c_G, r_t2); 640 ContractZ3d<NCOMP, P1d, Q1d>(data, r_t2, c_B, r_V + comp*Q1d + 1*NCOMP*Q1d); 641 ContractX3d<NCOMP, P1d, Q1d>(data, r_U + comp*P1d, c_B, r_t1); 642 ContractY3d<NCOMP, P1d, Q1d>(data, r_t1, c_B, r_t2); 643 ContractZ3d<NCOMP, P1d, Q1d>(data, r_t2, c_G, r_V + comp*Q1d + 2*NCOMP*Q1d); 644 } 645 } 646 647 //------------------------------------------------------------------------------ 648 // 3D derivatives transpose 649 //------------------------------------------------------------------------------ 650 template <int NCOMP, int P1d, int Q1d> 651 inline __device__ void gradTranspose3d(BackendData& data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 652 CeedScalar r_t1[T1d]; 653 CeedScalar r_t2[T1d]; 654 for (CeedInt comp = 0; comp < NCOMP; comp++) { 655 ContractTransposeZ3d<NCOMP, P1d, Q1d>(data, r_U + comp*Q1d + 0*NCOMP*Q1d, c_B, r_t1); 656 ContractTransposeY3d<NCOMP, P1d, Q1d>(data, r_t1, c_B, r_t2); 657 ContractTransposeX3d<NCOMP, P1d, Q1d>(data, r_t2, c_G, r_V + comp*P1d); 658 ContractTransposeZ3d<NCOMP, P1d, Q1d>(data, r_U + comp*Q1d + 1*NCOMP*Q1d, c_B, r_t1); 659 ContractTransposeY3d<NCOMP, P1d, Q1d>(data, r_t1, c_G, r_t2); 660 ContractTransposeAddX3d<NCOMP,P1d, Q1d>(data, r_t2, c_B, r_V + comp*P1d); 661 ContractTransposeZ3d<NCOMP, P1d, Q1d>(data, r_U + comp*Q1d + 2*NCOMP*Q1d, c_G, r_t1); 662 ContractTransposeY3d<NCOMP, P1d, Q1d>(data, r_t1, c_B, r_t2); 663 ContractTransposeAddX3d<NCOMP, P1d, Q1d>(data, r_t2, c_B, r_V + comp*P1d); 664 } 665 } 666 667 //------------------------------------------------------------------------------ 668 // 3D collocated derivatives computation 669 //------------------------------------------------------------------------------ 670 template <int NCOMP, int Q1d> 671 inline __device__ void gradCollo3d(BackendData& data, const CeedInt q, const CeedScalar *__restrict__ r_U, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 672 if (data.tidx < Q1d && data.tidy < Q1d) { 673 for (CeedInt comp = 0; comp < NCOMP; ++comp) { 674 data.slice[data.tidx + data.tidy*T1d] = r_U[q + comp*Q1d]; 675 __syncthreads(); 676 // X derivative 677 r_V[comp+0*NCOMP] = 0.0; 678 for (CeedInt i = 0; i < Q1d; ++i) 679 r_V[comp+0*NCOMP] += c_G[i + data.tidx*Q1d] * data.slice[i + data.tidy*T1d]; // Contract x direction (X derivative) 680 // Y derivative 681 r_V[comp+1*NCOMP] = 0.0; 682 for (CeedInt i = 0; i < Q1d; ++i) 683 r_V[comp+1*NCOMP] += c_G[i + data.tidy*Q1d] * data.slice[data.tidx + i*T1d]; // Contract y direction (Y derivative) 684 // Z derivative 685 r_V[comp+2*NCOMP] = 0.0; 686 for (CeedInt i = 0; i < Q1d; ++i) 687 r_V[comp+2*NCOMP] += c_G[i + q*Q1d] * r_U[i + comp*Q1d]; // Contract z direction (Z derivative) 688 __syncthreads(); 689 } 690 } 691 } 692 693 //------------------------------------------------------------------------------ 694 // 3D collocated derivatives transpose 695 //------------------------------------------------------------------------------ 696 template <int NCOMP, int Q1d> 697 inline __device__ void gradColloTranspose3d(BackendData& data, const CeedInt q, const CeedScalar *__restrict__ r_U, const CeedScalar *c_G, CeedScalar *__restrict__ r_V) { 698 if (data.tidx < Q1d && data.tidy < Q1d) { 699 for (CeedInt comp = 0; comp < NCOMP; ++comp) { 700 // X derivative 701 data.slice[data.tidx + data.tidy*T1d] = r_U[comp + 0*NCOMP]; 702 __syncthreads(); 703 for (CeedInt i = 0; i < Q1d; ++i) 704 r_V[q+comp*Q1d] += c_G[data.tidx + i*Q1d] * data.slice[i + data.tidy*T1d]; // Contract x direction (X derivative) 705 __syncthreads(); 706 // Y derivative 707 data.slice[data.tidx + data.tidy*T1d] = r_U[comp + 1*NCOMP]; 708 __syncthreads(); 709 for (CeedInt i = 0; i < Q1d; ++i) 710 r_V[q+comp*Q1d] += c_G[data.tidy + i*Q1d] * data.slice[data.tidx + i*T1d]; // Contract y direction (Y derivative) 711 __syncthreads(); 712 // Z derivative 713 for (CeedInt i = 0; i < Q1d; ++i) 714 r_V[i+comp*Q1d] += c_G[i + q*Q1d] * r_U[comp + 2*NCOMP]; // PARTIAL contract z direction (Z derivative) 715 } 716 } 717 } 718 719 //------------------------------------------------------------------------------ 720 // 1D quadrature weights 721 //------------------------------------------------------------------------------ 722 template <int Q1d> 723 inline __device__ void weight1d(BackendData& data, const CeedScalar *qweight1d, CeedScalar *w) { 724 *w = (data.tidx < Q1d) ? qweight1d[data.tidx] : 0.0; 725 } 726 727 //------------------------------------------------------------------------------ 728 // 2D quadrature weights 729 //------------------------------------------------------------------------------ 730 template <int Q1d> 731 inline __device__ void weight2d(BackendData& data, const CeedScalar *qweight1d, CeedScalar *w) { 732 *w = (data.tidx < Q1d && data.tidy < Q1d) ? 733 qweight1d[data.tidx]*qweight1d[data.tidy] : 0.0; 734 } 735 736 //------------------------------------------------------------------------------ 737 // 3D quadrature weights 738 //------------------------------------------------------------------------------ 739 template <int Q1d> 740 inline __device__ void weight3d(BackendData& data, const CeedScalar *qweight1d, CeedScalar *w) { 741 const bool quad = (data.tidx < Q1d && data.tidy < Q1d); 742 const CeedScalar pw = quad ? qweight1d[data.tidx]*qweight1d[data.tidy] : 0.0; 743 for (CeedInt z = 0; z < Q1d; ++z) 744 w[z] = quad ? pw*qweight1d[z] : 0.0; 745 } 746 747 ); 748 //------------------------------------------------------------------------------ 749 // Build singe operator kernel 750 //------------------------------------------------------------------------------ 751 extern "C" int CeedCudaGenOperatorBuild(CeedOperator op) { 752 753 using std::ostringstream; 754 using std::string; 755 int ierr; 756 bool setupdone; 757 ierr = CeedOperatorIsSetupDone(op, &setupdone); CeedChkBackend(ierr); 758 if (setupdone) return CEED_ERROR_SUCCESS; 759 Ceed ceed; 760 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 761 CeedOperator_Cuda_gen *data; 762 ierr = CeedOperatorGetData(op, &data); CeedChkBackend(ierr); 763 CeedQFunction qf; 764 CeedQFunction_Cuda_gen *qf_data; 765 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 766 ierr = CeedQFunctionGetData(qf, &qf_data); CeedChkBackend(ierr); 767 CeedInt Q, P1d, Q1d = 0, numelements, elemsize, numinputfields, 768 numoutputfields, ncomp, dim = 0, lsize; 769 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 770 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 771 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 772 CeedChkBackend(ierr); 773 CeedOperatorField *opinputfields, *opoutputfields; 774 ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields); 775 CeedChkBackend(ierr); 776 CeedQFunctionField *qfinputfields, *qfoutputfields; 777 ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields); 778 CeedChkBackend(ierr); 779 CeedEvalMode emode; 780 CeedBasis basis; 781 CeedBasis_Cuda_shared *basis_data; 782 CeedElemRestriction Erestrict; 783 CeedElemRestriction_Cuda *restr_data; 784 785 ostringstream code; 786 string devFunctions(deviceFunctions); 787 788 // Add atomicAdd function for old NVidia architectures 789 struct cudaDeviceProp prop; 790 Ceed_Cuda *ceed_data; 791 ierr = CeedGetData(ceed, &ceed_data); CeedChkBackend(ierr); 792 ierr = cudaGetDeviceProperties(&prop, ceed_data->deviceId); 793 if (prop.major<6){ 794 code << atomicAdd; 795 } 796 797 code << devFunctions; 798 799 string qFunction(qf_data->qFunctionSource); 800 string qFunctionName(qf_data->qFunctionName); 801 string oper; 802 oper = "CeedKernel_Cuda_gen_" + qFunctionName; 803 804 code << "\n#define CEED_QFUNCTION(name) inline __device__ int name\n"; 805 code << "#define CeedPragmaSIMD\n"; 806 code << "#define CEED_ERROR_SUCCESS 0\n\n"; 807 808 // Find dim and Q1d 809 bool useCollograd = true; 810 data->maxP1d = 0; 811 for (CeedInt i = 0; i < numinputfields; i++) { 812 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChkBackend(ierr); 813 if (basis != CEED_BASIS_COLLOCATED) { 814 ierr = CeedBasisGetData(basis, &basis_data); CeedChkBackend(ierr); 815 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 816 CeedChkBackend(ierr); 817 818 // Check for collocated gradient 819 useCollograd = useCollograd && basis_data->d_collograd1d; 820 821 // Collect dim and Q1d 822 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 823 bool isTensor; 824 ierr = CeedBasisIsTensor(basis, &isTensor); CeedChkBackend(ierr); 825 if (isTensor) { 826 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChkBackend(ierr); 827 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChkBackend(ierr); 828 if (P1d>data->maxP1d) data->maxP1d = P1d; 829 } else { 830 // LCOV_EXCL_START 831 return CeedError(ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); 832 // LCOV_EXCL_STOP 833 } 834 } 835 } 836 // Check output bases for Q1d, dim as well 837 // The only imput basis might be CEED_BASIS_COLLOCATED 838 for (CeedInt i = 0; i < numoutputfields; i++) { 839 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); CeedChkBackend(ierr); 840 841 if (basis != CEED_BASIS_COLLOCATED) { 842 ierr = CeedBasisGetData(basis, &basis_data); CeedChkBackend(ierr); 843 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 844 CeedChkBackend(ierr); 845 846 // Collect dim and Q1d 847 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 848 bool isTensor; 849 ierr = CeedBasisIsTensor(basis, &isTensor); CeedChkBackend(ierr); 850 if (isTensor) { 851 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChkBackend(ierr); 852 } else { 853 // LCOV_EXCL_START 854 return CeedError(ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); 855 // LCOV_EXCL_STOP 856 } 857 858 // Check for collocated gradient 859 useCollograd = useCollograd && basis_data->d_collograd1d; 860 } 861 } 862 data->dim = dim; 863 data->Q1d = Q1d; 864 865 // Define CEED_Q_VLA 866 if (dim != 3 || useCollograd) { 867 code << "\n#define CEED_Q_VLA 1\n\n"; 868 } else { 869 code << "\n#define CEED_Q_VLA "<<Q1d<<"\n\n"; 870 } 871 872 code << qFunction; 873 874 // Setup 875 code << "\n// -----------------------------------------------------------------------------\n"; 876 code << "\nextern \"C\" __global__ void "<<oper<<"(CeedInt nelem, void* ctx, CudaFieldsInt indices, CudaFields fields, CudaFields B, CudaFields G, CeedScalar* W) {\n"; 877 for (CeedInt i = 0; i < numinputfields; i++) { 878 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 879 CeedChkBackend(ierr); 880 if (emode != CEED_EVAL_WEIGHT) { // Skip CEED_EVAL_WEIGHT 881 code << " const CeedScalar* d_u" <<i<<" = fields.in["<<i<<"];\n"; 882 } 883 } 884 885 for (CeedInt i = 0; i < numoutputfields; i++) { 886 code << " CeedScalar* d_v"<<i<<" = fields.out["<<i<<"];\n"; 887 } 888 889 code << " const CeedInt Dim = "<<dim<<";\n"; 890 code << " const CeedInt Q1d = "<<Q1d<<";\n"; 891 892 code << " extern __shared__ CeedScalar slice[];\n"; 893 code << " BackendData data;\n"; 894 code << " data.tidx = threadIdx.x;\n"; 895 code << " data.tidy = threadIdx.y;\n"; 896 code << " data.tidz = threadIdx.z;\n"; 897 code << " data.tid = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; 898 code << " data.slice = slice+data.tidz*T1d"<<(dim>1?"*T1d":"")<<";\n"; 899 900 code << "\n // -- Input field constants and basis data --\n"; 901 //Initialize constants, and matrices B and G 902 for (CeedInt i = 0; i < numinputfields; i++) { 903 code << " // ---- Input field "<<i<<" ----\n"; 904 // Get elemsize, emode, ncomp 905 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 906 CeedChkBackend(ierr); 907 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 908 CeedChkBackend(ierr); 909 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 910 CeedChkBackend(ierr); 911 ierr = CeedElemRestrictionGetNumComponents(Erestrict, &ncomp); 912 CeedChkBackend(ierr); 913 914 // Set field constants 915 if (emode != CEED_EVAL_WEIGHT) { 916 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChkBackend(ierr); 917 if (basis != CEED_BASIS_COLLOCATED) { 918 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChkBackend(ierr); 919 code << " const CeedInt P_in_"<<i<<" = "<<P1d<<";\n"; 920 } else { 921 code << " const CeedInt P_in_"<<i<<" = "<<Q1d<<";\n"; 922 } 923 code << " const CeedInt ncomp_in_"<<i<<" = "<<ncomp<<";\n"; 924 } 925 926 // Load basis data 927 code << " // EvalMode: "<<CeedEvalModes[emode]<<"\n"; 928 switch (emode) { 929 case CEED_EVAL_NONE: 930 break; 931 case CEED_EVAL_INTERP: 932 ierr = CeedBasisGetData(basis, &basis_data); CeedChkBackend(ierr); 933 data->B.in[i] = basis_data->d_interp1d; 934 code << " __shared__ double s_B_in_"<<i<<"["<<P1d*Q1d<<"];\n"; 935 code << " loadMatrix<P_in_"<<i<<",Q1d>(data, B.in["<<i<<"], s_B_in_"<<i<<");\n"; 936 break; 937 case CEED_EVAL_GRAD: 938 ierr = CeedBasisGetData(basis, &basis_data); CeedChkBackend(ierr); 939 data->B.in[i] = basis_data->d_interp1d; 940 code << " __shared__ double s_B_in_"<<i<<"["<<P1d*Q1d<<"];\n"; 941 code << " loadMatrix<P_in_"<<i<<",Q1d>(data, B.in["<<i<<"], s_B_in_"<<i<<");\n"; 942 if (useCollograd) { 943 data->G.in[i] = basis_data->d_collograd1d; 944 code << " __shared__ double s_G_in_"<<i<<"["<<Q1d*Q1d<<"];\n"; 945 code << " loadMatrix<Q1d,Q1d>(data, G.in["<<i<<"], s_G_in_"<<i<<");\n"; 946 } else { 947 data->G.in[i] = basis_data->d_grad1d; 948 code << " __shared__ double s_G_in_"<<i<<"["<<P1d*Q1d<<"];\n"; 949 code << " loadMatrix<P_in_"<<i<<",Q1d>(data, G.in["<<i<<"], s_G_in_"<<i<<");\n"; 950 } 951 break; 952 case CEED_EVAL_WEIGHT: 953 break; // No action 954 case CEED_EVAL_DIV: 955 break; // TODO: Not implemented 956 case CEED_EVAL_CURL: 957 break; // TODO: Not implemented 958 } 959 } 960 961 code << "\n // -- Output field constants and basis data --\n"; 962 for (CeedInt i = 0; i < numoutputfields; i++) { 963 code << " // ---- Output field "<<i<<" ----\n"; 964 // Get elemsize, emode, ncomp 965 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 966 CeedChkBackend(ierr); 967 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 968 CeedChkBackend(ierr); 969 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 970 CeedChkBackend(ierr); 971 ierr = CeedElemRestrictionGetNumComponents(Erestrict, &ncomp); 972 CeedChkBackend(ierr); 973 974 // Set field constants 975 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); CeedChkBackend(ierr); 976 if (basis != CEED_BASIS_COLLOCATED) { 977 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChkBackend(ierr); 978 code << " const CeedInt P_out_"<<i<<" = "<<P1d<<";\n"; 979 } else { 980 code << " const CeedInt P_out_"<<i<<" = "<<Q1d<<";\n"; 981 } 982 code << " const CeedInt ncomp_out_"<<i<<" = "<<ncomp<<";\n"; 983 984 // Load basis data 985 code << " // EvalMode: "<<CeedEvalModes[emode]<<"\n"; 986 switch (emode) { 987 case CEED_EVAL_NONE: 988 break; // No action 989 case CEED_EVAL_INTERP: 990 ierr = CeedBasisGetData(basis, &basis_data); CeedChkBackend(ierr); 991 data->B.out[i] = basis_data->d_interp1d; 992 code << " __shared__ double s_B_out_"<<i<<"["<<P1d*Q1d<<"];\n"; 993 code << " loadMatrix<P_out_"<<i<<",Q1d>(data, B.out["<<i<<"], s_B_out_"<<i<<");\n"; 994 break; 995 case CEED_EVAL_GRAD: 996 ierr = CeedBasisGetData(basis, &basis_data); CeedChkBackend(ierr); 997 data->B.out[i] = basis_data->d_interp1d; 998 code << " __shared__ double s_B_out_"<<i<<"["<<P1d*Q1d<<"];\n"; 999 code << " loadMatrix<P_out_"<<i<<",Q1d>(data, B.out["<<i<<"], s_B_out_"<<i<<");\n"; 1000 if (useCollograd) { 1001 data->G.out[i] = basis_data->d_collograd1d; 1002 code << " __shared__ double s_G_out_"<<i<<"["<<Q1d*Q1d<<"];\n"; 1003 code << " loadMatrix<Q1d,Q1d>(data, G.out["<<i<<"], s_G_out_"<<i<<");\n"; 1004 } else { 1005 data->G.out[i] = basis_data->d_grad1d; 1006 code << " __shared__ double s_G_out_"<<i<<"["<<P1d*Q1d<<"];\n"; 1007 code << " loadMatrix<P_out_"<<i<<",Q1d>(data, G.out["<<i<<"], s_G_out_"<<i<<");\n"; 1008 } 1009 break; 1010 // LCOV_EXCL_START 1011 case CEED_EVAL_WEIGHT: { 1012 Ceed ceed; 1013 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1014 return CeedError(ceed, CEED_ERROR_BACKEND, 1015 "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1016 break; // Should not occur 1017 } 1018 case CEED_EVAL_DIV: 1019 break; // TODO: Not implemented 1020 case CEED_EVAL_CURL: 1021 break; // TODO: Not implemented 1022 // LCOV_EXCL_STOP 1023 } 1024 } 1025 code << "\n // -- Element loop --\n"; 1026 code << " __syncthreads();\n"; 1027 code << " for (CeedInt elem = blockIdx.x*blockDim.z + threadIdx.z; elem < nelem; elem += gridDim.x*blockDim.z) {\n"; 1028 // Input basis apply if needed 1029 // Generate the correct eval mode code for each input 1030 code << " // -- Input field restrictions and basis actions --\n"; 1031 for (CeedInt i = 0; i < numinputfields; i++) { 1032 code << " // ---- Input field "<<i<<" ----\n"; 1033 // Get elemsize, emode, ncomp 1034 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 1035 CeedChkBackend(ierr); 1036 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 1037 CeedChkBackend(ierr); 1038 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 1039 CeedChkBackend(ierr); 1040 ierr = CeedElemRestrictionGetNumComponents(Erestrict, &ncomp); 1041 CeedChkBackend(ierr); 1042 1043 // Restriction 1044 if (emode != CEED_EVAL_WEIGHT && 1045 !((emode == CEED_EVAL_NONE) && useCollograd)) { 1046 code << " CeedScalar r_u"<<i<<"[ncomp_in_"<<i<<"*P_in_"<<i<<"];\n"; 1047 1048 bool isStrided; 1049 ierr = CeedElemRestrictionIsStrided(Erestrict, &isStrided); CeedChkBackend(ierr); 1050 if (!isStrided) { 1051 ierr = CeedElemRestrictionGetLVectorSize(Erestrict, &lsize); 1052 CeedChkBackend(ierr); 1053 code << " const CeedInt lsize_in_"<<i<<" = "<<lsize<<";\n"; 1054 CeedInt compstride; 1055 ierr = CeedElemRestrictionGetCompStride(Erestrict, &compstride); CeedChkBackend(ierr); 1056 code << " // CompStride: "<<compstride<<"\n"; 1057 ierr = CeedElemRestrictionGetData(Erestrict, &restr_data); CeedChkBackend(ierr); 1058 data->indices.in[i] = restr_data->d_ind; 1059 code << " readDofsOffset"<<dim<<"d<ncomp_in_"<<i<<", "<<compstride<<", P_in_"<<i<<">(data, lsize_in_"<<i<<", elem, indices.in["<<i<<"], d_u"<<i<<", r_u"<<i<<");\n"; 1060 } else { 1061 bool backendstrides; 1062 ierr = CeedElemRestrictionHasBackendStrides(Erestrict, &backendstrides); 1063 CeedChkBackend(ierr); 1064 CeedInt nelem; 1065 ierr = CeedElemRestrictionGetNumElements(Erestrict, &nelem); 1066 CeedChkBackend(ierr); 1067 CeedInt strides[3] = {1, elemsize*nelem, elemsize}; 1068 if (!backendstrides) { 1069 ierr = CeedElemRestrictionGetStrides(Erestrict, &strides); 1070 CeedChkBackend(ierr); 1071 } 1072 code << " // Strides: {"<<strides[0]<<", "<<strides[1]<<", "<<strides[2]<<"}\n"; 1073 code << " readDofsStrided"<<dim<<"d<ncomp_in_"<<i<<",P_in_"<<i<<","<<strides[0]<<","<<strides[1]<<","<<strides[2]<<">(data, elem, d_u"<<i<<", r_u"<<i<<");\n"; 1074 } 1075 } 1076 1077 // Basis action 1078 code << " // EvalMode: "<<CeedEvalModes[emode]<<"\n"; 1079 switch (emode) { 1080 case CEED_EVAL_NONE: 1081 if (!useCollograd) { 1082 code << " CeedScalar* r_t"<<i<<" = r_u"<<i<<";\n"; 1083 } 1084 break; 1085 case CEED_EVAL_INTERP: 1086 code << " CeedScalar r_t"<<i<<"[ncomp_in_"<<i<<"*Q1d];\n"; 1087 code << " interp"<<dim<<"d<ncomp_in_"<<i<<",P_in_"<<i<<",Q1d>(data, r_u"<<i<<", s_B_in_"<<i<<", r_t"<<i<<");\n"; 1088 break; 1089 case CEED_EVAL_GRAD: 1090 if (useCollograd) { 1091 code << " CeedScalar r_t"<<i<<"[ncomp_in_"<<i<<"*Q1d];\n"; 1092 code << " interp"<<dim<<"d<ncomp_in_"<<i<<",P_in_"<<i<<",Q1d>(data, r_u"<<i<<", s_B_in_"<<i<<", r_t"<<i<<");\n"; 1093 } else { 1094 code << " CeedScalar r_t"<<i<<"[ncomp_in_"<<i<<"*Dim*Q1d];\n"; 1095 code << " grad"<<dim<<"d<ncomp_in_"<<i<<",P_in_"<<i<<",Q1d>(data, r_u"<<i<<", s_B_in_"<<i<<", s_G_in_"<<i<<", r_t"<<i<<");\n"; 1096 } 1097 break; 1098 case CEED_EVAL_WEIGHT: 1099 code << " CeedScalar r_t"<<i<<"[Q1d];\n"; 1100 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); CeedChkBackend(ierr); 1101 ierr = CeedBasisGetData(basis, &basis_data); CeedChkBackend(ierr); 1102 data->W = basis_data->d_qweight1d; 1103 code << " weight"<<dim<<"d<Q1d>(data, W, r_t"<<i<<");\n"; 1104 break; // No action 1105 case CEED_EVAL_DIV: 1106 break; // TODO: Not implemented 1107 case CEED_EVAL_CURL: 1108 break; // TODO: Not implemented 1109 } 1110 } 1111 1112 // Q function 1113 code << "\n // -- Output field setup --\n"; 1114 for (CeedInt i = 0; i < numoutputfields; i++) { 1115 code << "\n // ---- Output field "<<i<<" ----\n"; 1116 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 1117 CeedChkBackend(ierr); 1118 if (emode==CEED_EVAL_GRAD) 1119 { 1120 if (useCollograd) { 1121 //Accumulator for gradient slices 1122 code << " CeedScalar r_tt"<<i<<"[ncomp_out_"<<i<<"*Q1d];\n"; 1123 code << " for (CeedInt i = 0; i < ncomp_out_"<<i<<"; ++i) {\n"; 1124 code << " for (CeedInt j = 0; j < Q1d; ++j) {\n"; 1125 code << " r_tt"<<i<<"[j + i*Q1d] = 0.0;\n"; 1126 code << " }\n"; 1127 code << " }\n"; 1128 } else { 1129 code << " CeedScalar r_tt"<<i<<"[ncomp_out_"<<i<<"*Dim*Q1d];\n"; 1130 } 1131 } 1132 if (emode==CEED_EVAL_NONE || emode==CEED_EVAL_INTERP) 1133 { 1134 code << " CeedScalar r_tt"<<i<<"[ncomp_out_"<<i<<"*Q1d];\n"; 1135 } 1136 } 1137 // We treat quadrature points per slice in 3d to save registers 1138 if (useCollograd) { 1139 code << "\n // Note: Collocated Gradient\n"; 1140 code << "#pragma unroll\n"; 1141 code << " for (CeedInt q=0; q<Q1d; q++) {\n"; 1142 code << " // -- Input fields --\n"; 1143 for (CeedInt i = 0; i < numinputfields; i++) { 1144 code << " // ---- Input field "<<i<<" ----\n"; 1145 // Get elemsize, emode, ncomp 1146 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 1147 CeedChkBackend(ierr); 1148 // Basis action 1149 code << " // EvalMode: "<<CeedEvalModes[emode]<<"\n"; 1150 switch (emode) { 1151 case CEED_EVAL_NONE: 1152 code << " CeedScalar r_q"<<i<<"[ncomp_in_"<<i<<"];\n"; 1153 1154 bool isStrided; 1155 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); CeedChkBackend(ierr); 1156 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); CeedChkBackend(ierr); 1157 ierr = CeedElemRestrictionIsStrided(Erestrict, &isStrided); CeedChkBackend(ierr); 1158 if (!isStrided) { 1159 ierr = CeedElemRestrictionGetLVectorSize(Erestrict, &lsize); 1160 CeedChkBackend(ierr); 1161 code << " const CeedInt lsize_in_"<<i<<" = "<<lsize<<";\n"; 1162 CeedInt compstride; 1163 ierr = CeedElemRestrictionGetCompStride(Erestrict, &compstride); CeedChkBackend(ierr); 1164 code << " // CompStride: "<<compstride<<"\n"; 1165 ierr = CeedElemRestrictionGetData(Erestrict, &restr_data); CeedChkBackend(ierr); 1166 data->indices.in[i] = restr_data->d_ind; 1167 code << " readSliceQuadsOffset"<<"3d<ncomp_in_"<<i<<", "<<compstride<<", Q1d>(data, lsize_in_"<<i<<", elem, q, indices.in["<<i<<"], d_u"<<i<<", r_q"<<i<<");\n"; 1168 } else { 1169 bool backendstrides; 1170 ierr = CeedElemRestrictionHasBackendStrides(Erestrict, &backendstrides); 1171 CeedChkBackend(ierr); 1172 CeedInt nelem; 1173 ierr = CeedElemRestrictionGetNumElements(Erestrict, &nelem); 1174 CeedChkBackend(ierr); 1175 CeedInt strides[3] = {1, elemsize*nelem, elemsize}; 1176 if (!backendstrides) { 1177 ierr = CeedElemRestrictionGetStrides(Erestrict, &strides); 1178 CeedChkBackend(ierr); 1179 } 1180 code << " // Strides: {"<<strides[0]<<", "<<strides[1]<<", "<<strides[2]<<"}\n"; 1181 code << " readSliceQuadsStrided"<<"3d<ncomp_in_"<<i<<",Q1d"","<<strides[0]<<","<<strides[1]<<","<<strides[2]<<">(data, elem, q, d_u"<<i<<", r_q"<<i<<");\n"; 1182 } 1183 break; 1184 case CEED_EVAL_INTERP: 1185 code << " CeedScalar r_q"<<i<<"[ncomp_in_"<<i<<"];\n"; 1186 code << " for (CeedInt j = 0; j < ncomp_in_"<<i<<" ; ++j) {\n"; 1187 code << " r_q"<<i<<"[j] = r_t"<<i<<"[q + j*Q1d];\n"; 1188 code << " }\n"; 1189 break; 1190 case CEED_EVAL_GRAD: 1191 code << " CeedScalar r_q"<<i<<"[ncomp_in_"<<i<<"*Dim];\n"; 1192 code << " gradCollo3d<ncomp_in_"<<i<<",Q1d>(data, q, r_t"<<i<<", s_G_in_"<<i<<", r_q"<<i<<");\n"; 1193 break; 1194 case CEED_EVAL_WEIGHT: 1195 code << " CeedScalar r_q"<<i<<"[1];\n"; 1196 code << " r_q"<<i<<"[0] = r_t"<<i<<"[q];\n"; 1197 break; // No action 1198 case CEED_EVAL_DIV: 1199 break; // TODO: Not implemented 1200 case CEED_EVAL_CURL: 1201 break; // TODO: Not implemented 1202 } 1203 } 1204 code << "\n // -- Output fields --\n"; 1205 for (CeedInt i = 0; i < numoutputfields; i++) { 1206 code << " // ---- Output field "<<i<<" ----\n"; 1207 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 1208 CeedChkBackend(ierr); 1209 // Basis action 1210 switch (emode) { 1211 case CEED_EVAL_NONE: 1212 code << " CeedScalar r_qq"<<i<<"[ncomp_out_"<<i<<"];\n"; 1213 break; // No action 1214 case CEED_EVAL_INTERP: 1215 code << " CeedScalar r_qq"<<i<<"[ncomp_out_"<<i<<"];\n"; 1216 break; 1217 case CEED_EVAL_GRAD: 1218 code << " CeedScalar r_qq"<<i<<"[ncomp_out_"<<i<<"*Dim];\n"; 1219 break; 1220 case CEED_EVAL_WEIGHT: 1221 break; // Should not occur 1222 case CEED_EVAL_DIV: 1223 break; // TODO: Not implemented 1224 case CEED_EVAL_CURL: 1225 break; // TODO: Not implemented 1226 } 1227 } 1228 } else { 1229 code << "\n // Note: No Collocated Gradient\n"; 1230 code << " // -- Input fields --\n"; 1231 for (CeedInt i = 0; i < numinputfields; i++) { 1232 code << " // ---- Input field "<<i<<" ----\n"; 1233 code << " CeedScalar* r_q"<<i<<" = r_t"<<i<<";\n"; 1234 } 1235 code << " // -- Output fields --\n"; 1236 for (CeedInt i = 0; i < numoutputfields; i++) { 1237 code << " // ---- Output field "<<i<<" ----\n"; 1238 code << " CeedScalar* r_qq"<<i<<" = r_tt"<<i<<";\n"; 1239 } 1240 } 1241 code << "\n // -- QFunction Inputs and outputs --\n"; 1242 code << " CeedScalar* in["<<numinputfields<<"];\n"; 1243 for (CeedInt i = 0; i < numinputfields; i++) { 1244 code << " // ---- Input field "<<i<<" ----\n"; 1245 code << " in["<<i<<"] = r_q"<<i<<";\n"; 1246 } 1247 code << " CeedScalar* out["<<numoutputfields<<"];\n"; 1248 for (CeedInt i = 0; i < numoutputfields; i++) { 1249 code << " // ---- Output field "<<i<<" ----\n"; 1250 code << " out["<<i<<"] = r_qq"<<i<<";\n"; 1251 } 1252 code << "\n // -- Apply QFunction --\n"; 1253 code << " "<<qFunctionName<<"(ctx, "; 1254 if (dim != 3 || useCollograd) { 1255 code << "1"; 1256 } else { 1257 code << "Q1d"; 1258 } 1259 code << ", in, out);\n"; 1260 if (useCollograd) { 1261 code << "\n // Note: Collocated Gradient\n"; 1262 code << " // -- Output fields --\n"; 1263 for (CeedInt i = 0; i < numoutputfields; i++) { 1264 code << " // ---- Output field "<<i<<" ----\n"; 1265 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 1266 CeedChkBackend(ierr); 1267 // Basis action 1268 code << " // EvalMode: "<<CeedEvalModes[emode]<<"\n"; 1269 switch (emode) { 1270 case CEED_EVAL_NONE: 1271 code << " for (CeedInt j = 0; j < ncomp_out_"<<i<<" ; ++j) {\n"; 1272 code << " r_tt"<<i<<"[q + j*Q1d] = r_qq"<<i<<"[j];\n"; 1273 code << " }\n"; 1274 break; // No action 1275 case CEED_EVAL_INTERP: 1276 code << " for (CeedInt j = 0; j < ncomp_out_"<<i<<" ; ++j) {\n"; 1277 code << " r_tt"<<i<<"[q + j*Q1d] = r_qq"<<i<<"[j];\n"; 1278 code << " }\n"; 1279 break; 1280 case CEED_EVAL_GRAD: 1281 code << " gradColloTranspose3d<ncomp_out_"<<i<<",Q1d>(data, q, r_qq"<<i<<", s_G_out_"<<i<<", r_tt"<<i<<");\n"; 1282 break; 1283 case CEED_EVAL_WEIGHT: 1284 break; // Should not occur 1285 case CEED_EVAL_DIV: 1286 break; // TODO: Not implemented 1287 case CEED_EVAL_CURL: 1288 break; // TODO: Not implemented 1289 } 1290 } 1291 code << " }\n"; 1292 } 1293 1294 // Output basis apply if needed 1295 // Generate the correct eval mode code for each output 1296 code << "\n // -- Output field basis action and restrictions --\n"; 1297 for (CeedInt i = 0; i < numoutputfields; i++) { 1298 code << " // ---- Output field "<<i<<" ----\n"; 1299 // Get elemsize, emode, ncomp 1300 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 1301 CeedChkBackend(ierr); 1302 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 1303 CeedChkBackend(ierr); 1304 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 1305 CeedChkBackend(ierr); 1306 ierr = CeedElemRestrictionGetNumComponents(Erestrict, &ncomp); 1307 CeedChkBackend(ierr); 1308 // Basis action 1309 code << " // EvalMode: "<<CeedEvalModes[emode]<<"\n"; 1310 switch (emode) { 1311 case CEED_EVAL_NONE: 1312 code << " CeedScalar* r_v"<<i<<" = r_tt"<<i<<";\n"; 1313 break; // No action 1314 case CEED_EVAL_INTERP: 1315 code << " CeedScalar r_v"<<i<<"[ncomp_out_"<<i<<"*P_out_"<<i<<"];\n"; 1316 code << " interpTranspose"<<dim<<"d<ncomp_out_"<<i<<",P_out_"<<i<<",Q1d>(data, r_tt"<<i<<", s_B_out_"<<i<<", r_v"<<i<<");\n"; 1317 break; 1318 case CEED_EVAL_GRAD: 1319 code << " CeedScalar r_v"<<i<<"[ncomp_out_"<<i<<"*P_out_"<<i<<"];\n"; 1320 if (useCollograd) { 1321 code << " interpTranspose"<<dim<<"d<ncomp_out_"<<i<<",P_out_"<<i<<",Q1d>(data, r_tt"<<i<<", s_B_out_"<<i<<", r_v"<<i<<");\n"; 1322 } else { 1323 code << " gradTranspose"<<dim<<"d<ncomp_out_"<<i<<",P_out_"<<i<<",Q1d>(data, r_tt"<<i<<", s_B_out_"<<i<<", s_G_out_"<<i<<", r_v"<<i<<");\n"; 1324 } 1325 break; 1326 // LCOV_EXCL_START 1327 case CEED_EVAL_WEIGHT: { 1328 Ceed ceed; 1329 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1330 return CeedError(ceed, CEED_ERROR_BACKEND, 1331 "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1332 break; // Should not occur 1333 } 1334 case CEED_EVAL_DIV: 1335 break; // TODO: Not implemented 1336 case CEED_EVAL_CURL: 1337 break; // TODO: Not implemented 1338 // LCOV_EXCL_STOP 1339 } 1340 // Restriction 1341 bool isStrided; 1342 ierr = CeedElemRestrictionIsStrided(Erestrict, &isStrided); CeedChkBackend(ierr); 1343 if (!isStrided) { 1344 ierr = CeedElemRestrictionGetLVectorSize(Erestrict, &lsize); 1345 CeedChkBackend(ierr); 1346 code << " const CeedInt lsize_out_"<<i<<" = "<<lsize<<";\n"; 1347 CeedInt compstride; 1348 ierr = CeedElemRestrictionGetCompStride(Erestrict, &compstride); CeedChkBackend(ierr); 1349 code << " // CompStride: "<<compstride<<"\n"; 1350 ierr = CeedElemRestrictionGetData(Erestrict, &restr_data); CeedChkBackend(ierr); 1351 data->indices.out[i] = restr_data->d_ind; 1352 code << " writeDofsOffset"<<dim<<"d<ncomp_out_"<<i<<", "<<compstride<<", P_out_"<<i<<">(data, lsize_out_"<<i<<", elem, indices.out["<<i<<"], r_v"<<i<<", d_v"<<i<<");\n"; 1353 } else { 1354 bool backendstrides; 1355 ierr = CeedElemRestrictionHasBackendStrides(Erestrict, &backendstrides); 1356 CeedChkBackend(ierr); 1357 CeedInt nelem; 1358 ierr = CeedElemRestrictionGetNumElements(Erestrict, &nelem); 1359 CeedChkBackend(ierr); 1360 CeedInt strides[3] = {1, elemsize*nelem, elemsize}; 1361 if (!backendstrides) { 1362 ierr = CeedElemRestrictionGetStrides(Erestrict, &strides); 1363 CeedChkBackend(ierr); 1364 } 1365 code << " // Strides: {"<<strides[0]<<", "<<strides[1]<<", "<<strides[2]<<"}\n"; 1366 code << " writeDofsStrided"<<dim<<"d<ncomp_out_"<<i<<",P_out_"<<i<<","<<strides[0]<<","<<strides[1]<<","<<strides[2]<<">(data, elem, r_v"<<i<<", d_v"<<i<<");\n"; 1367 } 1368 } 1369 1370 code << " }\n"; 1371 code << "}\n"; 1372 code << "// -----------------------------------------------------------------------------\n\n"; 1373 1374 // View kernel for debugging 1375 CeedDebug(code.str().c_str()); 1376 1377 ierr = CeedCompileCuda(ceed, code.str().c_str(), &data->module, 1, 1378 "T1d", CeedIntMax(Q1d, data->maxP1d)); 1379 CeedChkBackend(ierr); 1380 ierr = CeedGetKernelCuda(ceed, data->module, oper.c_str(), &data->op); 1381 CeedChkBackend(ierr); 1382 1383 ierr = CeedOperatorSetSetupDone(op); CeedChkBackend(ierr); 1384 return CEED_ERROR_SUCCESS; 1385 } 1386 //------------------------------------------------------------------------------ 1387