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 #ifndef channel_h 12 #define channel_h 13 14 #include <ceed.h> 15 #include <math.h> 16 17 #include "newtonian_state.h" 18 #include "newtonian_types.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, const CeedScalar X[], CeedInt Nf, void *ctx) { 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 NewtonianIdealGasContext gas = &context->newtonian_ctx; 41 const CeedScalar cp = gas->cp; 42 const CeedScalar mu = gas->mu; 43 const CeedScalar k = gas->k; 44 // There is a gravity body force but it is excluded from 45 // the potential energy due to periodicity. 46 // g = (g, 0, 0) 47 // x = (0, x_2, x_3) 48 // e_potential = dot(g, x) = 0 49 const CeedScalar x[3] = {0, X[1], X[2]}; 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 - Square(Square((x[1] - center) / H)))); 54 CeedScalar Y[5] = {0.}; 55 Y[0] = P0; 56 Y[1] = umax * (1 - Square((x[1] - center) / H)); 57 Y[2] = 0.; 58 Y[3] = 0.; 59 Y[4] = theta; 60 61 return StateFromY(gas, Y, x); 62 } 63 64 // ***************************************************************************** 65 // This QFunction set the initial condition 66 // ***************************************************************************** 67 CEED_QFUNCTION(ICsChannel)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 68 // Inputs 69 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 70 71 // Outputs 72 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 73 74 // Context 75 const ChannelContext context = (ChannelContext)ctx; 76 77 // Quadrature Point Loop 78 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 79 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 80 State s = Exact_Channel(3, 0., x, 5, ctx); 81 CeedScalar q[5] = {0}; 82 switch (context->newtonian_ctx.state_var) { 83 case STATEVAR_CONSERVATIVE: 84 UnpackState_U(s.U, q); 85 break; 86 case STATEVAR_PRIMITIVE: 87 UnpackState_Y(s.Y, q); 88 break; 89 } 90 91 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 92 93 } // End of Quadrature Point Loop 94 return 0; 95 } 96 97 // ***************************************************************************** 98 // This QFunction set the inflow boundary condition for conservative variables 99 // ***************************************************************************** 100 CEED_QFUNCTION(Channel_Inflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 101 // Inputs 102 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 103 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 104 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 105 106 // Outputs 107 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 108 109 const ChannelContext context = (ChannelContext)ctx; 110 const bool implicit = context->implicit; 111 NewtonianIdealGasContext gas = &context->newtonian_ctx; 112 const CeedScalar cv = gas->cv; 113 const CeedScalar gamma = HeatCapacityRatio(&context->newtonian_ctx); 114 115 // Quadrature Point Loop 116 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 117 // Setup 118 // -- Interp-to-Interp q_data 119 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 120 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 121 // We can effect this by swapping the sign on this weight 122 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 123 124 // There is a gravity body force but it is excluded from 125 // the potential energy due to periodicity. 126 // g = (g, 0, 0) 127 // x = (0, x_2, x_3) 128 // e_potential = dot(g, x) = 0 129 const CeedScalar x[3] = {0, X[1][i], X[2][i]}; 130 131 // Calcualte prescribed inflow values 132 State s_exact = Exact_Channel(3, 0., x, 5, ctx); 133 CeedScalar q_exact[5] = {0.}; 134 UnpackState_U(s_exact.U, q_exact); 135 136 // Find pressure using state inside the domain 137 CeedScalar q_inside[5] = {0}; 138 for (CeedInt j = 0; j < 5; j++) q_inside[j] = q[j][i]; 139 State s_inside = StateFromU(gas, q_inside, x); 140 const CeedScalar P = s_inside.Y.pressure; 141 142 // Find inflow state using calculated P and prescribed velocity, theta0 143 const CeedScalar e_internal = cv * s_exact.Y.temperature; 144 const CeedScalar rho_in = P / ((gamma - 1) * e_internal); 145 const CeedScalar E_kinetic = .5 * rho_in * Dot3(s_exact.Y.velocity, s_exact.Y.velocity); 146 const CeedScalar E = rho_in * e_internal + E_kinetic; 147 148 // ---- Normal vect 149 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 150 // The Physics 151 // Zero v so all future terms can safely sum into it 152 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.; 153 154 const CeedScalar u_normal = Dot3(norm, s_exact.Y.velocity); 155 156 // The Physics 157 // -- Density 158 v[0][i] -= wdetJb * rho_in * u_normal; 159 160 // -- Momentum 161 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho_in * u_normal * s_exact.Y.velocity[j] + norm[j] * P); 162 163 // -- Total Energy Density 164 v[4][i] -= wdetJb * u_normal * (E + P); 165 166 } // End Quadrature Point Loop 167 return 0; 168 } 169 170 // ***************************************************************************** 171 // This QFunction set the outflow boundary condition for conservative variables 172 // ***************************************************************************** 173 CEED_QFUNCTION(Channel_Outflow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 174 // Inputs 175 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 176 const CeedScalar(*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 177 178 // Outputs 179 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 180 181 const ChannelContext context = (ChannelContext)ctx; 182 const bool implicit = context->implicit; 183 const CeedScalar P0 = context->P0; 184 185 // Quadrature Point Loop 186 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 187 // Setup 188 // -- Interp in 189 const CeedScalar rho = q[0][i]; 190 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 191 const CeedScalar E = q[4][i]; 192 193 // -- Interp-to-Interp q_data 194 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 195 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 196 // We can effect this by swapping the sign on this weight 197 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 198 199 // ---- Normal vect 200 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 201 // The Physics 202 // Zero v so all future terms can safely sum into it 203 for (CeedInt j = 0; j < 5; j++) v[j][i] = 0.; 204 205 // Implementing outflow condition 206 const CeedScalar P = P0; // pressure 207 const CeedScalar u_normal = Dot3(norm, u); // Normal velocity 208 // The Physics 209 // -- Density 210 v[0][i] -= wdetJb * rho * u_normal; 211 212 // -- Momentum 213 for (CeedInt j = 0; j < 3; j++) v[j + 1][i] -= wdetJb * (rho * u_normal * u[j] + norm[j] * P); 214 215 // -- Total Energy Density 216 v[4][i] -= wdetJb * u_normal * (E + P); 217 218 } // End Quadrature Point Loop 219 return 0; 220 } 221 222 // ***************************************************************************** 223 #endif // channel_h 224