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 gas->g[0] = 0.; 49 gas->g[1] = 0.; 50 gas->g[2] = 0.; 51 52 const CeedScalar y = X[1]; 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((y-center)/H)))); 57 CeedScalar Y[5] = {0.}; 58 Y[0] = P0; 59 Y[1] = umax*(1 - Square((y-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 if (context->newtonian_ctx.is_primitive) { 87 q0[0][i] = s.Y.pressure; 88 for (CeedInt j=0; j<3; j++) 89 q0[j+1][i] = s.Y.velocity[j]; 90 q0[4][i] = s.Y.temperature; 91 } else { 92 q0[0][i] = s.U.density; 93 for (CeedInt j=0; j<3; j++) 94 q0[j+1][i] = s.U.momentum[j]; 95 q0[4][i] = s.U.E_total; 96 } 97 98 } // End of Quadrature Point Loop 99 return 0; 100 } 101 102 // ***************************************************************************** 103 CEED_QFUNCTION(Channel_Inflow)(void *ctx, CeedInt Q, 104 const CeedScalar *const *in, 105 CeedScalar *const *out) { 106 // *INDENT-OFF* 107 // Inputs 108 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 109 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 110 (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 111 112 // Outputs 113 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 114 // *INDENT-ON* 115 const ChannelContext context = (ChannelContext)ctx; 116 const bool implicit = context->implicit; 117 const CeedScalar cv = context->newtonian_ctx.cv; 118 const CeedScalar cp = context->newtonian_ctx.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 // Calcualte prescribed inflow values 132 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 133 State s = Exact_Channel(3, 0., x, 5, ctx); 134 CeedScalar q_exact[5] = {0.}; 135 q_exact[0] = s.U.density; 136 for (CeedInt j=0; j<3; j++) 137 q_exact[j+1] = s.U.momentum[j]; 138 q_exact[4] = s.U.E_total; 139 const CeedScalar E_kinetic_exact = 0.5*Dot3(&q_exact[1], &q_exact[1]) 140 / q_exact[0]; 141 const CeedScalar velocity[3] = {q_exact[1]/q_exact[0], 142 q_exact[2]/q_exact[0], 143 q_exact[3]/q_exact[0] 144 }; 145 const CeedScalar theta = (q_exact[4] - E_kinetic_exact) / (q_exact[0]*cv); 146 147 // Find pressure using state inside the domain 148 const CeedScalar rho = q[0][i]; 149 const CeedScalar u[3] = {q[1][i]/rho, q[2][i]/rho, q[3][i]/rho}; 150 const CeedScalar E_internal = q[4][i] - .5 * rho * Dot3(u,u); 151 const CeedScalar P = E_internal * (gamma - 1.); 152 153 // Find inflow state using calculated P and prescribed velocity, theta0 154 const CeedScalar e_internal = cv * theta; 155 const CeedScalar rho_in = P / ((gamma - 1) * e_internal); 156 const CeedScalar E_kinetic = .5 * rho_in * Dot3(velocity, velocity); 157 const CeedScalar E = rho_in * e_internal + E_kinetic; 158 // ---- Normal vect 159 const CeedScalar norm[3] = {q_data_sur[1][i], 160 q_data_sur[2][i], 161 q_data_sur[3][i] 162 }; 163 164 // The Physics 165 // Zero v so all future terms can safely sum into it 166 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 167 168 const CeedScalar u_normal = Dot3(norm, velocity); 169 170 // The Physics 171 // -- Density 172 v[0][i] -= wdetJb * rho_in * u_normal; 173 174 // -- Momentum 175 for (CeedInt j=0; j<3; j++) 176 v[j+1][i] -= wdetJb * (rho_in * u_normal * velocity[j] + 177 norm[j] * P); 178 179 // -- Total Energy Density 180 v[4][i] -= wdetJb * u_normal * (E + P); 181 182 } // End Quadrature Point Loop 183 return 0; 184 } 185 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 227 // The Physics 228 // Zero v so all future terms can safely sum into it 229 for (CeedInt j=0; j<5; j++) v[j][i] = 0.; 230 231 // Implementing outflow condition 232 const CeedScalar P = P0; // pressure 233 const CeedScalar u_normal = Dot3(norm, u); // Normal velocity 234 // The Physics 235 // -- Density 236 v[0][i] -= wdetJb * rho * u_normal; 237 238 // -- Momentum 239 for (CeedInt j=0; j<3; j++) 240 v[j+1][i] -= wdetJb *(rho * u_normal * u[j] + norm[j] * P); 241 242 // -- Total Energy Density 243 v[4][i] -= wdetJb * u_normal * (E + P); 244 245 } // End Quadrature Point Loop 246 return 0; 247 } 248 249 // ***************************************************************************** 250 #endif // channel_h 251