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 mass operator example for a scalar field on the sphere 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 for integration and coordinate transformations when reference coordinates have a different
18 // dimension than the one of physical coordinates
19 //
20 // Reference (parent) 2D coordinates: X \in [-1, 1]^2
21 //
22 // Global 3D physical coordinates given by the mesh: xx \in [-R, R]^3 with R radius of the sphere
23 //
24 // Local 3D physical coordinates on the 2D manifold: x \in [-l, l]^3 with l half edge of the cube inscribed in the sphere
25 //
26 // Change of coordinates matrix computed by the library:
27 // (physical 3D coords relative to reference 2D coords)
28 // dxx_j/dX_i (indicial notation) [3 * 2]
29 //
30 // Change of coordinates x (on the 2D manifold) relative to xx (phyisical 3D):
31 // dx_i/dxx_j (indicial notation) [3 * 3]
32 //
33 // Change of coordinates x (on the 2D manifold) relative to X (reference 2D):
34 // (by chain rule)
35 // dx_i/dX_j [3 * 2] = dx_i/dxx_k [3 * 3] * dxx_k/dX_j [3 * 2]
36 //
37 // mod_J is given by the magnitude of the cross product of the columns of dx_i/dX_j
38 //
39 // The quadrature data is stored in the array q_data.
40 //
41 // We require the determinant of the Jacobian to properly compute integrals of the form: int( u v )
42 //
43 // Qdata: mod_J * w
44 // -----------------------------------------------------------------------------
SetupMassGeo(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)45 CEED_QFUNCTION(SetupMassGeo)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
46 // Inputs
47 const CeedScalar *X = in[0], *J = in[1], *w = in[2];
48 // Outputs
49 CeedScalar *q_data = out[0];
50
51 // Quadrature Point Loop
52 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
53 // Read global Cartesian coordinates
54 const CeedScalar xx[3] = {X[i + 0 * Q], X[i + 1 * Q], X[i + 2 * Q]};
55
56 // Read dxxdX Jacobian entries, stored as
57 // 0 3
58 // 1 4
59 // 2 5
60 const CeedScalar dxxdX[3][2] = {
61 {J[i + Q * 0], J[i + Q * 3]},
62 {J[i + Q * 1], J[i + Q * 4]},
63 {J[i + Q * 2], J[i + Q * 5]}
64 };
65
66 // Setup
67 // x = xx (xx^T xx)^{-1/2}
68 // dx/dxx = I (xx^T xx)^{-1/2} - xx xx^T (xx^T xx)^{-3/2}
69 const CeedScalar mod_xx_sq = xx[0] * xx[0] + xx[1] * xx[1] + xx[2] * xx[2];
70 CeedScalar xx_sq[3][3];
71 for (int j = 0; j < 3; j++) {
72 for (int k = 0; k < 3; k++) xx_sq[j][k] = xx[j] * xx[k] / (sqrt(mod_xx_sq) * mod_xx_sq);
73 }
74
75 const CeedScalar dxdxx[3][3] = {
76 {1. / sqrt(mod_xx_sq) - xx_sq[0][0], -xx_sq[0][1], -xx_sq[0][2] },
77 {-xx_sq[1][0], 1. / sqrt(mod_xx_sq) - xx_sq[1][1], -xx_sq[1][2] },
78 {-xx_sq[2][0], -xx_sq[2][1], 1. / sqrt(mod_xx_sq) - xx_sq[2][2]}
79 };
80
81 CeedScalar dxdX[3][2];
82 for (int j = 0; j < 3; j++) {
83 for (int k = 0; k < 2; k++) {
84 dxdX[j][k] = 0;
85 for (int l = 0; l < 3; l++) dxdX[j][k] += dxdxx[j][l] * dxxdX[l][k];
86 }
87 }
88
89 // J is given by the cross product of the columns of dxdX
90 const CeedScalar J[3] = {dxdX[1][0] * dxdX[2][1] - dxdX[2][0] * dxdX[1][1], dxdX[2][0] * dxdX[0][1] - dxdX[0][0] * dxdX[2][1],
91 dxdX[0][0] * dxdX[1][1] - dxdX[1][0] * dxdX[0][1]};
92
93 // Use the magnitude of J as our detJ (volume scaling factor)
94 const CeedScalar mod_J = sqrt(J[0] * J[0] + J[1] * J[1] + J[2] * J[2]);
95
96 // Interp-to-Interp q_data
97 q_data[i + Q * 0] = mod_J * w[i];
98 } // End of Quadrature Point Loop
99
100 return 0;
101 }
102
103 // -----------------------------------------------------------------------------
104 // This QFunction sets up the rhs and true solution for the problem
105 // -----------------------------------------------------------------------------
SetupMassRhs(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)106 CEED_QFUNCTION(SetupMassRhs)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
107 // Inputs
108 const CeedScalar *X = in[0], *q_data = in[1];
109 // Outputs
110 CeedScalar *true_soln = out[0], *rhs = out[1];
111
112 // Context
113 const CeedScalar *context = (const CeedScalar *)ctx;
114 const CeedScalar R = context[0];
115
116 // Quadrature Point Loop
117 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
118 // Compute latitude
119 const CeedScalar theta = asin(X[i + 2 * Q] / R);
120
121 // Use absolute value of latitude for true solution
122 true_soln[i] = fabs(theta);
123
124 rhs[i] = q_data[i] * true_soln[i];
125 } // End of Quadrature Point Loop
126
127 return 0;
128 }
129
130 // -----------------------------------------------------------------------------
131 // This QFunction applies the mass operator for a scalar field.
132 //
133 // Inputs:
134 // u - Input vector at quadrature points
135 // q_data - Geometric factors
136 //
137 // Output:
138 // v - Output vector (test functions) at quadrature points
139 // -----------------------------------------------------------------------------
Mass(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)140 CEED_QFUNCTION(Mass)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
141 // Inputs
142 const CeedScalar *u = in[0], *q_data = in[1];
143 // Outputs
144 CeedScalar *v = out[0];
145
146 // Quadrature Point Loop
147 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) v[i] = q_data[i] * u[i];
148
149 return 0;
150 }
151 // -----------------------------------------------------------------------------
152