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, CeedInt A, CeedInt B, CeedInt C, 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 const CeedInt dim = basis->dim; 53 const CeedInt ncomp = basis->ncomp; 54 CeedInt nqpt; 55 ierr = CeedBasisGetNumQuadraturePoints(basis, &nqpt); CeedChk(ierr); 56 CeedInt ndof; 57 ierr = CeedBasisGetNumNodes(basis, &ndof); CeedChk(ierr); 58 const CeedInt add = (tmode == CEED_TRANSPOSE); 59 60 if (nelem != 1) 61 return CeedError(basis->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 if (basis->tensorbasis) { 72 switch (emode) { 73 // Interpolate to/from quadrature points 74 case CEED_EVAL_INTERP: { 75 CeedInt P = basis->P1d, Q = basis->Q1d; 76 if (tmode == CEED_TRANSPOSE) { 77 P = basis->Q1d; Q = basis->P1d; 78 } 79 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1; 80 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 81 for (CeedInt d=0; d<dim; d++) { 82 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q, basis->interp1d, 83 tmode, add&&(d==dim-1), 84 d==0?u:tmp[d%2], d==dim-1?v:tmp[(d+1)%2]); 85 CeedChk(ierr); 86 pre /= P; 87 post *= Q; 88 } 89 } break; 90 // Evaluate the gradient to/from quadrature points 91 case CEED_EVAL_GRAD: { 92 CeedInt P = basis->P1d, Q = basis->Q1d; 93 // In CEED_NOTRANSPOSE mode: 94 // u has shape [dim, ncomp, P^dim, nelem], row-major layout 95 // v has shape [dim, ncomp, Q^dim, nelem], row-major layout 96 // In CEED_TRANSPOSE mode, the sizes of u and v are switched. 97 if (tmode == CEED_TRANSPOSE) { 98 P = basis->Q1d, Q = basis->P1d; 99 } 100 CeedScalar tmp[2][ncomp*Q*CeedIntPow(P>Q?P:Q, dim-1)]; 101 for (CeedInt p = 0; p < dim; p++) { 102 CeedInt pre = ncomp*CeedIntPow(P, dim-1), post = 1; 103 for (CeedInt d=0; d<dim; d++) { 104 ierr = CeedTensorContract_Ref(basis->ceed, pre, P, post, Q, 105 (p==d)?basis->grad1d:basis->interp1d, 106 tmode, add&&(d==dim-1), 107 (d == 0 108 ? (tmode==CEED_NOTRANSPOSE?u:u+p*ncomp*nqpt) 109 : tmp[d%2]), 110 (d == dim-1 111 ? (tmode==CEED_TRANSPOSE?v:v+p*ncomp*nqpt) 112 : tmp[(d+1)%2])); 113 CeedChk(ierr); 114 pre /= P; 115 post *= Q; 116 } 117 } 118 } break; 119 // Retrieve interpolation weights 120 case CEED_EVAL_WEIGHT: { 121 if (tmode == CEED_TRANSPOSE) 122 return CeedError(basis->ceed, 1, 123 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 124 CeedInt Q = basis->Q1d; 125 for (CeedInt d=0; d<dim; d++) { 126 CeedInt pre = CeedIntPow(Q, dim-d-1), post = CeedIntPow(Q, d); 127 for (CeedInt i=0; i<pre; i++) 128 for (CeedInt j=0; j<Q; j++) 129 for (CeedInt k=0; k<post; k++) 130 v[(i*Q + j)*post + k] = basis->qweight1d[j] 131 * (d == 0 ? 1 : v[(i*Q + j)*post + k]); 132 } 133 } break; 134 // Evaluate the divergence to/from the quadrature points 135 case CEED_EVAL_DIV: 136 return CeedError(basis->ceed, 1, "CEED_EVAL_DIV not supported"); 137 // Evaluate the curl to/from the quadrature points 138 case CEED_EVAL_CURL: 139 return CeedError(basis->ceed, 1, "CEED_EVAL_CURL not supported"); 140 // Take no action, BasisApply should not have been called 141 case CEED_EVAL_NONE: 142 return CeedError(basis->ceed, 1, 143 "CEED_EVAL_NONE does not make sense in this context"); 144 } 145 } else { 146 // Non-tensor basis 147 switch (emode) { 148 case CEED_EVAL_INTERP: { 149 CeedInt P = basis->P, Q = basis->Q; 150 if (tmode == CEED_TRANSPOSE) { 151 P = basis->Q; Q = basis->P; 152 } 153 ierr = CeedTensorContract_Ref(basis->ceed, ncomp, P, 1, Q, basis->interp1d, 154 tmode, add, u, v); 155 CeedChk(ierr); 156 } 157 break; 158 case CEED_EVAL_GRAD: { 159 CeedInt P = basis->P, Q = dim*basis->Q; 160 if (tmode == CEED_TRANSPOSE) { 161 P = dim*basis->Q; Q = basis->P; 162 } 163 ierr = CeedTensorContract_Ref(basis->ceed, ncomp, P, 1, Q, basis->grad1d, 164 tmode, add, u, v); 165 CeedChk(ierr); 166 } 167 break; 168 case CEED_EVAL_WEIGHT: { 169 if (tmode == CEED_TRANSPOSE) 170 return CeedError(basis->ceed, 1, 171 "CEED_EVAL_WEIGHT incompatible with CEED_TRANSPOSE"); 172 for (CeedInt i=0; i<nqpt; i++) 173 v[i] = basis->qweight1d[i]; 174 } break; 175 case CEED_EVAL_DIV: 176 return CeedError(basis->ceed, 1, "CEED_EVAL_DIV not supported"); 177 case CEED_EVAL_CURL: 178 return CeedError(basis->ceed, 1, "CEED_EVAL_CURL not supported"); 179 case CEED_EVAL_NONE: 180 return CeedError(basis->ceed, 1, 181 "CEED_EVAL_NONE does not make sense in this context"); 182 } 183 } 184 return 0; 185 } 186 187 static int CeedBasisDestroy_Ref(CeedBasis basis) { 188 return 0; 189 } 190 191 int CeedBasisCreateTensorH1_Ref(CeedInt dim, CeedInt P1d, 192 CeedInt Q1d, const CeedScalar *interp1d, 193 const CeedScalar *grad1d, 194 const CeedScalar *qref1d, 195 const CeedScalar *qweight1d, 196 CeedBasis basis) { 197 basis->Apply = CeedBasisApply_Ref; 198 basis->Destroy = CeedBasisDestroy_Ref; 199 return 0; 200 } 201 202 int CeedBasisCreateH1_Ref(CeedElemTopology topo, CeedInt dim, 203 CeedInt ndof, CeedInt nqpts, 204 const CeedScalar *interp, 205 const CeedScalar *grad, 206 const CeedScalar *qref, 207 const CeedScalar *qweight, 208 CeedBasis basis) { 209 basis->Apply = CeedBasisApply_Ref; 210 basis->Destroy = CeedBasisDestroy_Ref; 211 return 0; 212 } 213