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 vector 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 rhs and true solution for the problem
18 // -----------------------------------------------------------------------------
SetupDiffRhs3(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)19 CEED_QFUNCTION(SetupDiffRhs3)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
20 // Inputs
21 const CeedScalar *X = in[0], *q_data = in[1];
22 // Outputs
23 CeedScalar *true_soln = out[0], *rhs = out[1];
24
25 // Context
26 const CeedScalar *context = (const CeedScalar *)ctx;
27 const CeedScalar R = context[0];
28
29 // Quadrature Point Loop
30 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
31 // Read global Cartesian coordinates
32 CeedScalar x = X[i + Q * 0], y = X[i + Q * 1], z = X[i + Q * 2];
33 // Normalize quadrature point coordinates to sphere
34 CeedScalar rad = sqrt(x * x + y * y + z * z);
35 x *= R / rad;
36 y *= R / rad;
37 z *= R / rad;
38 // Compute latitude and longitude
39 const CeedScalar theta = asin(z / R); // latitude
40 const CeedScalar lambda = atan2(y, x); // longitude
41
42 // Use absolute value of latitude for true solution
43 // Component 1
44 true_soln[i + 0 * Q] = sin(lambda) * cos(theta);
45 // Component 2
46 true_soln[i + 1 * Q] = 2 * true_soln[i + 0 * Q];
47 // Component 3
48 true_soln[i + 2 * Q] = 3 * true_soln[i + 0 * Q];
49
50 // Component 1
51 rhs[i + 0 * Q] = q_data[i + Q * 0] * 2 * sin(lambda) * cos(theta) / (R * R);
52 // Component 2
53 rhs[i + 1 * Q] = 2 * rhs[i + 0 * Q];
54 // Component 3
55 rhs[i + 2 * Q] = 3 * rhs[i + 0 * Q];
56 } // End of Quadrature Point Loop
57
58 return 0;
59 }
60
61 // -----------------------------------------------------------------------------
62 // This QFunction applies the diffusion operator for a vector field of 3 components.
63 //
64 // Inputs:
65 // ug - Input vector Jacobian at quadrature points
66 // q_data - Geometric factors
67 //
68 // Output:
69 // vJ - Output vector (test functions) Jacobian at quadrature points
70 // -----------------------------------------------------------------------------
Diff3(void * ctx,const CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)71 CEED_QFUNCTION(Diff3)(void *ctx, const CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
72 const CeedScalar *ug = in[0], *q_data = in[1];
73 CeedScalar *vJ = out[0];
74
75 // Quadrature Point Loop
76 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
77 // Read spatial derivatives of u
78 const CeedScalar uJ[3][2] = {
79 {ug[i + (0 + 0 * 3) * Q], ug[i + (0 + 1 * 3) * Q]},
80 {ug[i + (1 + 0 * 3) * Q], ug[i + (1 + 1 * 3) * Q]},
81 {ug[i + (2 + 0 * 3) * Q], ug[i + (2 + 1 * 3) * Q]}
82 };
83 // Read q_data
84 const CeedScalar w_det_J = q_data[i + Q * 0];
85 // -- Grad-to-Grad q_data
86 // ---- dXdx_j,k * dXdx_k,j
87 const CeedScalar dXdxdXdx_T[2][2] = {
88 {q_data[i + Q * 1], q_data[i + Q * 3]},
89 {q_data[i + Q * 3], q_data[i + Q * 2]}
90 };
91
92 for (int k = 0; k < 3; k++) { // k = component
93 for (int j = 0; j < 2; j++) { // j = direction of vg
94 vJ[i + (k + j * 3) * Q] = w_det_J * (uJ[k][0] * dXdxdXdx_T[0][j] + uJ[k][1] * dXdxdXdx_T[1][j]);
95 }
96 }
97 } // End of Quadrature Point Loop
98
99 return 0;
100 }
101 // -----------------------------------------------------------------------------
102