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 CeedVector U, CeedVector 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 const CeedScalar *u; 61 CeedScalar *v; 62 if (U) { 63 ierr = CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); CeedChk(ierr); 64 } else if (emode != CEED_EVAL_WEIGHT) { 65 return CeedError(ceed, 1, 66 "An input vector is required for this CeedEvalMode"); 67 } 68 ierr = CeedVectorGetArray(V, CEED_MEM_HOST, &v); CeedChk(ierr); 69 70 if (nelem != 1) 71 return CeedError(ceed, 1, 72 "This backend does not support BasisApply for multiple elements"); 73 74 // Clear v if operating in transpose 75 if (tmode == CEED_TRANSPOSE) { 76 const CeedInt vsize = ncomp*ndof; 77 for (CeedInt i = 0; i < vsize; i++) 78 v[i] = (CeedScalar) 0; 79 } 80 // Tensor basis 81 bool tensorbasis; 82 ierr = CeedBasisGetTensorStatus(basis, &tensorbasis); CeedChk(ierr); 83 if (tensorbasis) { 84 CeedInt P1d, Q1d; 85 ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr); 86 ierr = CeedBasisGetNumQuadraturePoints1D(basis, &Q1d); CeedChk(ierr); 87 switch (emode) { 88 // Interpolate to/from quadrature points 89 case CEED_EVAL_INTERP: { 90 CeedInt P = P1d, Q = Q1d; 91 if (tmode == CEED_TRANSPOSE) { 92 P = Q1d; Q = P1d; 93 } 94 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1; 95 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 96 CeedScalar *interp1d; 97 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr); 98 for (CeedInt d=0; d<dim; d++) { 99 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q, interp1d, 100 tmode, add&&(d==dim-1), 101 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]); 102 CeedChk(ierr); 103 pre /= P; 104 post *= Q; 105 } 106 } break; 107 // Evaluate the gradient to/from quadrature points 108 case CEED_EVAL_GRAD: { 109 CeedInt P = P1d, Q = Q1d; 110 // In CEED_NOTRANSPOSE mode: 111 // u has shape [dim, ncomp, P^dim, nelem], row-major layout 112 // v has shape [dim, ncomp, Q^dim, nelem], row-major layout 113 // In CEED_TRANSPOSE mode, the sizes of u and v are switched. 114 if (tmode == CEED_TRANSPOSE) { 115 P = Q1d, Q = P1d; 116 } 117 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 118 CeedScalar *interp1d, *grad1d; 119 ierr = CeedBasisGetInterp(basis, &interp1d); CeedChk(ierr); 120 ierr = CeedBasisGetGrad(basis, &grad1d); CeedChk(ierr); 121 for (CeedInt p = 0; p < dim; p++) { 122 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1; 123 for (CeedInt d=0; d<dim; d++) { 124 ierr = CeedTensorContract_Ref(ceed, pre, P, post, Q, 125 (p==d)?grad1d:interp1d, 126 tmode, add&&(d==dim-1), 127 (d == 0 128 ? (tmode==CEED_NOTRANSPOSE?u:u+p*ncomp*nqpt) 129 : tmp[d%2]), 130 (d == dim-1 131 ? (tmode==CEED_TRANSPOSE?v:v+p*ncomp*nqpt) 132 : tmp[(d+1)%2])); 133 CeedChk(ierr); 134 pre /= P; 135 post *= Q; 136 } 137 } 138 } break; 139 // Retrieve interpolation weights 140 case CEED_EVAL_WEIGHT: { 141 if (tmode == CEED_TRANSPOSE) 142 return CeedError(ceed, 1, 143 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 144 CeedInt Q = Q1d; 145 CeedScalar *qweight1d; 146 ierr = CeedBasisGetQWeights(basis, &qweight1d); CeedChk(ierr); 147 for (CeedInt d=0; d<dim; d++) { 148 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d); 149 for (CeedInt i=0; i<pre; i++) 150 for (CeedInt j=0; j<Q; j++) 151 for (CeedInt k=0; k<post; k++) 152 v[(i*Q + j)*post + k] = qweight1d[j] 153 * (d == 0 ? 1 : v[(i*Q + j)*post + k]); 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 case CEED_EVAL_INTERP: { 171 CeedInt P = ndof, Q = nqpt; 172 CeedScalar *interp; 173 ierr = CeedBasisGetInterp(basis, &interp); CeedChk(ierr); 174 if (tmode == CEED_TRANSPOSE) { 175 P = nqpt; Q = ndof; 176 } 177 ierr = CeedTensorContract_Ref(ceed, ncomp, P, 1, Q, interp, 178 tmode, add, u, v); 179 CeedChk(ierr); 180 } 181 break; 182 case CEED_EVAL_GRAD: { 183 CeedInt P = ndof, Q = dim*nqpt; 184 CeedScalar *grad; 185 ierr = CeedBasisGetGrad(basis, &grad); CeedChk(ierr); 186 if (tmode == CEED_TRANSPOSE) { 187 P = dim*nqpt; Q = ndof; 188 } 189 ierr = CeedTensorContract_Ref(ceed, ncomp, P, 1, Q, grad, 190 tmode, add, u, v); 191 CeedChk(ierr); 192 } 193 break; 194 case CEED_EVAL_WEIGHT: { 195 if (tmode == CEED_TRANSPOSE) 196 return CeedError(ceed, 1, 197 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 198 CeedScalar *qweight; 199 ierr = CeedBasisGetQWeights(basis, &qweight); CeedChk(ierr); 200 for (CeedInt i=0; i<nqpt; i++) 201 v[i] = qweight[i]; 202 } break; 203 case CEED_EVAL_DIV: 204 return CeedError(ceed, 1, "CEED_EVAL_DIV not supported"); 205 case CEED_EVAL_CURL: 206 return CeedError(ceed, 1, "CEED_EVAL_CURL not supported"); 207 case CEED_EVAL_NONE: 208 return CeedError(ceed, 1, 209 "CEED_EVAL_NONE does not make sense in this context"); 210 } 211 } 212 if (U) { 213 ierr = CeedVectorRestoreArrayRead(U, &u); CeedChk(ierr); 214 } 215 ierr = CeedVectorRestoreArray(V, &v); CeedChk(ierr); 216 return 0; 217 } 218 219 static int CeedBasisDestroy_Ref(CeedBasis basis) { 220 return 0; 221 } 222 223 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P1d, 224 CeedInt Q1d, const CeedScalar *interp1d, 225 const CeedScalar *grad1d, 226 const CeedScalar *qref1d, 227 const CeedScalar *qweight1d, 228 CeedBasis basis) { 229 int ierr; 230 Ceed ceed; 231 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 232 233 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply", 234 CeedBasisApply_Ref); CeedChk(ierr); 235 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", 236 CeedBasisDestroy_Ref); CeedChk(ierr); 237 return 0; 238 } 239 240 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, 241 CeedInt ndof, CeedInt nqpts, 242 const CeedScalar *interp, 243 const CeedScalar *grad, 244 const CeedScalar *qref, 245 const CeedScalar *qweight, 246 CeedBasis basis) { 247 int ierr; 248 Ceed ceed; 249 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 250 251 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Apply", 252 CeedBasisApply_Ref); CeedChk(ierr); 253 ierr = CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", 254 CeedBasisDestroy_Ref); CeedChk(ierr); 255 return 0; 256 } 257