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