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 /// A structure used to pass additional data to f_build_diff 18 struct BuildContext { CeedInt dim, space_dim; }; 19 20 /// libCEED Q-function for building quadrature data for a diffusion operator 21 CEED_QFUNCTION(f_build_diff)(void *ctx, const CeedInt Q, 22 const CeedScalar *const *in, CeedScalar *const *out) { 23 struct BuildContext *bc = (struct BuildContext *)ctx; 24 // in[0] is Jacobians with shape [dim, nc=dim, Q] 25 // in[1] is quadrature weights, size (Q) 26 // 27 // At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store 28 // the symmetric part of the result. 29 const CeedScalar *J = in[0], *w = in[1]; 30 CeedScalar *qdata = out[0]; 31 32 switch (bc->dim + 10*bc->space_dim) { 33 case 11: 34 CeedPragmaSIMD 35 for (CeedInt i=0; i<Q; i++) { 36 qdata[i] = w[i] / J[i]; 37 } // End of Quadrature Point Loop 38 break; 39 case 22: 40 CeedPragmaSIMD 41 for (CeedInt i=0; i<Q; i++) { 42 // J: 0 2 qdata: 0 2 adj(J): J22 -J12 43 // 1 3 2 1 -J21 J11 44 const CeedScalar J11 = J[i+Q*0]; 45 const CeedScalar J21 = J[i+Q*1]; 46 const CeedScalar J12 = J[i+Q*2]; 47 const CeedScalar J22 = J[i+Q*3]; 48 const CeedScalar qw = w[i] / (J11*J22 - J21*J12); 49 qdata[i+Q*0] = qw * (J12*J12 + J22*J22); 50 qdata[i+Q*1] = qw * (J11*J11 + J21*J21); 51 qdata[i+Q*2] = - qw * (J11*J12 + J21*J22); 52 } // End of Quadrature Point Loop 53 break; 54 case 33: 55 CeedPragmaSIMD 56 for (CeedInt i=0; i<Q; i++) { 57 // Compute the adjoint 58 CeedScalar A[3][3]; 59 for (CeedInt j=0; j<3; j++) 60 for (CeedInt k=0; k<3; k++) 61 // Equivalent code with J as a VLA and no mod operations: 62 // A[k][j] = J[j+1][k+1]*J[j+2][k+2] - J[j+1][k+2]*J[j+2][k+1] 63 A[k][j] = J[i+Q*((j+1)%3+3*((k+1)%3))]*J[i+Q*((j+2)%3+3*((k+2)%3))] - 64 J[i+Q*((j+1)%3+3*((k+2)%3))]*J[i+Q*((j+2)%3+3*((k+1)%3))]; 65 66 // Compute quadrature weight / det(J) 67 const CeedScalar qw = w[i] / (J[i+Q*0]*A[0][0] + J[i+Q*1]*A[1][1] + 68 J[i+Q*2]*A[2][2]); 69 70 // Compute geometric factors 71 // Stored in Voigt convention 72 // 0 5 4 73 // 5 1 3 74 // 4 3 2 75 qdata[i+Q*0] = qw * (A[0][0]*A[0][0] + A[0][1]*A[0][1] + A[0][2]*A[0][2]); 76 qdata[i+Q*1] = qw * (A[1][0]*A[1][0] + A[1][1]*A[1][1] + A[1][2]*A[1][2]); 77 qdata[i+Q*2] = qw * (A[2][0]*A[2][0] + A[2][1]*A[2][1] + A[2][2]*A[2][2]); 78 qdata[i+Q*3] = qw * (A[1][0]*A[2][0] + A[1][1]*A[2][1] + A[1][2]*A[2][2]); 79 qdata[i+Q*4] = qw * (A[0][0]*A[2][0] + A[0][1]*A[2][1] + A[0][2]*A[2][2]); 80 qdata[i+Q*5] = qw * (A[0][0]*A[1][0] + A[0][1]*A[1][1] + A[0][2]*A[1][2]); 81 } // End of Quadrature Point Loop 82 break; 83 } 84 return 0; 85 } 86 87 /// libCEED Q-function for applying a diff operator 88 CEED_QFUNCTION(f_apply_diff)(void *ctx, const CeedInt Q, 89 const CeedScalar *const *in, CeedScalar *const *out) { 90 struct BuildContext *bc = (struct BuildContext *)ctx; 91 // in[0], out[0] have shape [dim, nc=1, Q] 92 const CeedScalar *ug = in[0], *qdata = in[1]; 93 CeedScalar *vg = out[0]; 94 95 switch (bc->dim) { 96 case 1: 97 CeedPragmaSIMD 98 for (CeedInt i=0; i<Q; i++) { 99 vg[i] = ug[i] * qdata[i]; 100 } // End of Quadrature Point Loop 101 break; 102 case 2: 103 CeedPragmaSIMD 104 for (CeedInt i=0; i<Q; i++) { 105 // Read spatial derivatives of u 106 const CeedScalar du[2] = {ug[i+Q*0], 107 ug[i+Q*1] 108 }; 109 110 // Read qdata (dXdxdXdxT symmetric matrix) 111 // Stored in Voigt convention 112 // 0 2 113 // 2 1 114 // *INDENT-OFF* 115 const CeedScalar dXdxdXdxT[2][2] = {{qdata[i+0*Q], 116 qdata[i+2*Q]}, 117 {qdata[i+2*Q], 118 qdata[i+1*Q]}}; 119 // *INDENT-ON* 120 // j = direction of vg 121 for (int j=0; j<2; j++) 122 vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] + 123 du[1] * dXdxdXdxT[1][j]); 124 } // End of Quadrature Point Loop 125 break; 126 case 3: 127 CeedPragmaSIMD 128 for (CeedInt i=0; i<Q; i++) { 129 // Read spatial derivatives of u 130 const CeedScalar du[3] = {ug[i+Q*0], 131 ug[i+Q*1], 132 ug[i+Q*2] 133 }; 134 135 // Read qdata (dXdxdXdxT symmetric matrix) 136 // Stored in Voigt convention 137 // 0 5 4 138 // 5 1 3 139 // 4 3 2 140 // *INDENT-OFF* 141 const CeedScalar dXdxdXdxT[3][3] = {{qdata[i+0*Q], 142 qdata[i+5*Q], 143 qdata[i+4*Q]}, 144 {qdata[i+5*Q], 145 qdata[i+1*Q], 146 qdata[i+3*Q]}, 147 {qdata[i+4*Q], 148 qdata[i+3*Q], 149 qdata[i+2*Q]} 150 }; 151 // *INDENT-ON* 152 // j = direction of vg 153 for (int j=0; j<3; j++) 154 vg[i+j*Q] = (du[0] * dXdxdXdxT[0][j] + 155 du[1] * dXdxdXdxT[1][j] + 156 du[2] * dXdxdXdxT[2][j]); 157 } // End of Quadrature Point Loop 158 break; 159 } 160 return 0; 161 } 162