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 <string.h> 18 #include "ceed-ref.h" 19 20 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt nelem, 21 CeedTransposeMode tmode, CeedEvalMode emode, 22 CeedVector U, CeedVector V) { 23 int ierr; 24 Ceed ceed; 25 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 26 CeedInt dim, ncomp, ndof, nqpt; 27 ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr); 28 ierr = CeedBasisGetNumComponents(basis, &ncomp); CeedChk(ierr); 29 ierr = CeedBasisGetNumNodes(basis, &ndof); CeedChk(ierr); 30 ierr = CeedBasisGetNumQuadraturePoints(basis, &nqpt); CeedChk(ierr); 31 CeedTensorContract contract; 32 ierr = CeedBasisGetTensorContract(basis, &contract); CeedChk(ierr); 33 const CeedInt add = (tmode == CEED_TRANSPOSE); 34 const CeedScalar *u; 35 CeedScalar *v; 36 if (U) { 37 ierr = CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); CeedChk(ierr); 38 } else if (emode != CEED_EVAL_WEIGHT) { 39 return CeedError(ceed, 1, 40 "An input vector is required for this CeedEvalMode"); 41 } 42 ierr = CeedVectorGetArray(V, CEED_MEM_HOST, &v); CeedChk(ierr); 43 44 // Clear v if operating in transpose 45 if (tmode == CEED_TRANSPOSE) { 46 const CeedInt vsize = nelem*ncomp*ndof; 47 for (CeedInt i = 0; i < vsize; i++) 48 v[i] = (CeedScalar) 0.0; 49 } 50 bool tensorbasis; 51 ierr = CeedBasisGetTensorStatus(basis, &tensorbasis); CeedChk(ierr); 52 // Tensor basis 53 if (tensorbasis) { 54 CeedInt P1d, Q1d; 55 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr); 56 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChk(ierr); 57 switch (emode) { 58 // Interpolate to/from quadrature points 59 case CEED_EVAL_INTERP: { 60 CeedInt P = P1d, Q = Q1d; 61 if (tmode == CEED_TRANSPOSE) { 62 P = Q1d; Q = P1d; 63 } 64 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem; 65 CeedScalar tmp[2][nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 66 CeedScalar *interp1d; 67 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr); 68 for (CeedInt d=0; d<dim; d++) { 69 ierr = CeedTensorContractApply(contract, pre, P, post, Q, 70 interp1d, tmode, add&&(d==dim-1), 71 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]); 72 CeedChk(ierr); 73 pre /= P; 74 post *= Q; 75 } 76 } break; 77 // Evaluate the gradient to/from quadrature points 78 case CEED_EVAL_GRAD: { 79 // In CEED_NOTRANSPOSE mode: 80 // u has shape [dim, ncomp, P^dim, nelem], row-major layout 81 // v has shape [dim, ncomp, Q^dim, nelem], row-major layout 82 // In CEED_TRANSPOSE mode, the sizes of u and v are switched. 83 CeedInt P = P1d, Q = Q1d; 84 if (tmode == CEED_TRANSPOSE) { 85 P = Q1d, Q = Q1d; 86 } 87 CeedBasis_Ref *impl; 88 ierr = CeedBasisGetData(basis, (void *)&impl); CeedChk(ierr); 89 CeedScalar interp[nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 90 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = nelem; 91 CeedScalar tmp[2][nelem*ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 92 CeedScalar *interp1d; 93 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr); 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 } break; 136 // Retrieve interpolation weights 137 case CEED_EVAL_WEIGHT: { 138 if (tmode == CEED_TRANSPOSE) 139 return CeedError(ceed, 1, 140 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 141 CeedInt Q = Q1d; 142 CeedScalar *qweight1d; 143 ierr = CeedBasisGetQWeights(basis, &qweight1d); CeedChk(ierr); 144 for (CeedInt d=0; d<dim; d++) { 145 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d); 146 for (CeedInt i=0; i<pre; i++) 147 for (CeedInt j=0; j<Q; j++) 148 for (CeedInt k=0; k<post; k++) { 149 CeedScalar w = qweight1d[j] 150 * (d == 0 ? 1 : v[((i*Q + j)*post + k)*nelem]); 151 for (CeedInt e=0; e<nelem; e++) 152 v[((i*Q + j)*post + k)*nelem + e] = w; 153 } 154 } 155 } break; 156 // Evaluate the divergence to/from the quadrature points 157 case CEED_EVAL_DIV: 158 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported"); 159 // Evaluate the curl to/from the quadrature points 160 case CEED_EVAL_CURL: 161 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported"); 162 // Take no action, BasisApply should not have been called 163 case CEED_EVAL_NONE: 164 return CeedError(ceed, 1, 165 "CEED_EVAL_NONE does not make sense in this context"); 166 } 167 } else { 168 // Non-tensor basis 169 switch (emode) { 170 // Interpolate to/from quadrature points 171 case CEED_EVAL_INTERP: { 172 CeedInt P = ndof, Q = nqpt; 173 CeedScalar *interp; 174 ierr = CeedBasisGetInterp(basis, &interp); CeedChk(ierr); 175 if (tmode == CEED_TRANSPOSE) { 176 P = nqpt; Q = ndof; 177 } 178 ierr = CeedTensorContractApply(contract, ncomp, P, nelem, Q, 179 interp, tmode, add, u, v); 180 CeedChk(ierr); 181 } 182 break; 183 // Evaluate the gradient to/from quadrature points 184 case CEED_EVAL_GRAD: { 185 CeedInt P = ndof, Q = dim*nqpt; 186 CeedScalar *grad; 187 ierr = CeedBasisGetGrad(basis, &grad); CeedChk(ierr); 188 if (tmode == CEED_TRANSPOSE) { 189 P = dim*nqpt; Q = ndof; 190 } 191 ierr = CeedTensorContractApply(contract, ncomp, P, nelem, Q, 192 grad, tmode, add, u, v); 193 CeedChk(ierr); 194 } 195 break; 196 // Retrieve interpolation weights 197 case CEED_EVAL_WEIGHT: { 198 if (tmode == CEED_TRANSPOSE) 199 return CeedError(ceed, 1, 200 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 201 CeedScalar *qweight; 202 ierr = CeedBasisGetQWeights(basis, &qweight); CeedChk(ierr); 203 for (CeedInt i=0; i<nqpt; i++) 204 for (CeedInt e=0; e<nelem; e++) 205 v[i*nelem + e] = qweight[i]; 206 } break; 207 // Evaluate the divergence to/from the quadrature points 208 case CEED_EVAL_DIV: 209 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported"); 210 // Evaluate the curl to/from the quadrature points 211 case CEED_EVAL_CURL: 212 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported"); 213 // Take no action, BasisApply should not have been called 214 case CEED_EVAL_NONE: 215 return CeedError(ceed, 1, 216 "CEED_EVAL_NONE does not make sense in this context"); 217 } 218 } 219 if (U) { 220 ierr = CeedVectorRestoreArrayRead(U, &u); CeedChk(ierr); 221 } 222 ierr = CeedVectorRestoreArray(V, &v); CeedChk(ierr); 223 return 0; 224 } 225 226 static int CeedBasisDestroyNonTensor_Ref(CeedBasis basis) { 227 int ierr; 228 CeedTensorContract contract; 229 ierr = CeedBasisGetTensorContract(basis, &contract); CeedChk(ierr); 230 ierr = CeedTensorContractDestroy(&contract); CeedChk(ierr); 231 return 0; 232 } 233 234 static int CeedBasisDestroyTensor_Ref(CeedBasis basis) { 235 int ierr; 236 CeedTensorContract contract; 237 ierr = CeedBasisGetTensorContract(basis, &contract); CeedChk(ierr); 238 ierr = CeedTensorContractDestroy(&contract); CeedChk(ierr); 239 240 CeedBasis_Ref *impl; 241 ierr = CeedBasisGetData(basis, (void *)&impl); CeedChk(ierr); 242 ierr = CeedFree(&impl->colograd1d); CeedChk(ierr); 243 ierr = CeedFree(&impl); CeedChk(ierr); 244 245 return 0; 246 } 247 248 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P1d, 249 CeedInt Q1d, const CeedScalar *interp1d, 250 const CeedScalar *grad1d, 251 const CeedScalar *qref1d, 252 const CeedScalar *qweight1d, 253 CeedBasis basis) { 254 int ierr; 255 Ceed ceed; 256 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 257 CeedBasis_Ref *impl; 258 ierr = CeedCalloc(1, &impl); CeedChk(ierr); 259 ierr = CeedMalloc(Q1d*Q1d, &impl->colograd1d); CeedChk(ierr); 260 ierr = CeedBasisGetCollocatedGrad(basis, impl->colograd1d); CeedChk(ierr); 261 ierr = CeedBasisSetData(basis, (void *)&impl); CeedChk(ierr); 262 263 Ceed parent; 264 ierr = CeedGetParent(ceed, &parent); CeedChk(ierr); 265 CeedTensorContract contract; 266 ierr = CeedTensorContractCreate(parent, &contract); CeedChk(ierr); 267 ierr = CeedBasisSetTensorContract(basis, &contract); CeedChk(ierr); 268 269 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply", 270 CeedBasisApply_Ref); CeedChk(ierr); 271 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", 272 CeedBasisDestroyTensor_Ref); CeedChk(ierr); 273 return 0; 274 } 275 276 277 278 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, 279 CeedInt ndof, CeedInt nqpts, 280 const CeedScalar *interp, 281 const CeedScalar *grad, 282 const CeedScalar *qref, 283 const CeedScalar *qweight, 284 CeedBasis basis) { 285 int ierr; 286 Ceed ceed; 287 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 288 289 CeedTensorContract contract; 290 ierr = CeedTensorContractCreate(ceed, &contract); CeedChk(ierr); 291 ierr = CeedBasisSetTensorContract(basis, &contract); CeedChk(ierr); 292 293 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply", 294 CeedBasisApply_Ref); CeedChk(ierr); 295 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", 296 CeedBasisDestroyNonTensor_Ref); CeedChk(ierr); 297 298 return 0; 299 } 300