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 /// @file 9 /// QFunctions for the `bc_slip` boundary conditions 10 11 #include "bc_freestream_type.h" 12 #include "newtonian_state.h" 13 #include "newtonian_types.h" 14 #include "riemann_solver.h" 15 16 CEED_QFUNCTION_HELPER int Slip(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 17 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 18 const CeedScalar(*q_data_sur) = in[2]; 19 20 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 21 CeedScalar(*jac_data_sur) = out[1]; 22 23 const NewtonianIdealGasContext newt_ctx = (const NewtonianIdealGasContext)ctx; 24 const CeedScalar zeros[6] = {0.}; 25 26 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 27 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 28 State s = StateFromQ(newt_ctx, qi, state_var); 29 30 CeedScalar wdetJb, norm[3]; 31 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 32 wdetJb *= newt_ctx->is_implicit ? -1. : 1.; 33 34 CeedScalar vel_reflect[3]; 35 const CeedScalar vel_normal = Dot3(s.Y.velocity, norm); 36 for (CeedInt j = 0; j < 3; j++) vel_reflect[j] = s.Y.velocity[j] - 2. * norm[j] * vel_normal; 37 const CeedScalar Y_reflect[5] = {s.Y.pressure, vel_reflect[0], vel_reflect[1], vel_reflect[2], s.Y.temperature}; 38 State s_reflect = StateFromY(newt_ctx, Y_reflect); 39 40 StateConservative flux = RiemannFlux_HLLC(newt_ctx, s, s_reflect, norm); 41 42 CeedScalar Flux[5]; 43 UnpackState_U(flux, Flux); 44 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * Flux[j]; 45 46 StoredValuesPack(Q, i, 0, 5, qi, jac_data_sur); 47 StoredValuesPack(Q, i, 5, 6, zeros, jac_data_sur); // Every output value must be set 48 } 49 return 0; 50 } 51 52 CEED_QFUNCTION(Slip_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 53 return Slip(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 54 } 55 56 CEED_QFUNCTION(Slip_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 57 return Slip(ctx, Q, in, out, STATEVAR_PRIMITIVE); 58 } 59 60 CEED_QFUNCTION_HELPER int Slip_Jacobian(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 61 const CeedScalar(*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 62 const CeedScalar(*q_data_sur) = in[2]; 63 const CeedScalar(*jac_data_sur) = in[4]; 64 65 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 66 67 const NewtonianIdealGasContext newt_ctx = (const NewtonianIdealGasContext)ctx; 68 69 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 70 CeedScalar wdetJb, norm[3]; 71 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 72 wdetJb *= newt_ctx->is_implicit ? -1. : 1.; 73 74 CeedScalar qi[5], dqi[5]; 75 StoredValuesUnpack(Q, i, 0, 5, jac_data_sur, qi); 76 for (int j = 0; j < 5; j++) dqi[j] = dq[j][i]; 77 State s = StateFromQ(newt_ctx, qi, state_var); 78 State ds = StateFromQ_fwd(newt_ctx, s, dqi, state_var); 79 80 CeedScalar vel_reflect[3]; 81 const CeedScalar vel_normal = Dot3(s.Y.velocity, norm); 82 for (CeedInt j = 0; j < 3; j++) vel_reflect[j] = s.Y.velocity[j] - 2. * norm[j] * vel_normal; 83 const CeedScalar Y_reflect[5] = {s.Y.pressure, vel_reflect[0], vel_reflect[1], vel_reflect[2], s.Y.temperature}; 84 State s_reflect = StateFromY(newt_ctx, Y_reflect); 85 86 CeedScalar dvel_reflect[3]; 87 const CeedScalar dvel_normal = Dot3(ds.Y.velocity, norm); 88 for (CeedInt j = 0; j < 3; j++) dvel_reflect[j] = ds.Y.velocity[j] - 2. * norm[j] * dvel_normal; 89 const CeedScalar dY_reflect[5] = {ds.Y.pressure, dvel_reflect[0], dvel_reflect[1], dvel_reflect[2], ds.Y.temperature}; 90 State ds_reflect = StateFromY_fwd(newt_ctx, s_reflect, dY_reflect); 91 92 StateConservative dflux = RiemannFlux_HLLC_fwd(newt_ctx, s, ds, s_reflect, ds_reflect, norm); 93 94 CeedScalar dFlux[5]; 95 UnpackState_U(dflux, dFlux); 96 for (CeedInt j = 0; j < 5; j++) v[j][i] = -wdetJb * dFlux[j]; 97 } 98 return 0; 99 } 100 101 CEED_QFUNCTION(Slip_Jacobian_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 102 return Slip_Jacobian(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 103 } 104 105 CEED_QFUNCTION(Slip_Jacobian_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 106 return Slip_Jacobian(ctx, Q, in, out, STATEVAR_PRIMITIVE); 107 } 108