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)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])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 const CeedScalar wdetJ = q_data[0][i]; 31 const CeedScalar dXdx[3][3] = { 32 {q_data[1][i], q_data[2][i], q_data[3][i]}, 33 {q_data[4][i], q_data[5][i], q_data[6][i]}, 34 {q_data[7][i], q_data[8][i], q_data[9][i]} 35 }; 36 37 const State s = StateFromQ(context, qi, x_i, state_var); 38 State grad_s[3]; 39 for (CeedInt k = 0; k < 3; k++) { 40 CeedScalar dx_i[3] = {0}, dqi[5]; 41 for (CeedInt j = 0; j < 5; j++) { 42 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]; 43 } 44 dx_i[k] = 1.; 45 grad_s[k] = StateFromQ_fwd(context, s, dqi, x_i, dx_i, state_var); 46 } 47 48 CeedScalar grad_velocity[3][3]; 49 VelocityGradient(grad_s, grad_velocity); 50 51 for (CeedInt j = 0; j < 3; j++) { 52 for (CeedInt k = 0; k < 3; k++) { 53 v[j * 3 + k][i] = wdetJ * grad_velocity[j][k]; 54 } 55 } 56 } 57 return 0; 58 } 59 60 CEED_QFUNCTION(VelocityGradientProjectionRHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 61 return VelocityGradientProjectionRHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 62 } 63 64 CEED_QFUNCTION(VelocityGradientProjectionRHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 65 return VelocityGradientProjectionRHS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 66 } 67 #endif // velocity_gradient_projection_h 68