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 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 // ----------------------------------------------------------------------------- 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 // ----------------------------------------------------------------------------- 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