xref: /honee/qfunctions/channel.h (revision ea615d4cc464aa6ad650c06fae6d120cc2465bc4) 
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_ newtonian_ctx;
22 };
23 
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   NewtonianIdealGasContext gas     = &context->newtonian_ctx;
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 // *****************************************************************************
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   const NewtonianIdealGasContext gas     = &context->newtonian_ctx;
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, gas->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 // *****************************************************************************
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   NewtonianIdealGasContext gas         = &context->newtonian_ctx;
87   const CeedScalar         gamma       = HeatCapacityRatio(&context->newtonian_ctx);
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 // *****************************************************************************
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