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