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