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 12 #ifndef channel_h 13 #define channel_h 14 15 #include <ceed.h> 16 #include <math.h> 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, 34 const CeedScalar X[], CeedInt Nf, void *ctx) { 35 36 const ChannelContext context = (ChannelContext)ctx; 37 const CeedScalar theta0 = context->theta0; 38 const CeedScalar P0 = context->P0; 39 const CeedScalar umax = context->umax; 40 const CeedScalar center = context->center; 41 const CeedScalar H = context->H; 42 NewtonianIdealGasContext gas = &context->newtonian_ctx; 43 const CeedScalar cp = gas->cp; 44 const CeedScalar mu = gas->mu; 45 const CeedScalar k = gas->k; 46 // There is a gravity body force but it is excluded from 47 // the potential energy due to periodicity. 48 // g = (g, 0, 0) 49 // x = (0, x_2, x_3) 50 // e_potential = dot(g, x) = 0 51 const CeedScalar x[3] = {0, X[1], X[2]}; 52 53 const CeedScalar Pr = mu / (cp*k); 54 const CeedScalar Ec = (umax*umax) / (cp*theta0); 55 const CeedScalar theta = theta0*(1 + (Pr*Ec/3) 56 * (1 - Square(Square((x[1]-center)/H)))); 57 CeedScalar Y[5] = {0.}; 58 Y[0] = P0; 59 Y[1] = umax*(1 - Square((x[1]-center)/H)); 60 Y[2] = 0.; 61 Y[3] = 0.; 62 Y[4] = theta; 63 64 return StateFromY(gas, Y, x); 65 } 66 67 // ***************************************************************************** 68 // This QFunction set the initial condition 69 // ***************************************************************************** 70 CEED_QFUNCTION(ICsChannel)(void *ctx, CeedInt Q, 71 const CeedScalar *const *in, CeedScalar *const *out) { 72 // Inputs 73 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 74 75 // Outputs 76 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 77 78 // Context 79 const ChannelContext context = (ChannelContext)ctx; 80 81 // Quadrature Point Loop 82 CeedPragmaSIMD 83 for (CeedInt i=0; i<Q; i++) { 84 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 85 State s = Exact_Channel(3, 0., x, 5, ctx); 86 CeedScalar q[5] = {0}; 87 switch (context->newtonian_ctx.state_var) { 88 case STATEVAR_CONSERVATIVE: 89 UnpackState_U(s.U, q); 90 break; 91 case STATEVAR_PRIMITIVE: 92 UnpackState_Y(s.Y, q); 93 break; 94 } 95 96 for (CeedInt j=0; j<5; j++) 97 q0[j][i] = q[j]; 98 99 } // End of Quadrature Point Loop 100 return 0; 101 } 102 103 // ***************************************************************************** 104 // This QFunction set the inflow boundary condition for conservative variables 105 // ***************************************************************************** 106 CEED_QFUNCTION(Channel_Inflow)(void *ctx, CeedInt Q, 107 const CeedScalar *const *in, 108 CeedScalar *const *out) { 109 // *INDENT-OFF* 110 // Inputs 111 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 112 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 113 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 114 115 // Outputs 116 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 117 // *INDENT-ON* 118 const ChannelContext context = (ChannelContext)ctx; 119 const bool implicit = context->implicit; 120 NewtonianIdealGasContext gas = &context->newtonian_ctx; 121 const CeedScalar cv = gas->cv; 122 const CeedScalar cp = gas->cp; 123 const CeedScalar gamma = cp / cv; 124 125 CeedPragmaSIMD 126 // Quadrature Point Loop 127 for (CeedInt i=0; i<Q; i++) { 128 // Setup 129 // -- Interp-to-Interp q_data 130 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 131 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 132 // We can effect this by swapping the sign on this weight 133 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 134 135 // There is a gravity body force but it is excluded from 136 // the potential energy due to periodicity. 137 // g = (g, 0, 0) 138 // x = (0, x_2, x_3) 139 // e_potential = dot(g, x) = 0 140 const CeedScalar x[3] = {0, X[1][i], X[2][i]}; 141 142 // Calcualte prescribed inflow values 143 State s_exact = Exact_Channel(3, 0., x, 5, ctx); 144 CeedScalar q_exact[5] = {0.}; 145 UnpackState_U(s_exact.U, q_exact); 146 147 // Find pressure using state inside the domain 148 CeedScalar q_inside[5] = {0}; 149 for (CeedInt j=0; j<5; j++) 150 q_inside[j] = q[j][i]; 151 State s_inside = StateFromU(gas, q_inside, x); 152 const CeedScalar P = s_inside.Y.pressure; 153 154 // Find inflow state using calculated P and prescribed velocity, theta0 155 const CeedScalar e_internal = cv * s_exact.Y.temperature; 156 const CeedScalar rho_in = P / ((gamma - 1) * e_internal); 157 const CeedScalar E_kinetic = .5 * rho_in * Dot3(s_exact.Y.velocity, 158 s_exact.Y.velocity); 159 const CeedScalar E = rho_in * e_internal + E_kinetic; 160 161 // ---- Normal vect 162 const CeedScalar norm[3] = {q_data_sur[1][i], 163 q_data_sur[2][i], 164 q_data_sur[3][i] 165 }; 166 // The Physics 167 // Zero v so all future terms can safely sum into it 168 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 169 170 const CeedScalar u_normal = Dot3(norm, s_exact.Y.velocity); 171 172 // The Physics 173 // -- Density 174 v[0][i] -= wdetJb * rho_in * u_normal; 175 176 // -- Momentum 177 for (CeedInt j=0; j<3; j++) 178 v[j+1][i] -= wdetJb * (rho_in * u_normal * s_exact.Y.velocity[j] + 179 norm[j] * P); 180 181 // -- Total Energy Density 182 v[4][i] -= wdetJb * u_normal * (E + P); 183 184 } // End Quadrature Point Loop 185 return 0; 186 } 187 188 // ***************************************************************************** 189 // This QFunction set the outflow boundary condition for conservative variables 190 // ***************************************************************************** 191 CEED_QFUNCTION(Channel_Outflow)(void *ctx, CeedInt Q, 192 const CeedScalar *const *in, 193 CeedScalar *const *out) { 194 // *INDENT-OFF* 195 // Inputs 196 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 197 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 198 199 // Outputs 200 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 201 // *INDENT-ON* 202 203 const ChannelContext context = (ChannelContext)ctx; 204 const bool implicit = context->implicit; 205 const CeedScalar P0 = context->P0; 206 207 CeedPragmaSIMD 208 // Quadrature Point Loop 209 for (CeedInt i=0; i<Q; i++) { 210 // Setup 211 // -- Interp in 212 const CeedScalar rho = q[0][i]; 213 const CeedScalar u[3] = {q[1][i] / rho, 214 q[2][i] / rho, 215 q[3][i] / rho 216 }; 217 const CeedScalar E = q[4][i]; 218 219 // -- Interp-to-Interp q_data 220 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 221 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 222 // We can effect this by swapping the sign on this weight 223 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 224 225 // ---- Normal vect 226 const CeedScalar norm[3] = {q_data_sur[1][i], 227 q_data_sur[2][i], 228 q_data_sur[3][i] 229 }; 230 // The Physics 231 // Zero v so all future terms can safely sum into it 232 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 233 234 // Implementing outflow condition 235 const CeedScalar P = P0; // pressure 236 const CeedScalar u_normal = Dot3(norm, u); // Normal velocity 237 // The Physics 238 // -- Density 239 v[0][i] -= wdetJb * rho * u_normal; 240 241 // -- Momentum 242 for (CeedInt j=0; j<3; j++) 243 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 244 245 // -- Total Energy Density 246 v[4][i] -= wdetJb * u_normal * (E + P); 247 248 } // End Quadrature Point Loop 249 return 0; 250 } 251 252 // ***************************************************************************** 253 #endif // channel_h 254