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