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