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 #include "ceed-ref.h" 18 19 //------------------------------------------------------------------------------ 20 // Basis Apply 21 //------------------------------------------------------------------------------ 22 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt nelem, 23 CeedTransposeMode tmode, CeedEvalMode emode, 24 CeedVector U, CeedVector V) { 25 int ierr; 26 Ceed ceed; 27 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 28 CeedInt dim, ncomp, nnodes, nqpt; 29 ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr); 30 ierr = CeedBasisGetNumComponents(basis, &ncomp); CeedChk(ierr); 31 ierr = CeedBasisGetNumNodes(basis, &nnodes); CeedChk(ierr); 32 ierr = CeedBasisGetNumQuadraturePoints(basis, &nqpt); CeedChk(ierr); 33 CeedTensorContract contract; 34 ierr = CeedBasisGetTensorContract(basis, &contract); CeedChk(ierr); 35 const CeedInt add = (tmode == CEED_TRANSPOSE); 36 const CeedScalar *u; 37 CeedScalar *v; 38 if (U != CEED_VECTOR_NONE) { 39 ierr = CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); CeedChk(ierr); 40 } else if (emode != CEED_EVAL_WEIGHT) { 41 // LCOV_EXCL_START 42 return CeedError(ceed, 1, 43 "An input vector is required for this CeedEvalMode"); 44 // LCOV_EXCL_STOP 45 } 46 ierr = CeedVectorGetArray(V, CEED_MEM_HOST, &v); CeedChk(ierr); 47 48 // Clear v if operating in transpose 49 if (tmode == CEED_TRANSPOSE) { 50 const CeedInt vsize = nelem*ncomp*nnodes; 51 for (CeedInt i = 0; i < vsize; i++) 52 v[i] = (CeedScalar) 0.0; 53 } 54 bool tensorbasis; 55 ierr = CeedBasisGetTensorStatus(basis, &tensorbasis); CeedChk(ierr); 56 // Tensor basis 57 if (tensorbasis) { 58 CeedInt P1d, Q1d; 59 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr); 60 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChk(ierr); 61 switch (emode) { 62 // Interpolate to/from quadrature points 63 case CEED_EVAL_INTERP: { 64 CeedBasis_Ref *impl; 65 ierr = CeedBasisGetData(basis, (void *)&impl); CeedChk(ierr); 66 if (impl->collointerp) { 67 memcpy(v, u, nelem*ncomp*nnodes*sizeof(u[0])); 68 } else { 69 CeedInt P = P1d, Q = Q1d; 70 if (tmode == CEED_TRANSPOSE) { 71 P = Q1d; Q = P1d; 72 } 73 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem; 74 CeedScalar tmp[2][nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 75 const CeedScalar *interp1d; 76 ierr = CeedBasisGetInterp1D(basis, &interp1d); CeedChk(ierr); 77 for (CeedInt d=0; d<dim; d++) { 78 ierr = CeedTensorContractApply(contract, pre, P, post, Q, 79 interp1d, tmode, add&&(d==dim-1), 80 d==0?u:tmp[d%2], 81 d==dim-1?v:tmp[(d+1)%2]); 82 CeedChk(ierr); 83 pre /= P; 84 post *= Q; 85 } 86 } 87 } break; 88 // Evaluate the gradient to/from quadrature points 89 case CEED_EVAL_GRAD: { 90 // In CEED_NOTRANSPOSE mode: 91 // u has shape [dim, ncomp, P^dim, nelem], row-major layout 92 // v has shape [dim, ncomp, Q^dim, nelem], row-major layout 93 // In CEED_TRANSPOSE mode, the sizes of u and v are switched. 94 CeedInt P = P1d, Q = Q1d; 95 if (tmode == CEED_TRANSPOSE) { 96 P = Q1d, Q = Q1d; 97 } 98 CeedBasis_Ref *impl; 99 ierr = CeedBasisGetData(basis, (void *)&impl); CeedChk(ierr); 100 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem; 101 const CeedScalar *interp1d; 102 ierr = CeedBasisGetInterp1D(basis, &interp1d); CeedChk(ierr); 103 if (impl->collograd1d) { 104 CeedScalar tmp[2][nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 105 CeedScalar interp[nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 106 // Interpolate to quadrature points (NoTranspose) 107 // or Grad to quadrature points (Transpose) 108 for (CeedInt d=0; d<dim; d++) { 109 ierr = CeedTensorContractApply(contract, pre, P, post, Q, 110 (tmode == CEED_NOTRANSPOSE 111 ? interp1d 112 : impl->collograd1d), 113 tmode, add&&(d>0), 114 (tmode == CEED_NOTRANSPOSE 115 ? (d==0?u:tmp[d%2]) 116 : u + d*nqpt*ncomp*nelem), 117 (tmode == CEED_NOTRANSPOSE 118 ? (d==dim-1?interp:tmp[(d+1)%2]) 119 : interp)); 120 CeedChk(ierr); 121 pre /= P; 122 post *= Q; 123 } 124 // Grad to quadrature points (NoTranspose) 125 // or Interpolate to nodes (Transpose) 126 P = Q1d, Q = Q1d; 127 if (tmode == CEED_TRANSPOSE) { 128 P = Q1d, Q = P1d; 129 } 130 pre = ncomp*CeedIntPow(P, dim-1), post = nelem; 131 for (CeedInt d=0; d<dim; d++) { 132 ierr = CeedTensorContractApply(contract, pre, P, post, Q, 133 (tmode == CEED_NOTRANSPOSE 134 ? impl->collograd1d 135 : interp1d), 136 tmode, add&&(d==dim-1), 137 (tmode == CEED_NOTRANSPOSE 138 ? interp 139 : (d==0?interp:tmp[d%2])), 140 (tmode == CEED_NOTRANSPOSE 141 ? v + d*nqpt*ncomp*nelem 142 : (d==dim-1?v:tmp[(d+1)%2]))); 143 CeedChk(ierr); 144 pre /= P; 145 post *= Q; 146 } 147 } else if (impl->collointerp) { // Qpts collocated with nodes 148 const CeedScalar *grad1d; 149 ierr = CeedBasisGetGrad1D(basis, &grad1d); CeedChk(ierr); 150 151 // Dim contractions, identity in other directions 152 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem; 153 for (CeedInt d=0; d<dim; d++) { 154 ierr = CeedTensorContractApply(contract, pre, P, post, Q, 155 grad1d, tmode, add&&(d>0), 156 tmode == CEED_NOTRANSPOSE 157 ? u : u+d*ncomp*nqpt*nelem, 158 tmode == CEED_TRANSPOSE 159 ? v : v+d*ncomp*nqpt*nelem); 160 CeedChk(ierr); 161 pre /= P; 162 post *= Q; 163 } 164 } else { // Underintegration, P > Q 165 const CeedScalar *grad1d; 166 ierr = CeedBasisGetGrad1D(basis, &grad1d); CeedChk(ierr); 167 168 if (tmode == CEED_TRANSPOSE) { 169 P = Q1d, Q = P1d; 170 } 171 CeedScalar tmp[2][nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 172 173 // Dim**2 contractions, apply grad when pass == dim 174 for (CeedInt p=0; p<dim; p++) { 175 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem; 176 for (CeedInt d=0; d<dim; d++) { 177 ierr = CeedTensorContractApply(contract, pre, P, post, Q, 178 (p==d)? grad1d : interp1d, 179 tmode, add&&(d==dim-1), 180 (d == 0 181 ? (tmode == CEED_NOTRANSPOSE 182 ? u : u+p*ncomp*nqpt*nelem) 183 : tmp[d%2]), 184 (d == dim-1 185 ? (tmode == CEED_TRANSPOSE 186 ? v : v+p*ncomp*nqpt*nelem) 187 : tmp[(d+1)%2])); 188 CeedChk(ierr); 189 pre /= P; 190 post *= Q; 191 } 192 } 193 } 194 } break; 195 // Retrieve interpolation weights 196 case CEED_EVAL_WEIGHT: { 197 if (tmode == CEED_TRANSPOSE) 198 // LCOV_EXCL_START 199 return CeedError(ceed, 1, 200 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 201 // LCOV_EXCL_STOP 202 CeedInt Q = Q1d; 203 const CeedScalar *qweight1d; 204 ierr = CeedBasisGetQWeights(basis, &qweight1d); CeedChk(ierr); 205 for (CeedInt d=0; d<dim; d++) { 206 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d); 207 for (CeedInt i=0; i<pre; i++) 208 for (CeedInt j=0; j<Q; j++) 209 for (CeedInt k=0; k<post; k++) { 210 CeedScalar w = qweight1d[j] 211 * (d == 0 ? 1 : v[((i*Q + j)*post + k)*nelem]); 212 for (CeedInt e=0; e<nelem; e++) 213 v[((i*Q + j)*post + k)*nelem + e] = w; 214 } 215 } 216 } break; 217 // LCOV_EXCL_START 218 // Evaluate the divergence to/from the quadrature points 219 case CEED_EVAL_DIV: 220 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported"); 221 // Evaluate the curl to/from the quadrature points 222 case CEED_EVAL_CURL: 223 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported"); 224 // Take no action, BasisApply should not have been called 225 case CEED_EVAL_NONE: 226 return CeedError(ceed, 1, 227 "CEED_EVAL_NONE does not make sense in this context"); 228 // LCOV_EXCL_STOP 229 } 230 } else { 231 // Non-tensor basis 232 switch (emode) { 233 // Interpolate to/from quadrature points 234 case CEED_EVAL_INTERP: { 235 CeedInt P = nnodes, Q = nqpt; 236 const CeedScalar *interp; 237 ierr = CeedBasisGetInterp(basis, &interp); CeedChk(ierr); 238 if (tmode == CEED_TRANSPOSE) { 239 P = nqpt; Q = nnodes; 240 } 241 ierr = CeedTensorContractApply(contract, ncomp, P, nelem, Q, 242 interp, tmode, add, u, v); 243 CeedChk(ierr); 244 } 245 break; 246 // Evaluate the gradient to/from quadrature points 247 case CEED_EVAL_GRAD: { 248 CeedInt P = nnodes, Q = nqpt; 249 CeedInt dimstride = nqpt * ncomp * nelem; 250 CeedInt gradstride = nqpt * nnodes; 251 const CeedScalar *grad; 252 ierr = CeedBasisGetGrad(basis, &grad); CeedChk(ierr); 253 if (tmode == CEED_TRANSPOSE) { 254 P = nqpt; Q = nnodes; 255 for (CeedInt d = 0; d < dim; d++) { 256 ierr = CeedTensorContractApply(contract, ncomp, P, nelem, Q, 257 grad + d * gradstride, tmode, add, 258 u + d * dimstride, v); CeedChk(ierr); 259 } 260 } else { 261 for (CeedInt d = 0; d < dim; d++) { 262 ierr = CeedTensorContractApply(contract, ncomp, P, nelem, Q, 263 grad + d * gradstride, tmode, add, 264 u, v + d * dimstride); CeedChk(ierr); 265 } 266 } 267 } 268 break; 269 // Retrieve interpolation weights 270 case CEED_EVAL_WEIGHT: { 271 if (tmode == CEED_TRANSPOSE) 272 // LCOV_EXCL_START 273 return CeedError(ceed, 1, 274 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 275 // LCOV_EXCL_STOP 276 const CeedScalar *qweight; 277 ierr = CeedBasisGetQWeights(basis, &qweight); CeedChk(ierr); 278 for (CeedInt i=0; i<nqpt; i++) 279 for (CeedInt e=0; e<nelem; e++) 280 v[i*nelem + e] = qweight[i]; 281 } break; 282 // LCOV_EXCL_START 283 // Evaluate the divergence to/from the quadrature points 284 case CEED_EVAL_DIV: 285 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported"); 286 // Evaluate the curl to/from the quadrature points 287 case CEED_EVAL_CURL: 288 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported"); 289 // Take no action, BasisApply should not have been called 290 case CEED_EVAL_NONE: 291 return CeedError(ceed, 1, 292 "CEED_EVAL_NONE does not make sense in this context"); 293 // LCOV_EXCL_STOP 294 } 295 } 296 if (U != CEED_VECTOR_NONE) { 297 ierr = CeedVectorRestoreArrayRead(U, &u); CeedChk(ierr); 298 } 299 ierr = CeedVectorRestoreArray(V, &v); CeedChk(ierr); 300 return 0; 301 } 302 303 //------------------------------------------------------------------------------ 304 // Basis Destroy Non-Tensor 305 //------------------------------------------------------------------------------ 306 static int CeedBasisDestroyNonTensor_Ref(CeedBasis basis) { 307 int ierr; 308 CeedTensorContract contract; 309 ierr = CeedBasisGetTensorContract(basis, &contract); CeedChk(ierr); 310 ierr = CeedTensorContractDestroy(&contract); CeedChk(ierr); 311 return 0; 312 } 313 314 //------------------------------------------------------------------------------ 315 // Basis Create Non-Tensor 316 //------------------------------------------------------------------------------ 317 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, 318 CeedInt nnodes, CeedInt nqpts, 319 const CeedScalar *interp, 320 const CeedScalar *grad, 321 const CeedScalar *qref, 322 const CeedScalar *qweight, 323 CeedBasis basis) { 324 int ierr; 325 Ceed ceed; 326 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 327 328 Ceed parent; 329 ierr = CeedGetParent(ceed, &parent); CeedChk(ierr); 330 CeedTensorContract contract; 331 ierr = CeedTensorContractCreate(parent, basis, &contract); CeedChk(ierr); 332 ierr = CeedBasisSetTensorContract(basis, &contract); CeedChk(ierr); 333 334 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply", 335 CeedBasisApply_Ref); CeedChk(ierr); 336 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", 337 CeedBasisDestroyNonTensor_Ref); CeedChk(ierr); 338 339 return 0; 340 } 341 342 //------------------------------------------------------------------------------ 343 // Basis Destroy Tensor 344 //------------------------------------------------------------------------------ 345 static int CeedBasisDestroyTensor_Ref(CeedBasis basis) { 346 int ierr; 347 CeedTensorContract contract; 348 ierr = CeedBasisGetTensorContract(basis, &contract); CeedChk(ierr); 349 ierr = CeedTensorContractDestroy(&contract); CeedChk(ierr); 350 351 CeedBasis_Ref *impl; 352 ierr = CeedBasisGetData(basis, (void *)&impl); CeedChk(ierr); 353 ierr = CeedFree(&impl->collograd1d); CeedChk(ierr); 354 ierr = CeedFree(&impl); CeedChk(ierr); 355 356 return 0; 357 } 358 359 //------------------------------------------------------------------------------ 360 // Basis Create Tensor 361 //------------------------------------------------------------------------------ 362 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P1d, 363 CeedInt Q1d, const CeedScalar *interp1d, 364 const CeedScalar *grad1d, 365 const CeedScalar *qref1d, 366 const CeedScalar *qweight1d, 367 CeedBasis basis) { 368 int ierr; 369 Ceed ceed; 370 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 371 CeedBasis_Ref *impl; 372 ierr = CeedCalloc(1, &impl); CeedChk(ierr); 373 // Check for collocated interp 374 if (Q1d == P1d) { 375 bool collocated = 1; 376 for (CeedInt i=0; i<P1d; i++) { 377 collocated = collocated && (fabs(interp1d[i+P1d*i] - 1.0) < 1e-14); 378 for (CeedInt j=0; j<P1d; j++) 379 if (j != i) 380 collocated = collocated && (fabs(interp1d[j+P1d*i]) < 1e-14); 381 } 382 impl->collointerp = collocated; 383 } 384 // Calculate collocated grad 385 if (Q1d >= P1d && !impl->collointerp) { 386 ierr = CeedMalloc(Q1d*Q1d, &impl->collograd1d); CeedChk(ierr); 387 ierr = CeedBasisGetCollocatedGrad(basis, impl->collograd1d); CeedChk(ierr); 388 } 389 ierr = CeedBasisSetData(basis, (void *)&impl); CeedChk(ierr); 390 391 Ceed parent; 392 ierr = CeedGetParent(ceed, &parent); CeedChk(ierr); 393 CeedTensorContract contract; 394 ierr = CeedTensorContractCreate(parent, basis, &contract); CeedChk(ierr); 395 ierr = CeedBasisSetTensorContract(basis, &contract); CeedChk(ierr); 396 397 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply", 398 CeedBasisApply_Ref); CeedChk(ierr); 399 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", 400 CeedBasisDestroyTensor_Ref); CeedChk(ierr); 401 return 0; 402 } 403 404 //------------------------------------------------------------------------------ 405