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