1 // Copyright (c) 2017-2026, 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 // -----------------------------------------------------------------------------
SetupDiffGeo(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)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 // -----------------------------------------------------------------------------
SetupDiffRhs(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)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 return 0;
93 }
94
95 // -----------------------------------------------------------------------------
96 // This QFunction applies the diffusion operator for a scalar field.
97 //
98 // Inputs:
99 // ug - Input vector gradient at quadrature points
100 // q_data - Geometric factors
101 //
102 // Output:
103 // vg - Output vector (test functions) gradient at quadrature points
104 // -----------------------------------------------------------------------------
Diff(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)105 CEED_QFUNCTION(Diff)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
106 const CeedScalar *ug = in[0], *q_data = in[1];
107 CeedScalar *vg = out[0];
108
109 // Quadrature Point Loop
110 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
111 // Read spatial derivatives of u
112 const CeedScalar du[3] = {ug[i + Q * 0], ug[i + Q * 1], ug[i + Q * 2]};
113 // Read q_data (dXdxdXdx_T symmetric matrix)
114 const CeedScalar dXdxdXdx_T[3][3] = {
115 {q_data[i + 1 * Q], q_data[i + 2 * Q], q_data[i + 3 * Q]},
116 {q_data[i + 2 * Q], q_data[i + 4 * Q], q_data[i + 5 * Q]},
117 {q_data[i + 3 * Q], q_data[i + 5 * Q], q_data[i + 6 * Q]}
118 };
119
120 for (int j = 0; j < 3; j++) { // j = direction of vg
121 vg[i + j * Q] = (du[0] * dXdxdXdx_T[0][j] + du[1] * dXdxdXdx_T[1][j] + du[2] * dXdxdXdx_T[2][j]);
122 }
123 } // End of Quadrature Point Loop
124 return 0;
125 }
126 // -----------------------------------------------------------------------------
127