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