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 19 typedef struct ChannelContext_ *ChannelContext; 20 struct ChannelContext_ { 21 bool implicit; // !< Using implicit timesteping or not 22 CeedScalar theta0; // !< Reference temperature 23 CeedScalar P0; // !< Reference Pressure 24 CeedScalar umax; // !< Centerline velocity 25 CeedScalar center; // !< Y Coordinate for center of channel 26 CeedScalar H; // !< Channel half-height 27 CeedScalar B; // !< Body-force driving the flow 28 struct NewtonianIdealGasContext_ newtonian_ctx; 29 }; 30 31 CEED_QFUNCTION_HELPER CeedInt Exact_Channel(CeedInt dim, CeedScalar time, 32 const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) { 33 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 const CeedScalar cv = context->newtonian_ctx.cv; 41 const CeedScalar cp = context->newtonian_ctx.cp; 42 const CeedScalar Rd = cp - cv; 43 const CeedScalar mu = context->newtonian_ctx.mu; 44 const CeedScalar k = context->newtonian_ctx.k; 45 46 const CeedScalar y=X[1]; 47 48 const CeedScalar Pr = mu / (cp*k); 49 const CeedScalar Ec = (umax*umax) / (cp*theta0); 50 const CeedScalar theta = theta0*(1 + (Pr*Ec/3) 51 * (1 - Square(Square((y-center)/H)))); 52 53 const CeedScalar p = P0; 54 55 const CeedScalar rho = p / (Rd*theta); 56 57 q[0] = rho; 58 q[1] = rho * umax*(1 - Square((y-center)/H)); 59 q[2] = 0; 60 q[3] = 0; 61 q[4] = rho * (cv*theta) + .5 * (q[1]*q[1] + q[2]*q[2] + q[3]*q[3]) / rho; 62 63 return 0; 64 } 65 66 // ***************************************************************************** 67 // This QFunction sets the initial condition 68 // ***************************************************************************** 69 CEED_QFUNCTION(ICsChannel)(void *ctx, CeedInt Q, 70 const CeedScalar *const *in, CeedScalar *const *out) { 71 // Inputs 72 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 73 74 // Outputs 75 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 76 77 // Quadrature Point Loop 78 CeedPragmaSIMD 79 for (CeedInt i=0; i<Q; i++) { 80 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 81 CeedScalar q[5] = {0.}; 82 Exact_Channel(3, 0., x, 5, q, ctx); 83 84 for (CeedInt j=0; j<5; j++) 85 q0[j][i] = q[j]; 86 } // End of Quadrature Point Loop 87 return 0; 88 } 89 90 // ***************************************************************************** 91 CEED_QFUNCTION(Channel_Inflow)(void *ctx, CeedInt Q, 92 const CeedScalar *const *in, 93 CeedScalar *const *out) { 94 // *INDENT-OFF* 95 // Inputs 96 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 97 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1], 98 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 99 100 // Outputs 101 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 102 // *INDENT-ON* 103 const ChannelContext context = (ChannelContext)ctx; 104 const bool implicit = context->implicit; 105 const CeedScalar cv = context->newtonian_ctx.cv; 106 const CeedScalar cp = context->newtonian_ctx.cp; 107 const CeedScalar gamma = cp/cv; 108 109 CeedPragmaSIMD 110 // Quadrature Point Loop 111 for (CeedInt i=0; i<Q; i++) { 112 // Setup 113 // -- Interp-to-Interp q_data 114 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 115 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 116 // We can effect this by swapping the sign on this weight 117 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 118 119 // Calcualte prescribed inflow values 120 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 121 CeedScalar q_exact[5] = {0.}; 122 Exact_Channel(3, 0., x, 5, q_exact, ctx); 123 const CeedScalar E_kinetic_exact = 0.5*(q_exact[1]*q_exact[1] + 124 q_exact[2]*q_exact[2] + 125 q_exact[3]*q_exact[3]) / q_exact[0]; 126 const CeedScalar velocity[3] = {q_exact[1]/q_exact[0], 127 q_exact[2]/q_exact[0], 128 q_exact[3]/q_exact[0] 129 }; 130 const CeedScalar theta = (q_exact[4] - E_kinetic_exact) / (q_exact[0]*cv); 131 132 // Find pressure using state inside the domain 133 const CeedScalar rho = q[0][i]; 134 const CeedScalar u[3] = {q[1][i]/rho, q[2][i]/rho, q[3][i]/rho}; 135 const CeedScalar E_internal = q[4][i] - .5 * rho * (u[0]*u[0] + u[1]*u[1] + 136 u[2]*u[2]); 137 const CeedScalar P = E_internal * (gamma - 1.); 138 139 // Find inflow state using calculated P and prescribed velocity, theta0 140 const CeedScalar e_internal = cv * theta; 141 const CeedScalar rho_in = P / ((gamma - 1) * e_internal); 142 const CeedScalar E_kinetic = .5 * rho_in * (velocity[0]*velocity[0] + 143 velocity[1]*velocity[1] + 144 velocity[2]*velocity[2]); 145 const CeedScalar E = rho_in * e_internal + E_kinetic; 146 // ---- Normal vect 147 const CeedScalar norm[3] = {q_data_sur[1][i], 148 q_data_sur[2][i], 149 q_data_sur[3][i] 150 }; 151 152 // The Physics 153 // Zero v so all future terms can safely sum into it 154 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 155 156 const CeedScalar u_normal = norm[0]*velocity[0] + 157 norm[1]*velocity[1] + 158 norm[2]*velocity[2]; 159 160 // The Physics 161 // -- Density 162 v[0][i] -= wdetJb * rho_in * u_normal; 163 164 // -- Momentum 165 for (CeedInt j=0; j<3; j++) 166 v[j+1][i] -= wdetJb * (rho_in * u_normal * velocity[j] + 167 norm[j] * P); 168 169 // -- Total Energy Density 170 v[4][i] -= wdetJb * u_normal * (E + P); 171 172 } // End Quadrature Point Loop 173 return 0; 174 } 175 176 // ***************************************************************************** 177 CEED_QFUNCTION(Channel_Outflow)(void *ctx, CeedInt Q, 178 const CeedScalar *const *in, 179 CeedScalar *const *out) { 180 // *INDENT-OFF* 181 // Inputs 182 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 183 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 184 // Outputs 185 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 186 // *INDENT-ON* 187 188 const ChannelContext context = (ChannelContext)ctx; 189 const bool implicit = context->implicit; 190 const CeedScalar P0 = context->P0; 191 192 CeedPragmaSIMD 193 // Quadrature Point Loop 194 for (CeedInt i=0; i<Q; i++) { 195 // Setup 196 // -- Interp in 197 const CeedScalar rho = q[0][i]; 198 const CeedScalar u[3] = {q[1][i] / rho, 199 q[2][i] / rho, 200 q[3][i] / rho 201 }; 202 const CeedScalar E = q[4][i]; 203 204 // -- Interp-to-Interp q_data 205 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 206 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 207 // We can effect this by swapping the sign on this weight 208 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 209 210 // ---- Normal vect 211 const CeedScalar norm[3] = {q_data_sur[1][i], 212 q_data_sur[2][i], 213 q_data_sur[3][i] 214 }; 215 216 // The Physics 217 // Zero v so all future terms can safely sum into it 218 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 219 220 // Implementing outflow condition 221 const CeedScalar P = P0; // pressure 222 const CeedScalar u_normal = norm[0]*u[0] + norm[1]*u[1] + 223 norm[2]*u[2]; // Normal velocity 224 // The Physics 225 // -- Density 226 v[0][i] -= wdetJb * rho * u_normal; 227 228 // -- Momentum 229 for (CeedInt j=0; j<3; j++) 230 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 231 232 // -- Total Energy Density 233 v[4][i] -= wdetJb * u_normal * (E + P); 234 235 } // End Quadrature Point Loop 236 return 0; 237 } 238 #endif // channel_h 239