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 <math.h> 16 #include <ceed/ceed.h> 17 #include "newtonian_types.h" 18 #include "newtonian_state.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 if (context->newtonian_ctx.is_primitive) 88 UnpackState_Y(s.Y, q); 89 else 90 UnpackState_U(s.U, q); 91 92 for (CeedInt j=0; j<5; j++) 93 q0[j][i] = q[j]; 94 95 } // End of Quadrature Point Loop 96 return 0; 97 } 98 99 // ***************************************************************************** 100 // This QFunction set the inflow boundary condition for conservative variables 101 // ***************************************************************************** 102 CEED_QFUNCTION(Channel_Inflow)(void *ctx, CeedInt Q, 103 const CeedScalar *const *in, 104 CeedScalar *const *out) { 105 // *INDENT-OFF* 106 // Inputs 107 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 108 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 109 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 110 111 // Outputs 112 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 113 // *INDENT-ON* 114 const ChannelContext context = (ChannelContext)ctx; 115 const bool implicit = context->implicit; 116 NewtonianIdealGasContext gas = &context->newtonian_ctx; 117 const CeedScalar cv = gas->cv; 118 const CeedScalar cp = gas->cp; 119 const CeedScalar gamma = cp / cv; 120 121 CeedPragmaSIMD 122 // Quadrature Point Loop 123 for (CeedInt i=0; i<Q; i++) { 124 // Setup 125 // -- Interp-to-Interp q_data 126 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 127 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 128 // We can effect this by swapping the sign on this weight 129 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 130 131 // There is a gravity body force but it is excluded from 132 // the potential energy due to periodicity. 133 // g = (g, 0, 0) 134 // x = (0, x_2, x_3) 135 // e_potential = dot(g, x) = 0 136 const CeedScalar x[3] = {0, X[1][i], X[2][i]}; 137 138 // Calcualte prescribed inflow values 139 State s_exact = Exact_Channel(3, 0., x, 5, ctx); 140 CeedScalar q_exact[5] = {0.}; 141 UnpackState_U(s_exact.U, q_exact); 142 143 // Find pressure using state inside the domain 144 CeedScalar q_inside[5] = {0}; 145 for (CeedInt j; j<5; j++) 146 q_inside[j] = q[j][i]; 147 State s_inside = StateFromU(gas, q_inside, x); 148 const CeedScalar P = s_inside.Y.pressure; 149 150 // Find inflow state using calculated P and prescribed velocity, theta0 151 const CeedScalar e_internal = cv * s_exact.Y.temperature; 152 const CeedScalar rho_in = P / ((gamma - 1) * e_internal); 153 const CeedScalar E_kinetic = .5 * rho_in * Dot3(s_exact.Y.velocity, 154 s_exact.Y.velocity); 155 const CeedScalar E = rho_in * e_internal + E_kinetic; 156 157 // ---- Normal vect 158 const CeedScalar norm[3] = {q_data_sur[1][i], 159 q_data_sur[2][i], 160 q_data_sur[3][i] 161 }; 162 // The Physics 163 // Zero v so all future terms can safely sum into it 164 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 165 166 const CeedScalar u_normal = Dot3(norm, s_exact.Y.velocity); 167 168 // The Physics 169 // -- Density 170 v[0][i] -= wdetJb * rho_in * u_normal; 171 172 // -- Momentum 173 for (CeedInt j=0; j<3; j++) 174 v[j+1][i] -= wdetJb * (rho_in * u_normal * s_exact.Y.velocity[j] + 175 norm[j] * P); 176 177 // -- Total Energy Density 178 v[4][i] -= wdetJb * u_normal * (E + P); 179 180 } // End Quadrature Point Loop 181 return 0; 182 } 183 184 // ***************************************************************************** 185 // This QFunction set the outflow boundary condition for conservative variables 186 // ***************************************************************************** 187 CEED_QFUNCTION(Channel_Outflow)(void *ctx, CeedInt Q, 188 const CeedScalar *const *in, 189 CeedScalar *const *out) { 190 // *INDENT-OFF* 191 // Inputs 192 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 193 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 194 195 // Outputs 196 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 197 // *INDENT-ON* 198 199 const ChannelContext context = (ChannelContext)ctx; 200 const bool implicit = context->implicit; 201 const CeedScalar P0 = context->P0; 202 203 CeedPragmaSIMD 204 // Quadrature Point Loop 205 for (CeedInt i=0; i<Q; i++) { 206 // Setup 207 // -- Interp in 208 const CeedScalar rho = q[0][i]; 209 const CeedScalar u[3] = {q[1][i] / rho, 210 q[2][i] / rho, 211 q[3][i] / rho 212 }; 213 const CeedScalar E = q[4][i]; 214 215 // -- Interp-to-Interp q_data 216 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 217 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 218 // We can effect this by swapping the sign on this weight 219 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 220 221 // ---- Normal vect 222 const CeedScalar norm[3] = {q_data_sur[1][i], 223 q_data_sur[2][i], 224 q_data_sur[3][i] 225 }; 226 // The Physics 227 // Zero v so all future terms can safely sum into it 228 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 229 230 // Implementing outflow condition 231 const CeedScalar P = P0; // pressure 232 const CeedScalar u_normal = Dot3(norm, u); // Normal velocity 233 // The Physics 234 // -- Density 235 v[0][i] -= wdetJb * rho * u_normal; 236 237 // -- Momentum 238 for (CeedInt j=0; j<3; j++) 239 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 240 241 // -- Total Energy Density 242 v[4][i] -= wdetJb * u_normal * (E + P); 243 244 } // End Quadrature Point Loop 245 return 0; 246 } 247 248 // ***************************************************************************** 249 #endif // channel_h 250