1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Operator for HONEE
6 #include <ceed/types.h>
7
8 #include "newtonian_state.h"
9 #include "newtonian_types.h"
10 #include "utils.h"
11
12 typedef struct ChannelContext_ *ChannelContext;
13 struct ChannelContext_ {
14 bool implicit; // !< Using implicit timesteping or not
15 CeedScalar theta0; // !< Reference temperature
16 CeedScalar P0; // !< Reference Pressure
17 CeedScalar umax; // !< Centerline velocity
18 CeedScalar center; // !< Y Coordinate for center of channel
19 CeedScalar H; // !< Channel half-height
20 CeedScalar B; // !< Body-force driving the flow
21 struct NewtonianIdealGasContext_ newt_ctx;
22 };
23
Exact_Channel(CeedInt dim,CeedScalar time,const CeedScalar X[],CeedInt Nf,void * ctx)24 CEED_QFUNCTION_HELPER State Exact_Channel(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, void *ctx) {
25 const ChannelContext context = (ChannelContext)ctx;
26 const CeedScalar theta0 = context->theta0;
27 const CeedScalar P0 = context->P0;
28 const CeedScalar umax = context->umax;
29 const CeedScalar center = context->center;
30 const CeedScalar H = context->H;
31 NewtonianIGProperties gas = context->newt_ctx.gas;
32 const CeedScalar cp = gas.cp;
33 const CeedScalar mu = gas.mu;
34 const CeedScalar k = gas.k;
35 // There is a gravity body force but it is excluded from
36 // the potential energy due to periodicity.
37 // g = (g, 0, 0)
38 // x = (0, x_2, x_3)
39 // e_potential = dot(g, x) = 0
40 const CeedScalar x[3] = {0, X[1], X[2]};
41
42 const CeedScalar Pr = mu / (cp * k);
43 const CeedScalar Ec = (umax * umax) / (cp * theta0);
44 const CeedScalar theta = theta0 * (1 + (Pr * Ec / 3) * (1 - Square(Square((x[1] - center) / H))));
45 CeedScalar Y[5] = {0.};
46 Y[0] = P0;
47 Y[1] = umax * (1 - Square((x[1] - center) / H));
48 Y[2] = 0.;
49 Y[3] = 0.;
50 Y[4] = theta;
51
52 return StateFromY(gas, Y);
53 }
54
55 // *****************************************************************************
56 // This QFunction set the initial condition
57 // *****************************************************************************
ICsChannel(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)58 CEED_QFUNCTION(ICsChannel)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
59 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
60 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
61
62 const ChannelContext context = (ChannelContext)ctx;
63 NewtonianIGProperties gas = context->newt_ctx.gas;
64
65 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
66 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]};
67 State s = Exact_Channel(3, 0., x, 5, ctx);
68 CeedScalar q[5];
69 StateToQ(gas, s, q, context->newt_ctx.state_var);
70 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j];
71 }
72 return 0;
73 }
74
75 // *****************************************************************************
76 // This QFunction set the inflow boundary condition for conservative variables
77 // *****************************************************************************
Channel_Inflow(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)78 CEED_QFUNCTION(Channel_Inflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
79 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
80 const CeedScalar(*q_data_sur) = in[2];
81 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3];
82 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
83
84 const ChannelContext context = (ChannelContext)ctx;
85 const bool is_implicit = context->implicit;
86 NewtonianIGProperties gas = context->newt_ctx.gas;
87 const CeedScalar gamma = HeatCapacityRatio(gas);
88
89 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
90 CeedScalar wdetJb, norm[3];
91 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm);
92 wdetJb *= is_implicit ? -1. : 1.;
93
94 // There is a gravity body force but it is excluded from
95 // the potential energy due to periodicity.
96 // g = (g, 0, 0)
97 // x = (0, x_2, x_3)
98 // e_potential = dot(g, x) = 0
99 const CeedScalar x[3] = {0, X[1][i], X[2][i]};
100
101 // Calculate prescribed inflow values
102 State s_exact = Exact_Channel(3, 0., x, 5, ctx);
103 CeedScalar q_exact[5] = {0.};
104 UnpackState_U(s_exact.U, q_exact);
105
106 // Find pressure using state inside the domain
107 CeedScalar q_inside[5] = {0};
108 for (CeedInt j = 0; j < 5; j++) q_inside[j] = q[j][i];
109 State s_inside = StateFromU(gas, q_inside);
110 const CeedScalar P = s_inside.Y.pressure;
111
112 // Find inflow state using calculated P and prescribed velocity, theta0
113 const CeedScalar e_internal = gas.cv * s_exact.Y.temperature;
114 const CeedScalar rho_in = P / ((gamma - 1) * e_internal);
115 const CeedScalar E_kinetic = .5 * rho_in * Dot3(s_exact.Y.velocity, s_exact.Y.velocity);
116 const CeedScalar E = rho_in * e_internal + E_kinetic;
117
118 // The Physics
119 // Zero v so all future terms can safely sum into it
120 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.;
121
122 const CeedScalar u_normal = Dot3(norm, s_exact.Y.velocity);
123
124 // The Physics
125 // -- Density
126 v[0][i] -= wdetJb * rho_in * u_normal;
127
128 // -- Momentum
129 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho_in * u_normal * s_exact.Y.velocity[j] + norm[j] * P);
130
131 // -- Total Energy Density
132 v[4][i] -= wdetJb * u_normal * (E + P);
133 }
134 return 0;
135 }
136
137 // *****************************************************************************
138 // This QFunction set the outflow boundary condition for conservative variables
139 // *****************************************************************************
Channel_Outflow(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)140 CEED_QFUNCTION(Channel_Outflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
141 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
142 const CeedScalar(*q_data_sur) = in[2];
143 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
144
145 const ChannelContext context = (ChannelContext)ctx;
146 const bool is_implicit = context->implicit;
147
148 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
149 CeedScalar wdetJb, norm[3];
150 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm);
151 wdetJb *= is_implicit ? -1. : 1.;
152
153 const CeedScalar rho = q[0][i];
154 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho};
155 const CeedScalar E = q[4][i];
156
157 // The Physics
158 // Zero v so all future terms can safely sum into it
159 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.;
160
161 // Implementing outflow condition
162 const CeedScalar P = context->P0; // pressure
163 const CeedScalar u_normal = Dot3(norm, u); // Normal velocity
164 // The Physics
165 // -- Density
166 v[0][i] -= wdetJb * rho * u_normal;
167
168 // -- Momentum
169 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho * u_normal * u[j] + norm[j] * P);
170
171 // -- Total Energy Density
172 v[4][i] -= wdetJb * u_normal * (E + P);
173 }
174 return 0;
175 }
176