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