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 using PETSc 19 20 #ifndef bp1_h 21 #define bp1_h 22 23 #ifndef __CUDACC__ 24 # include <math.h> 25 #endif 26 27 // ***************************************************************************** 28 // This QFunction sets up the geometric factors required to apply the 29 // mass operator 30 // 31 // The quadrature data is stored in the array qdata. 32 // 33 // We require the determinant of the Jacobian to properly compute integrals of 34 // the form: int( u v ) 35 // 36 // Qdata: detJ * w 37 // 38 // ***************************************************************************** 39 40 // ----------------------------------------------------------------------------- 41 CEED_QFUNCTION(SetupMassGeo)(void *ctx, const CeedInt Q, 42 const CeedScalar *const *in, 43 CeedScalar *const *out) { 44 const CeedScalar *J = in[0], *w = in[1]; 45 CeedScalar *qdata = out[0]; 46 47 // Quadrature Point Loop 48 CeedPragmaSIMD 49 for (CeedInt i=0; i<Q; i++) { 50 const CeedScalar detJ = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) - 51 J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) + 52 J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])); 53 qdata[i] = detJ * w[i]; 54 } // End of Quadrature Point Loop 55 return 0; 56 } 57 58 // ***************************************************************************** 59 // This QFunction sets up the rhs and true solution for the problem 60 // ***************************************************************************** 61 62 // ----------------------------------------------------------------------------- 63 CEED_QFUNCTION(SetupMassRhs)(void *ctx, const CeedInt Q, 64 const CeedScalar *const *in, 65 CeedScalar *const *out) { 66 const CeedScalar *x = in[0], *J = in[1], *w = in[2]; 67 CeedScalar *true_soln = out[0], *rhs = out[1]; 68 69 // Quadrature Point Loop 70 CeedPragmaSIMD 71 for (CeedInt i=0; i<Q; i++) { 72 const CeedScalar det = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) - 73 J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) + 74 J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])); 75 76 true_soln[i] = sqrt(x[i]*x[i] + x[i+Q]*x[i+Q] + x[i+2*Q]*x[i+2*Q]); 77 78 rhs[i] = det * w[i] * true_soln[i]; 79 } // End of Quadrature Point Loop 80 return 0; 81 } 82 83 84 // ***************************************************************************** 85 // This QFunction applies the mass operator for a scalar field. 86 // 87 // Inputs: 88 // u - Input vector at quadrature points 89 // qdata - Geometric factors 90 // 91 // Output: 92 // v - Output vector (test functions) at quadrature points 93 // 94 // ***************************************************************************** 95 96 // ----------------------------------------------------------------------------- 97 CEED_QFUNCTION(Mass)(void *ctx, const CeedInt Q, 98 const CeedScalar *const *in, CeedScalar *const *out) { 99 const CeedScalar *u = in[0], *qdata = in[1]; 100 CeedScalar *v = out[0]; 101 102 // Quadrature Point Loop 103 CeedPragmaSIMD 104 for (CeedInt i=0; i<Q; i++) 105 v[i] = qdata[i] * u[i]; 106 107 return 0; 108 } 109 // ----------------------------------------------------------------------------- 110 111 #endif // bp1_h 112