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