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