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