xref: /libCEED/examples/fluids/qfunctions/velocity_gradient_projection.h (revision 5cd6c1fb67d52eb6a42b887bb79c183682dd86ca)
1 // Copyright (c) 2017-2023, 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 #ifndef velocity_gradient_projection_h
9 #define velocity_gradient_projection_h
10 
11 #include <ceed.h>
12 
13 #include "newtonian_state.h"
14 #include "newtonian_types.h"
15 #include "utils.h"
16 
17 CEED_QFUNCTION_HELPER int VelocityGradientProjectionRHS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out,
18                                                         StateVariable state_var) {
19   const CeedScalar(*q)[CEED_Q_VLA]         = (const CeedScalar(*)[CEED_Q_VLA])in[0];
20   const CeedScalar(*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1];
21   const CeedScalar(*q_data)                = in[2];
22   const CeedScalar(*x)[CEED_Q_VLA]         = (const CeedScalar(*)[CEED_Q_VLA])in[3];
23   CeedScalar(*v)[CEED_Q_VLA]               = (CeedScalar(*)[CEED_Q_VLA])out[0];
24 
25   NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx;
26 
27   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
28     const CeedScalar qi[5]  = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
29     const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]};
30     CeedScalar       wdetJ, dXdx[3][3];
31     QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
32 
33     const State s = StateFromQ(context, qi, x_i, state_var);
34     State       grad_s[3];
35     for (CeedInt k = 0; k < 3; k++) {
36       CeedScalar dx_i[3] = {0}, dqi[5];
37       for (CeedInt j = 0; j < 5; j++) {
38         dqi[j] = Grad_q[0][j][i] * dXdx[0][k] + Grad_q[1][j][i] * dXdx[1][k] + Grad_q[2][j][i] * dXdx[2][k];
39       }
40       dx_i[k]   = 1.;
41       grad_s[k] = StateFromQ_fwd(context, s, dqi, x_i, dx_i, state_var);
42     }
43 
44     CeedScalar grad_velocity[3][3];
45     VelocityGradient(grad_s, grad_velocity);
46 
47     for (CeedInt j = 0; j < 3; j++) {
48       for (CeedInt k = 0; k < 3; k++) {
49         v[j * 3 + k][i] = wdetJ * grad_velocity[j][k];
50       }
51     }
52   }
53   return 0;
54 }
55 
56 CEED_QFUNCTION(VelocityGradientProjectionRHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
57   return VelocityGradientProjectionRHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE);
58 }
59 
60 CEED_QFUNCTION(VelocityGradientProjectionRHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
61   return VelocityGradientProjectionRHS(ctx, Q, in, out, STATEVAR_PRIMITIVE);
62 }
63 #endif  // velocity_gradient_projection_h
64