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