1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// libCEED QFunctions for mass operator example for a vector field on the sphere using PETSc 19 20 #ifndef bp4sphere_h 21 #define bp4sphere_h 22 #include <ceed.h> 23 24 #ifndef __CUDACC__ 25 # include <math.h> 26 #endif 27 28 // ***************************************************************************** 29 // This QFunction sets up the rhs and true solution for the problem 30 // ***************************************************************************** 31 32 // ----------------------------------------------------------------------------- 33 CEED_QFUNCTION(SetupDiffRhs3)(void *ctx, const CeedInt Q, 34 const CeedScalar *const *in, 35 CeedScalar *const *out) { 36 // Inputs 37 const CeedScalar *X = in[0], *qdata = in[1]; 38 // Outputs 39 CeedScalar *true_soln = out[0], *rhs = out[1]; 40 41 // Context 42 const CeedScalar *context = (const CeedScalar*)ctx; 43 const CeedScalar R = context[0]; 44 45 // Quadrature Point Loop 46 CeedPragmaSIMD 47 for (CeedInt i=0; i<Q; i++) { 48 // Read global Cartesian coordinates 49 CeedScalar x = X[i+Q*0], y = X[i+Q*1], z = X[i+Q*2]; 50 // Normalize quadrature point coordinates to sphere 51 CeedScalar rad = sqrt(x*x + y*y + z*z); 52 x *= R / rad; 53 y *= R / rad; 54 z *= R / rad; 55 // Compute latitude and longitude 56 const CeedScalar theta = asin(z / R); // latitude 57 const CeedScalar lambda = atan2(y, x); // longitude 58 59 // Use absolute value of latitute for true solution 60 // Component 1 61 true_soln[i+0*Q] = sin(lambda) * cos(theta); 62 // Component 2 63 true_soln[i+1*Q] = 2 * true_soln[i+0*Q]; 64 // Component 3 65 true_soln[i+2*Q] = 3 * true_soln[i+0*Q]; 66 67 // Component 1 68 rhs[i+0*Q] = qdata[i+Q*0] * 2 * sin(lambda)*cos(theta) / (R*R); 69 // Component 2 70 rhs[i+1*Q] = 2 * rhs[i+0*Q]; 71 // Component 3 72 rhs[i+2*Q] = 3 * rhs[i+0*Q]; 73 } // End of Quadrature Point Loop 74 75 return 0; 76 } 77 78 // ***************************************************************************** 79 // This QFunction applies the diffusion operator for a vector field of 3 components. 80 // 81 // Inputs: 82 // ug - Input vector Jacobian at quadrature points 83 // qdata - Geometric factors 84 // 85 // Output: 86 // vJ - Output vector (test functions) Jacobian at quadrature points 87 // 88 // ***************************************************************************** 89 90 // ----------------------------------------------------------------------------- 91 CEED_QFUNCTION(Diff3)(void *ctx, const CeedInt Q, 92 const CeedScalar *const *in, CeedScalar *const *out) { 93 const CeedScalar *ug = in[0], *qdata = in[1]; 94 CeedScalar *vJ = out[0]; 95 96 // Quadrature Point Loop 97 CeedPragmaSIMD 98 for (CeedInt i=0; i<Q; i++) { 99 // Read spatial derivatives of u 100 const CeedScalar uJ[3][2] = {{ug[i+(0+0*3)*Q], 101 ug[i+(0+1*3)*Q]}, 102 {ug[i+(1+0*3)*Q], 103 ug[i+(1+1*3)*Q]}, 104 {ug[i+(2+0*3)*Q], 105 ug[i+(2+1*3)*Q]} 106 }; 107 // Read qdata 108 const CeedScalar wJ = qdata[i+Q*0]; 109 // -- Grad-to-Grad qdata 110 // ---- dXdx_j,k * dXdx_k,j 111 const CeedScalar dXdxdXdxT[2][2] = {{qdata[i+Q*1], 112 qdata[i+Q*3]}, 113 {qdata[i+Q*3], 114 qdata[i+Q*2]} 115 }; 116 117 for (int k=0; k<3; k++) // k = component 118 for (int j=0; j<2; j++) // j = direction of vg 119 vJ[i+(k+j*3)*Q] = wJ * (uJ[k][0] * dXdxdXdxT[0][j] + 120 uJ[k][1] * dXdxdXdxT[1][j]); 121 122 } // End of Quadrature Point Loop 123 124 return 0; 125 } 126 // ----------------------------------------------------------------------------- 127 128 #endif // bp4sphere_h 129