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