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