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-impl.h> 18 #include <string.h> 19 #include "ceed-ref.h" 20 21 // Contracts on the middle index 22 // NOTRANSPOSE: V_ajc = T_jb U_abc 23 // TRANSPOSE: V_ajc = T_bj U_abc 24 // If Add != 0, "=" is replaced by "+=" 25 static int CeedTensorContract_Ref(Ceed ceed, 26 CeedInt A, CeedInt B, CeedInt C, CeedInt J, 27 const CeedScalar *restrict t, CeedTransposeMode tmode, 28 const CeedInt Add, 29 const CeedScalar *restrict u, CeedScalar *restrict v) { 30 CeedInt tstride0 = B, tstride1 = 1; 31 if (tmode == CEED_TRANSPOSE) { 32 tstride0 = 1; tstride1 = J; 33 } 34 35 if (!Add) 36 for (CeedInt q=0; q<A*J*C; q++) 37 v[q] = (CeedScalar) 0.0; 38 39 for (CeedInt a=0; a<A; a++) 40 for (CeedInt b=0; b<B; b++) 41 for (CeedInt j=0; j<J; j++) { 42 CeedScalar tq = t[j*tstride0 + b*tstride1]; 43 for (CeedInt c=0; c<C; c++) 44 v[(a*J+j)*C+c] += tq * u[(a*B+b)*C+c]; 45 } 46 return 0; 47 } 48 49 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt nelem, 50 CeedTransposeMode tmode, CeedEvalMode emode, 51 const CeedScalar *u, CeedScalar *v) { 52 int ierr; 53 const CeedInt dim = basis->dim; 54 const CeedInt ncomp = basis->ncomp; 55 const CeedInt nqpt = CeedPowInt(basis->Q1d, dim); 56 const CeedInt add = (tmode == CEED_TRANSPOSE); 57 58 if (nelem != 1) 59 return CeedError(basis->ceed, 1, 60 "This backend does not support BasisApply for multiple elements"); 61 62 // Clear v if operating in transpose 63 if (tmode == CEED_TRANSPOSE) { 64 const CeedInt vsize = ncomp*CeedPowInt(basis->P1d, dim); 65 for (CeedInt i = 0; i < vsize; i++) 66 v[i] = (CeedScalar) 0; 67 } 68 switch (emode) { 69 // Interpolate to/from quadrature points 70 case CEED_EVAL_INTERP: { 71 CeedInt P = basis->P1d, Q = basis->Q1d; 72 if (tmode == CEED_TRANSPOSE) { 73 P = basis->Q1d; Q = basis->P1d; 74 } 75 CeedInt pre = ncomp*CeedPowInt(P, dim-1), post = 1; 76 CeedScalar tmp[2][ncomp*Q*CeedPowInt(P>Q?P:Q, dim-1)]; 77 for (CeedInt d=0; d<dim; d++) { 78 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q, basis->interp1d, 79 tmode, add&&(d==dim-1), 80 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]); 81 CeedChk(ierr); 82 pre /= P; 83 post *= Q; 84 } 85 } break; 86 // Evaluate the gradient to/from quadrature points 87 case CEED_EVAL_GRAD: { 88 CeedInt P = basis->P1d, Q = basis->Q1d; 89 // In CEED_NOTRANSPOSE mode: 90 // u is [dim, ncomp, P^dim, nelem], row-major layout 91 // v is [dim, ncomp, Q^dim, nelem], row-major layout 92 // In CEED_TRANSPOSE mode, the sizes of u and v are switched. 93 if (tmode == CEED_TRANSPOSE) { 94 P = basis->Q1d, Q = basis->P1d; 95 } 96 CeedScalar tmp[2][ncomp*Q*CeedPowInt(P>Q?P:Q, dim-1)]; 97 for (CeedInt p = 0; p < dim; p++) { 98 CeedInt pre = ncomp*CeedPowInt(P, dim-1), post = 1; 99 for (CeedInt d=0; d<dim; d++) { 100 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q, 101 (p==d)?basis->grad1d:basis->interp1d, 102 tmode, add&&(d==dim-1), 103 d==0 104 ? (tmode==CEED_NOTRANSPOSE?u:u+p*ncomp*nqpt) 105 : tmp[d%2], 106 d==dim-1 107 ? (tmode==CEED_TRANSPOSE?v:v+p*ncomp*nqpt) 108 : tmp[(d+1)%2]); 109 CeedChk(ierr); 110 pre /= P; 111 post *= Q; 112 } 113 } 114 } break; 115 // Retrieve interpolation weights 116 case CEED_EVAL_WEIGHT: { 117 if (tmode == CEED_TRANSPOSE) 118 return CeedError(basis->ceed, 1, 119 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 120 CeedInt Q = basis->Q1d; 121 for (CeedInt d=0; d<dim; d++) { 122 CeedInt pre = CeedPowInt(Q, dim-d-1), post = CeedPowInt(Q, d); 123 for (CeedInt i=0; i<pre; i++) 124 for (CeedInt j=0; j<Q; j++) 125 for (CeedInt k=0; k<post; k++) 126 v[(i*Q + j)*post + k] = basis->qweight1d[j] 127 * (d == 0 ? 1 : v[(i*Q + j)*post + k]); 128 } 129 } break; 130 // Evaluate the divergence to/from the quadrature points 131 case CEED_EVAL_DIV: 132 return CeedError(basis->ceed, 1, "CEED_EVAL_DIV not supported"); 133 // Evaluate the curl to/from the quadrature points 134 case CEED_EVAL_CURL: 135 return CeedError(basis->ceed, 1, "CEED_EVAL_CURL not supported"); 136 // Take no action, BasisApply should not have been called 137 case CEED_EVAL_NONE: 138 return CeedError(basis->ceed, 1, "CEED_EVAL_NONE does not make sense in this context"); 139 } 140 return 0; 141 } 142 143 static int CeedBasisDestroy_Ref(CeedBasis basis) { 144 return 0; 145 } 146 147 int CeedBasisCreateTensorH1_Ref(Ceed ceed, CeedInt dim, CeedInt P1d, 148 CeedInt Q1d, const CeedScalar *interp1d, 149 const CeedScalar *grad1d, 150 const CeedScalar *qref1d, 151 const CeedScalar *qweight1d, 152 CeedBasis basis) { 153 basis->Apply = CeedBasisApply_Ref; 154 basis->Destroy = CeedBasisDestroy_Ref; 155 return 0; 156 } 157