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