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