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