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