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 // Contracts on the middle index 21 // NOTRANSPOSE: V_ajc = T_jb U_abc 22 // TRANSPOSE: V_ajc = T_bj U_abc 23 // If Add != 0, "=" is replaced by "+=" 24 static int CeedTensorContract_Ref(Ceed ceed, CeedInt A, CeedInt B, CeedInt C, CeedInt J, 25 const CeedScalar *restrict t, CeedTransposeMode tmode, 26 const CeedInt Add, 27 const CeedScalar *restrict u, CeedScalar *restrict v) { 28 CeedInt tstride0 = B, tstride1 = 1; 29 if (tmode == CEED_TRANSPOSE) { 30 tstride0 = 1; tstride1 = J; 31 } 32 33 if (!Add) 34 for (CeedInt q=0; q<A*J*C; q++) 35 v[q] = (CeedScalar) 0.0; 36 37 for (CeedInt a=0; a<A; a++) 38 for (CeedInt b=0; b<B; b++) 39 for (CeedInt j=0; j<J; j++) { 40 CeedScalar tq = t[j*tstride0 + b*tstride1]; 41 for (CeedInt c=0; c<C; c++) 42 v[(a*J+j)*C+c] += tq * u[(a*B+b)*C+c]; 43 } 44 return 0; 45 } 46 47 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt nelem, 48 CeedTransposeMode tmode, CeedEvalMode emode, 49 const CeedScalar *u, CeedScalar *v) { 50 int ierr; 51 Ceed ceed; 52 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 53 CeedInt dim, ncomp, ndof, nqpt; 54 ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr); 55 ierr = CeedBasisGetNumComponents(basis, &ncomp); CeedChk(ierr); 56 ierr = CeedBasisGetNumNodes(basis, &ndof); CeedChk(ierr); 57 ierr = CeedBasisGetNumQuadraturePoints(basis, &nqpt); CeedChk(ierr); 58 const CeedInt add = (tmode == CEED_TRANSPOSE); 59 60 if (nelem != 1) 61 return CeedError(ceed, 1, 62 "This backend does not support BasisApply for multiple elements"); 63 64 // Clear v if operating in transpose 65 if (tmode == CEED_TRANSPOSE) { 66 const CeedInt vsize = ncomp*ndof; 67 for (CeedInt i = 0; i < vsize; i++) 68 v[i] = (CeedScalar) 0; 69 } 70 // Tensor basis 71 bool tensorbasis; 72 ierr = CeedBasisGetTensorStatus(basis, &tensorbasis); CeedChk(ierr); 73 if (tensorbasis) { 74 CeedInt P1d, Q1d; 75 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr); 76 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChk(ierr); 77 switch (emode) { 78 // Interpolate to/from quadrature points 79 case CEED_EVAL_INTERP: { 80 CeedInt P = P1d, Q = Q1d; 81 if (tmode == CEED_TRANSPOSE) { 82 P = Q1d; Q = P1d; 83 } 84 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1; 85 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 86 CeedScalar *interp1d; 87 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr); 88 for (CeedInt d=0; d<dim; d++) { 89 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q, interp1d, 90 tmode, add&&(d==dim-1), 91 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]); 92 CeedChk(ierr); 93 pre /= P; 94 post *= Q; 95 } 96 } break; 97 // Evaluate the gradient to/from quadrature points 98 case CEED_EVAL_GRAD: { 99 CeedInt P = P1d, Q = Q1d; 100 // In CEED_NOTRANSPOSE mode: 101 // u has shape [dim, ncomp, P^dim, nelem], row-major layout 102 // v has shape [dim, ncomp, Q^dim, nelem], row-major layout 103 // In CEED_TRANSPOSE mode, the sizes of u and v are switched. 104 if (tmode == CEED_TRANSPOSE) { 105 P = Q1d, Q = P1d; 106 } 107 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 108 CeedScalar *interp1d, *grad1d; 109 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr); 110 ierr = CeedBasisGetGrad(basis, &grad1d); CeedChk(ierr); 111 for (CeedInt p = 0; p < dim; p++) { 112 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1; 113 for (CeedInt d=0; d<dim; d++) { 114 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q, 115 (p==d)?grad1d:interp1d, 116 tmode, add&&(d==dim-1), 117 (d == 0 118 ? (tmode==CEED_NOTRANSPOSE?u:u+p*ncomp*nqpt) 119 : tmp[d%2]), 120 (d == dim-1 121 ? (tmode==CEED_TRANSPOSE?v:v+p*ncomp*nqpt) 122 : tmp[(d+1)%2])); 123 CeedChk(ierr); 124 pre /= P; 125 post *= Q; 126 } 127 } 128 } break; 129 // Retrieve interpolation weights 130 case CEED_EVAL_WEIGHT: { 131 if (tmode == CEED_TRANSPOSE) 132 return CeedError(ceed, 1, 133 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 134 CeedInt Q = Q1d; 135 CeedScalar *qweight1d; 136 ierr = CeedBasisGetQWeights(basis, &qweight1d); CeedChk(ierr); 137 for (CeedInt d=0; d<dim; d++) { 138 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d); 139 for (CeedInt i=0; i<pre; i++) 140 for (CeedInt j=0; j<Q; j++) 141 for (CeedInt k=0; k<post; k++) 142 v[(i*Q + j)*post + k] = qweight1d[j] 143 * (d == 0 ? 1 : v[(i*Q + j)*post + k]); 144 } 145 } break; 146 // Evaluate the divergence to/from the quadrature points 147 case CEED_EVAL_DIV: 148 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported"); 149 // Evaluate the curl to/from the quadrature points 150 case CEED_EVAL_CURL: 151 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported"); 152 // Take no action, BasisApply should not have been called 153 case CEED_EVAL_NONE: 154 return CeedError(ceed, 1, 155 "CEED_EVAL_NONE does not make sense in this context"); 156 } 157 } else { 158 // Non-tensor basis 159 switch (emode) { 160 case CEED_EVAL_INTERP: { 161 CeedInt P = ndof, Q = nqpt; 162 CeedScalar *interp; 163 ierr = CeedBasisGetInterp(basis, &interp); CeedChk(ierr); 164 if (tmode == CEED_TRANSPOSE) { 165 P = nqpt; Q = ndof; 166 } 167 ierr = CeedTensorContract_Ref(ceed, ncomp, P, 1, Q, interp, 168 tmode, add, u, v); 169 CeedChk(ierr); 170 } 171 break; 172 case CEED_EVAL_GRAD: { 173 CeedInt P = ndof, Q = dim*nqpt; 174 CeedScalar *grad; 175 ierr = CeedBasisGetGrad(basis, &grad); CeedChk(ierr); 176 if (tmode == CEED_TRANSPOSE) { 177 P = dim*nqpt; Q = ndof; 178 } 179 ierr = CeedTensorContract_Ref(ceed, ncomp, P, 1, Q, grad, 180 tmode, add, u, v); 181 CeedChk(ierr); 182 } 183 break; 184 case CEED_EVAL_WEIGHT: { 185 if (tmode == CEED_TRANSPOSE) 186 return CeedError(ceed, 1, 187 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 188 CeedScalar *qweight; 189 ierr = CeedBasisGetQWeights(basis, &qweight); CeedChk(ierr); 190 for (CeedInt i=0; i<nqpt; i++) 191 v[i] = qweight[i]; 192 } break; 193 case CEED_EVAL_DIV: 194 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported"); 195 case CEED_EVAL_CURL: 196 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported"); 197 case CEED_EVAL_NONE: 198 return CeedError(ceed, 1, 199 "CEED_EVAL_NONE does not make sense in this context"); 200 } 201 } 202 return 0; 203 } 204 205 static int CeedBasisDestroy_Ref(CeedBasis basis) { 206 return 0; 207 } 208 209 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P1d, 210 CeedInt Q1d, const CeedScalar *interp1d, 211 const CeedScalar *grad1d, 212 const CeedScalar *qref1d, 213 const CeedScalar *qweight1d, 214 CeedBasis basis) { 215 basis->Apply = CeedBasisApply_Ref; 216 basis->Destroy = CeedBasisDestroy_Ref; 217 return 0; 218 } 219 220 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, 221 CeedInt ndof, CeedInt nqpts, 222 const CeedScalar *interp, 223 const CeedScalar *grad, 224 const CeedScalar *qref, 225 const CeedScalar *qweight, 226 CeedBasis basis) { 227 basis->Apply = CeedBasisApply_Ref; 228 basis->Destroy = CeedBasisDestroy_Ref; 229 return 0; 230 } 231