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