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