1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed/types.h> 9 10 /// A structure used to pass additional data to f_build_mass_diff 11 struct BuildContext { 12 CeedInt dim, space_dim; 13 }; 14 15 /// libCEED Q-function for building quadrature data for a mass + diffusion operator 16 CEED_QFUNCTION(build_mass_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 17 struct BuildContext *build_data = (struct BuildContext *)ctx; 18 // in[0] is Jacobians with shape [dim, nc=dim, Q] 19 // in[1] is quadrature weights, size (Q) 20 // 21 // At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store 22 // the symmetric part of the result. 23 const CeedScalar *J = in[0], *w = in[1]; 24 CeedScalar *q_data = out[0]; 25 26 switch (build_data->dim + 10 * build_data->space_dim) { 27 case 11: 28 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 29 // Mass 30 q_data[i + Q * 0] = w[i] * J[i]; 31 // Diffusion 32 q_data[i + Q * 1] = w[i] / J[i]; 33 } // End of Quadrature Point Loop 34 break; 35 case 22: 36 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 37 // J: 0 2 q_data: 1 3 adj(J): J22 -J12 38 // 1 3 3 2 -J21 J11 39 const CeedScalar J11 = J[i + Q * 0]; 40 const CeedScalar J21 = J[i + Q * 1]; 41 const CeedScalar J12 = J[i + Q * 2]; 42 const CeedScalar J22 = J[i + Q * 3]; 43 const CeedScalar qw = w[i] / (J11 * J22 - J21 * J12); 44 45 // Mass 46 q_data[i + Q * 0] = w[i] * (J11 * J22 - J21 * J12); 47 // Diffusion 48 q_data[i + Q * 1] = qw * (J12 * J12 + J22 * J22); 49 q_data[i + Q * 2] = qw * (J11 * J11 + J21 * J21); 50 q_data[i + Q * 3] = -qw * (J11 * J12 + J21 * J22); 51 } // End of Quadrature Point Loop 52 break; 53 case 33: 54 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 55 // Compute the adjoint 56 CeedScalar A[3][3]; 57 for (CeedInt j = 0; j < 3; j++) { 58 for (CeedInt k = 0; k < 3; k++) { 59 // Equivalent code with J as a VLA and no mod operations: 60 // A[k][j] = J[j+1][k+1]*J[j+2][k+2] - J[j+1][k+2]*J[j+2][k+1] 61 A[k][j] = J[i + Q * ((j + 1) % 3 + 3 * ((k + 1) % 3))] * J[i + Q * ((j + 2) % 3 + 3 * ((k + 2) % 3))] - 62 J[i + Q * ((j + 1) % 3 + 3 * ((k + 2) % 3))] * J[i + Q * ((j + 2) % 3 + 3 * ((k + 1) % 3))]; 63 } 64 } 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[0][1] + J[i + Q * 2] * A[0][2]); 68 69 // Mass 70 q_data[i + Q * 0] = w[i] * (J[i + Q * 0] * A[0][0] + J[i + Q * 1] * A[0][1] + J[i + Q * 2] * A[0][2]); 71 // Diffusion 72 // Stored in Voigt convention 73 // 1 6 5 74 // 6 2 4 75 // 5 4 3 76 q_data[i + Q * 1] = qw * (A[0][0] * A[0][0] + A[0][1] * A[0][1] + A[0][2] * A[0][2]); 77 q_data[i + Q * 2] = qw * (A[1][0] * A[1][0] + A[1][1] * A[1][1] + A[1][2] * A[1][2]); 78 q_data[i + Q * 3] = qw * (A[2][0] * A[2][0] + A[2][1] * A[2][1] + A[2][2] * A[2][2]); 79 q_data[i + Q * 4] = qw * (A[1][0] * A[2][0] + A[1][1] * A[2][1] + A[1][2] * A[2][2]); 80 q_data[i + Q * 5] = qw * (A[0][0] * A[2][0] + A[0][1] * A[2][1] + A[0][2] * A[2][2]); 81 q_data[i + Q * 6] = qw * (A[0][0] * A[1][0] + A[0][1] * A[1][1] + A[0][2] * A[1][2]); 82 } // End of Quadrature Point Loop 83 break; 84 } 85 return CEED_ERROR_SUCCESS; 86 } 87 88 /// libCEED Q-function for applying a mass + diffusion operator 89 CEED_QFUNCTION(apply_mass_diff)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 90 struct BuildContext *build_data = (struct BuildContext *)ctx; 91 // in[0], out[0] have shape [1, nc=1, Q] 92 // in[1], out[1] have shape [dim, nc=1, Q] 93 const CeedScalar *u = in[0], *ug = in[1], *q_data = in[2]; 94 CeedScalar *v = out[0], *vg = out[1]; 95 96 switch (build_data->dim) { 97 case 1: 98 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 99 // Mass 100 v[i] = q_data[i + Q * 0] * u[i]; 101 // Diffusion 102 vg[i] = ug[i] * q_data[i + Q * 1]; 103 } // End of Quadrature Point Loop 104 break; 105 case 2: 106 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 107 // Mass 108 v[i] = q_data[i + Q * 0] * u[i]; 109 110 // Diffusion 111 // Read spatial derivatives of u 112 const CeedScalar du[2] = {ug[i + Q * 0], ug[i + Q * 1]}; 113 114 // Read q_data (dXdxdXdx_T symmetric matrix) 115 // Stored in Voigt convention 116 // 1 3 117 // 3 2 118 const CeedScalar dXdxdXdx_T[2][2] = { 119 {q_data[i + 1 * Q], q_data[i + 3 * Q]}, 120 {q_data[i + 3 * Q], q_data[i + 2 * Q]} 121 }; 122 // j = direction of vg 123 for (int j = 0; j < 2; j++) vg[i + j * Q] = (du[0] * dXdxdXdx_T[0][j] + du[1] * dXdxdXdx_T[1][j]); 124 } // End of Quadrature Point Loop 125 break; 126 case 3: 127 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 128 // Mass 129 v[i] = q_data[i + Q * 0] * u[i]; 130 131 // Diffusion 132 // Read spatial derivatives of u 133 const CeedScalar du[3] = {ug[i + Q * 0], ug[i + Q * 1], ug[i + Q * 2]}; 134 135 // Read q_data (dXdxdXdx_T symmetric matrix) 136 // Stored in Voigt convention 137 // 0 5 4 138 // 5 1 3 139 // 4 3 2 140 const CeedScalar dXdxdXdx_T[3][3] = { 141 {q_data[i + 1 * Q], q_data[i + 6 * Q], q_data[i + 5 * Q]}, 142 {q_data[i + 6 * Q], q_data[i + 2 * Q], q_data[i + 4 * Q]}, 143 {q_data[i + 5 * Q], q_data[i + 4 * Q], q_data[i + 3 * Q]} 144 }; 145 // j = direction of vg 146 for (int j = 0; j < 3; j++) vg[i + j * Q] = (du[0] * dXdxdXdx_T[0][j] + du[1] * dXdxdXdx_T[1][j] + du[2] * dXdxdXdx_T[2][j]); 147 } // End of Quadrature Point Loop 148 break; 149 } 150 return CEED_ERROR_SUCCESS; 151 } 152