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