1 // Copyright (c) 2017-2025, 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 /// @file 9 /// libCEED QFunctions for diffusion operator example using PETSc 10 11 #include <ceed/types.h> 12 #ifndef CEED_RUNNING_JIT_PASS 13 #include <math.h> 14 #endif 15 16 // ----------------------------------------------------------------------------- 17 // This QFunction sets up the geometric factors required to apply the diffusion operator 18 // 19 // We require the product of the inverse of the Jacobian and its transpose to properly compute integrals of the form: int( gradv gradu) 20 // 21 // Determinant of Jacobian: 22 // detJ = J11*A11 + J21*A12 + J31*A13 23 // Jij = Jacobian entry ij 24 // Aij = Adjoint ij 25 // 26 // Inverse of Jacobian: 27 // Bij = Aij / detJ 28 // 29 // Product of Inverse and Transpose: 30 // BBij = sum( Bik Bkj ) 31 // 32 // Stored: w B^T B detJ = w A^T A / detJ 33 // Note: This matrix is symmetric, so we only store 6 distinct entries 34 // qd: 1 4 7 35 // 2 5 8 36 // 3 6 9 37 // ----------------------------------------------------------------------------- 38 CEED_QFUNCTION(SetupDiffGeo)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 39 // Inputs 40 const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1]; 41 const CeedScalar(*w) = in[2]; // Note: *X = in[0] 42 // Outputs 43 CeedScalar(*qd) = out[0]; 44 45 const CeedInt dim = 3; 46 // Quadrature Point Loop 47 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 48 // Setup 49 CeedScalar A[3][3]; 50 for (CeedInt j = 0; j < dim; j++) { 51 for (CeedInt k = 0; k < dim; k++) { 52 // Equivalent code with no mod operations: 53 // A[k][j] = J[k+1][j+1]*J[k+2][j+2] - J[k+1][j+2]*J[k+2][j+1] 54 A[k][j] = J[(k + 1) % dim][(j + 1) % dim][i] * J[(k + 2) % dim][(j + 2) % dim][i] - 55 J[(k + 1) % dim][(j + 2) % dim][i] * J[(k + 2) % dim][(j + 1) % dim][i]; 56 } 57 } 58 const CeedScalar detJ = J[0][0][i] * A[0][0] + J[0][1][i] * A[0][1] + J[0][2][i] * A[0][2]; 59 60 const CeedScalar qw = w[i] / detJ; 61 qd[i + Q * 0] = w[i] * detJ; 62 qd[i + Q * 1] = qw * (A[0][0] * A[0][0] + A[0][1] * A[0][1] + A[0][2] * A[0][2]); 63 qd[i + Q * 2] = qw * (A[0][0] * A[1][0] + A[0][1] * A[1][1] + A[0][2] * A[1][2]); 64 qd[i + Q * 3] = qw * (A[0][0] * A[2][0] + A[0][1] * A[2][1] + A[0][2] * A[2][2]); 65 qd[i + Q * 4] = qw * (A[1][0] * A[1][0] + A[1][1] * A[1][1] + A[1][2] * A[1][2]); 66 qd[i + Q * 5] = qw * (A[1][0] * A[2][0] + A[1][1] * A[2][1] + A[1][2] * A[2][2]); 67 qd[i + Q * 6] = qw * (A[2][0] * A[2][0] + A[2][1] * A[2][1] + A[2][2] * A[2][2]); 68 } // End of Quadrature Point Loop 69 70 return 0; 71 } 72 73 // ----------------------------------------------------------------------------- 74 // This QFunction sets up the rhs and true solution for the problem 75 // ----------------------------------------------------------------------------- 76 CEED_QFUNCTION(SetupDiffRhs)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 77 #ifndef M_PI 78 #define M_PI 3.14159265358979323846 79 #endif 80 const CeedScalar *x = in[0], *w = in[1]; 81 CeedScalar *true_soln = out[0], *rhs = out[1]; 82 83 // Quadrature Point Loop 84 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 85 const CeedScalar c[3] = {0, 1., 2.}; 86 const CeedScalar k[3] = {1., 2., 3.}; 87 88 true_soln[i] = sin(M_PI * (c[0] + k[0] * x[i + Q * 0])) * sin(M_PI * (c[1] + k[1] * x[i + Q * 1])) * sin(M_PI * (c[2] + k[2] * x[i + Q * 2])); 89 90 rhs[i] = w[i + Q * 0] * M_PI * M_PI * (k[0] * k[0] + k[1] * k[1] + k[2] * k[2]) * true_soln[i]; 91 } // End of Quadrature Point Loop 92 93 return 0; 94 } 95 96 // ----------------------------------------------------------------------------- 97 // This QFunction applies the diffusion operator for a scalar field. 98 // 99 // Inputs: 100 // ug - Input vector gradient at quadrature points 101 // q_data - Geometric factors 102 // 103 // Output: 104 // vg - Output vector (test functions) gradient at quadrature points 105 // ----------------------------------------------------------------------------- 106 CEED_QFUNCTION(Diff)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 107 const CeedScalar *ug = in[0], *q_data = in[1]; 108 CeedScalar *vg = out[0]; 109 110 // Quadrature Point Loop 111 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 112 // Read spatial derivatives of u 113 const CeedScalar du[3] = {ug[i + Q * 0], ug[i + Q * 1], ug[i + Q * 2]}; 114 // Read q_data (dXdxdXdx_T symmetric matrix) 115 const CeedScalar dXdxdXdx_T[3][3] = { 116 {q_data[i + 1 * Q], q_data[i + 2 * Q], q_data[i + 3 * Q]}, 117 {q_data[i + 2 * Q], q_data[i + 4 * Q], q_data[i + 5 * Q]}, 118 {q_data[i + 3 * Q], q_data[i + 5 * Q], q_data[i + 6 * Q]} 119 }; 120 121 for (int j = 0; j < 3; j++) { // j = direction of vg 122 vg[i + j * Q] = (du[0] * dXdxdXdx_T[0][j] + du[1] * dXdxdXdx_T[1][j] + du[2] * dXdxdXdx_T[2][j]); 123 } 124 } // End of Quadrature Point Loop 125 return 0; 126 } 127 // ----------------------------------------------------------------------------- 128