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 newtonian_h 13 #define newtonian_h 14 15 #include <ceed.h> 16 #include <math.h> 17 #include "newtonian_state.h" 18 #include "newtonian_types.h" 19 #include "stabilization.h" 20 #include "utils.h" 21 22 // ***************************************************************************** 23 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 24 // ***************************************************************************** 25 CEED_QFUNCTION(ICsNewtonianIG)(void *ctx, CeedInt Q, 26 const CeedScalar *const *in, CeedScalar *const *out) { 27 // Inputs 28 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 29 30 // Outputs 31 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 32 33 // Context 34 const SetupContext context = (SetupContext)ctx; 35 const CeedScalar theta0 = context->theta0; 36 const CeedScalar P0 = context->P0; 37 const CeedScalar cv = context->cv; 38 const CeedScalar cp = context->cp; 39 const CeedScalar *g = context->g; 40 const CeedScalar Rd = cp - cv; 41 42 // Quadrature Point Loop 43 CeedPragmaSIMD 44 for (CeedInt i=0; i<Q; i++) { 45 CeedScalar q[5] = {0.}; 46 47 // Setup 48 // -- Coordinates 49 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 50 const CeedScalar e_potential = -Dot3(g, x); 51 52 // -- Density 53 const CeedScalar rho = P0 / (Rd*theta0); 54 55 // Initial Conditions 56 q[0] = rho; 57 q[1] = 0.0; 58 q[2] = 0.0; 59 q[3] = 0.0; 60 q[4] = rho * (cv*theta0 + e_potential); 61 62 for (CeedInt j=0; j<5; j++) 63 q0[j][i] = q[j]; 64 65 } // End of Quadrature Point Loop 66 return 0; 67 } 68 69 // ***************************************************************************** 70 // This QFunction sets a "still" initial condition for generic Newtonian IG 71 // problems in primitive variables 72 // ***************************************************************************** 73 CEED_QFUNCTION(ICsNewtonianIG_Prim)(void *ctx, CeedInt Q, 74 const CeedScalar *const *in, CeedScalar *const *out) { 75 // Outputs 76 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 77 78 // Context 79 const SetupContext context = (SetupContext)ctx; 80 const CeedScalar theta0 = context->theta0; 81 const CeedScalar P0 = context->P0; 82 83 // Quadrature Point Loop 84 CeedPragmaSIMD 85 for (CeedInt i=0; i<Q; i++) { 86 CeedScalar q[5] = {0.}; 87 88 // Initial Conditions 89 q[0] = P0; 90 q[1] = 0.0; 91 q[2] = 0.0; 92 q[3] = 0.0; 93 q[4] = theta0; 94 95 for (CeedInt j=0; j<5; j++) 96 q0[j][i] = q[j]; 97 98 } // End of Quadrature Point Loop 99 return 0; 100 } 101 102 // ***************************************************************************** 103 // This QFunction implements the following formulation of Navier-Stokes with 104 // explicit time stepping method 105 // 106 // This is 3D compressible Navier-Stokes in conservation form with state 107 // variables of density, momentum density, and total energy density. 108 // 109 // State Variables: q = ( rho, U1, U2, U3, E ) 110 // rho - Mass Density 111 // Ui - Momentum Density, Ui = rho ui 112 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 113 // 114 // Navier-Stokes Equations: 115 // drho/dt + div( U ) = 0 116 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 117 // dE/dt + div( (E + P) u ) = div( Fe ) 118 // 119 // Viscous Stress: 120 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 121 // 122 // Thermal Stress: 123 // Fe = u Fu + k grad( T ) 124 // Equation of State 125 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 126 // 127 // Stabilization: 128 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 129 // f1 = rho sqrt(ui uj gij) 130 // gij = dXi/dX * dXi/dX 131 // TauC = Cc f1 / (8 gii) 132 // TauM = min( 1 , 1 / f1 ) 133 // TauE = TauM / (Ce cv) 134 // 135 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 136 // 137 // Constants: 138 // lambda = - 2 / 3, From Stokes hypothesis 139 // mu , Dynamic viscosity 140 // k , Thermal conductivity 141 // cv , Specific heat, constant volume 142 // cp , Specific heat, constant pressure 143 // g , Gravity 144 // gamma = cp / cv, Specific heat ratio 145 // 146 // We require the product of the inverse of the Jacobian (dXdx_j,k) and 147 // its transpose (dXdx_k,j) to properly compute integrals of the form: 148 // int( gradv gradu ) 149 // 150 // ***************************************************************************** 151 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, 152 const CeedScalar *const *in, CeedScalar *const *out) { 153 // *INDENT-OFF* 154 // Inputs 155 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 156 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 157 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 158 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 159 // Outputs 160 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 161 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 162 // *INDENT-ON* 163 164 // Context 165 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 166 const CeedScalar *g = context->g; 167 const CeedScalar dt = context->dt; 168 169 CeedPragmaSIMD 170 // Quadrature Point Loop 171 for (CeedInt i=0; i<Q; i++) { 172 CeedScalar U[5]; 173 for (int j=0; j<5; j++) U[j] = q[j][i]; 174 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 175 State s = StateFromU(context, U, x_i); 176 177 // -- Interp-to-Interp q_data 178 const CeedScalar wdetJ = q_data[0][i]; 179 // -- Interp-to-Grad q_data 180 // ---- Inverse of change of coordinate matrix: X_i,j 181 // *INDENT-OFF* 182 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 183 {q_data[4][i], q_data[5][i], q_data[6][i]}, 184 {q_data[7][i], q_data[8][i], q_data[9][i]} 185 }; 186 // *INDENT-ON* 187 State grad_s[3]; 188 for (CeedInt j=0; j<3; j++) { 189 CeedScalar dx_i[3] = {0}, dU[5]; 190 for (CeedInt k=0; k<5; k++) 191 dU[k] = Grad_q[0][k][i] * dXdx[0][j] + 192 Grad_q[1][k][i] * dXdx[1][j] + 193 Grad_q[2][k][i] * dXdx[2][j]; 194 dx_i[j] = 1.; 195 grad_s[j] = StateFromU_fwd(context, s, dU, x_i, dx_i); 196 } 197 198 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 199 KMStrainRate(grad_s, strain_rate); 200 NewtonianStress(context, strain_rate, kmstress); 201 KMUnpack(kmstress, stress); 202 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 203 204 StateConservative F_inviscid[3]; 205 FluxInviscid(context, s, F_inviscid); 206 207 // Total flux 208 CeedScalar Flux[5][3]; 209 FluxTotal(F_inviscid, stress, Fe, Flux); 210 211 for (CeedInt j=0; j<3; j++) 212 for (CeedInt k=0; k<5; k++) 213 Grad_v[j][k][i] = wdetJ * (dXdx[j][0] * Flux[k][0] + 214 dXdx[j][1] * Flux[k][1] + 215 dXdx[j][2] * Flux[k][2]); 216 217 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 218 for (int j=0; j<5; j++) 219 v[j][i] = wdetJ * body_force[j]; 220 221 // -- Stabilization method: none (Galerkin), SU, or SUPG 222 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}; 223 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 224 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 225 226 for (CeedInt j=0; j<5; j++) 227 for (CeedInt k=0; k<3; k++) 228 Grad_v[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 229 stab[j][1] * dXdx[k][1] + 230 stab[j][2] * dXdx[k][2]); 231 232 } // End Quadrature Point Loop 233 234 // Return 235 return 0; 236 } 237 238 // ***************************************************************************** 239 // This QFunction implements the Navier-Stokes equations (mentioned above) with 240 // implicit time stepping method 241 // 242 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 243 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 244 // (diffussive terms will be added later) 245 // 246 // ***************************************************************************** 247 CEED_QFUNCTION_HELPER int IFunction_Newtonian(void *ctx, CeedInt Q, 248 const CeedScalar *const *in, CeedScalar *const *out, 249 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 250 // *INDENT-OFF* 251 // Inputs 252 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 253 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 254 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 255 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 256 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 257 // Outputs 258 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 259 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1], 260 (*jac_data)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[2]; 261 // *INDENT-ON* 262 // Context 263 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 264 const CeedScalar *g = context->g; 265 const CeedScalar dt = context->dt; 266 267 CeedPragmaSIMD 268 // Quadrature Point Loop 269 for (CeedInt i=0; i<Q; i++) { 270 CeedScalar qi[5]; 271 for (CeedInt j=0; j<5; j++) qi[j] = q[j][i]; 272 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 273 State s = StateFromQi(context, qi, x_i); 274 275 // -- Interp-to-Interp q_data 276 const CeedScalar wdetJ = q_data[0][i]; 277 // -- Interp-to-Grad q_data 278 // ---- Inverse of change of coordinate matrix: X_i,j 279 // *INDENT-OFF* 280 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 281 {q_data[4][i], q_data[5][i], q_data[6][i]}, 282 {q_data[7][i], q_data[8][i], q_data[9][i]} 283 }; 284 // *INDENT-ON* 285 State grad_s[3]; 286 for (CeedInt j=0; j<3; j++) { 287 CeedScalar dx_i[3] = {0}, dqi[5]; 288 for (CeedInt k=0; k<5; k++) 289 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 290 Grad_q[1][k][i] * dXdx[1][j] + 291 Grad_q[2][k][i] * dXdx[2][j]; 292 dx_i[j] = 1.; 293 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 294 } 295 296 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 297 KMStrainRate(grad_s, strain_rate); 298 NewtonianStress(context, strain_rate, kmstress); 299 KMUnpack(kmstress, stress); 300 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 301 302 StateConservative F_inviscid[3]; 303 FluxInviscid(context, s, F_inviscid); 304 305 // Total flux 306 CeedScalar Flux[5][3]; 307 FluxTotal(F_inviscid, stress, Fe, Flux); 308 309 for (CeedInt j=0; j<3; j++) 310 for (CeedInt k=0; k<5; k++) 311 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 312 dXdx[j][1] * Flux[k][1] + 313 dXdx[j][2] * Flux[k][2]); 314 315 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 316 317 // -- Stabilization method: none (Galerkin), SU, or SUPG 318 CeedScalar Tau_d[3], stab[5][3], U_dot[5] = {0}, qi_dot[5], dx0[3] = {0}; 319 for (int j=0; j<5; j++) qi_dot[j] = q_dot[j][i]; 320 State s_dot = StateFromQi_fwd(context, s, qi_dot, x_i, dx0); 321 UnpackState_U(s_dot.U, U_dot); 322 323 for (CeedInt j=0; j<5; j++) 324 v[j][i] = wdetJ * (U_dot[j] - body_force[j]); 325 Tau_diagPrim(context, s, dXdx, dt, Tau_d); 326 Stabilization(context, s, Tau_d, grad_s, U_dot, body_force, x_i, stab); 327 328 for (CeedInt j=0; j<5; j++) 329 for (CeedInt k=0; k<3; k++) 330 Grad_v[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 331 stab[j][1] * dXdx[k][1] + 332 stab[j][2] * dXdx[k][2]); 333 334 for (CeedInt j=0; j<5; j++) jac_data[j][i] = qi[j]; 335 for (CeedInt j=0; j<6; j++) jac_data[5+j][i] = kmstress[j]; 336 for (CeedInt j=0; j<3; j++) jac_data[5+6+j][i] = Tau_d[j]; 337 338 } // End Quadrature Point Loop 339 340 // Return 341 return 0; 342 } 343 344 CEED_QFUNCTION(IFunction_Newtonian_Conserv)(void *ctx, CeedInt Q, 345 const CeedScalar *const *in, CeedScalar *const *out) { 346 return IFunction_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 347 } 348 349 CEED_QFUNCTION(IFunction_Newtonian_Prim)(void *ctx, CeedInt Q, 350 const CeedScalar *const *in, CeedScalar *const *out) { 351 return IFunction_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 352 } 353 354 // ***************************************************************************** 355 // This QFunction implements the jacobian of the Navier-Stokes equations 356 // for implicit time stepping method. 357 // ***************************************************************************** 358 CEED_QFUNCTION_HELPER int IJacobian_Newtonian(void *ctx, CeedInt Q, 359 const CeedScalar *const *in, CeedScalar *const *out, 360 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 361 // *INDENT-OFF* 362 // Inputs 363 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 364 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 365 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 366 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 367 (*jac_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 368 // Outputs 369 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 370 (*Grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 371 // *INDENT-ON* 372 // Context 373 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 374 const CeedScalar *g = context->g; 375 376 CeedPragmaSIMD 377 // Quadrature Point Loop 378 for (CeedInt i=0; i<Q; i++) { 379 // -- Interp-to-Interp q_data 380 const CeedScalar wdetJ = q_data[0][i]; 381 // -- Interp-to-Grad q_data 382 // ---- Inverse of change of coordinate matrix: X_i,j 383 // *INDENT-OFF* 384 const CeedScalar dXdx[3][3] = {{q_data[1][i], q_data[2][i], q_data[3][i]}, 385 {q_data[4][i], q_data[5][i], q_data[6][i]}, 386 {q_data[7][i], q_data[8][i], q_data[9][i]} 387 }; 388 // *INDENT-ON* 389 390 CeedScalar qi[5], kmstress[6], Tau_d[3] __attribute((unused)); 391 for (int j=0; j<5; j++) qi[j] = jac_data[j][i]; 392 for (int j=0; j<6; j++) kmstress[j] = jac_data[5+j][i]; 393 for (int j=0; j<3; j++) Tau_d[j] = jac_data[5+6+j][i]; 394 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 395 State s = StateFromQi(context, qi, x_i); 396 397 CeedScalar dqi[5], dx0[3] = {0}; 398 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 399 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx0); 400 401 State grad_ds[3]; 402 for (int j=0; j<3; j++) { 403 CeedScalar dqi_j[5]; 404 for (int k=0; k<5; k++) 405 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 406 Grad_dq[1][k][i] * dXdx[1][j] + 407 Grad_dq[2][k][i] * dXdx[2][j]; 408 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx0); 409 } 410 411 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 412 KMStrainRate(grad_ds, dstrain_rate); 413 NewtonianStress(context, dstrain_rate, dkmstress); 414 KMUnpack(dkmstress, dstress); 415 KMUnpack(kmstress, stress); 416 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 417 418 StateConservative dF_inviscid[3]; 419 FluxInviscid_fwd(context, s, ds, dF_inviscid); 420 421 // Total flux 422 CeedScalar dFlux[5][3]; 423 FluxTotal(dF_inviscid, dstress, dFe, dFlux); 424 425 for (int j=0; j<3; j++) 426 for (int k=0; k<5; k++) 427 Grad_v[j][k][i] = -wdetJ * (dXdx[j][0] * dFlux[k][0] + 428 dXdx[j][1] * dFlux[k][1] + 429 dXdx[j][2] * dFlux[k][2]); 430 431 const CeedScalar dbody_force[5] = {0, ds.U.density *g[0], ds.U.density *g[1], ds.U.density *g[2], 0}; 432 CeedScalar dU[5] = {0.}; 433 UnpackState_U(ds.U, dU); 434 for (int j=0; j<5; j++) 435 v[j][i] = wdetJ * (context->ijacobian_time_shift * dU[j] - dbody_force[j]); 436 437 // -- Stabilization method: none (Galerkin), SU, or SUPG 438 CeedScalar dstab[5][3], U_dot[5] = {0}; 439 for (CeedInt j=0; j<5; j++) U_dot[j] = context->ijacobian_time_shift * dU[j]; 440 Stabilization(context, s, Tau_d, grad_ds, U_dot, dbody_force, x_i, dstab); 441 442 for (int j=0; j<5; j++) 443 for (int k=0; k<3; k++) 444 Grad_v[k][j][i] += wdetJ*(dstab[j][0] * dXdx[k][0] + 445 dstab[j][1] * dXdx[k][1] + 446 dstab[j][2] * dXdx[k][2]); 447 448 } // End Quadrature Point Loop 449 return 0; 450 } 451 452 CEED_QFUNCTION(IJacobian_Newtonian_Conserv)(void *ctx, CeedInt Q, 453 const CeedScalar *const *in, CeedScalar *const *out) { 454 return IJacobian_Newtonian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 455 } 456 457 CEED_QFUNCTION(IJacobian_Newtonian_Prim)(void *ctx, CeedInt Q, 458 const CeedScalar *const *in, CeedScalar *const *out) { 459 return IJacobian_Newtonian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 460 } 461 462 // ***************************************************************************** 463 // Compute boundary integral (ie. for strongly set inflows) 464 // ***************************************************************************** 465 CEED_QFUNCTION_HELPER int BoundaryIntegral(void *ctx, CeedInt Q, 466 const CeedScalar *const *in, CeedScalar *const *out, 467 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 468 469 //*INDENT-OFF* 470 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 471 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 472 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 473 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 474 475 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[0], 476 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA]) out[1]; 477 478 //*INDENT-ON* 479 480 const NewtonianIdealGasContext context = (NewtonianIdealGasContext) ctx; 481 const bool is_implicit = context->is_implicit; 482 483 CeedPragmaSIMD 484 for(CeedInt i=0; i<Q; i++) { 485 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 486 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 487 State s = StateFromQi(context, qi, x_i); 488 489 const CeedScalar wdetJb = (is_implicit ? -1. : 1.) * q_data_sur[0][i]; 490 // ---- Normal vector 491 const CeedScalar norm[3] = {q_data_sur[1][i], 492 q_data_sur[2][i], 493 q_data_sur[3][i] 494 }; 495 496 const CeedScalar dXdx[2][3] = { 497 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 498 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 499 }; 500 501 State grad_s[3]; 502 for (CeedInt j=0; j<3; j++) { 503 CeedScalar dx_i[3] = {0}, dqi[5]; 504 for (CeedInt k=0; k<5; k++) 505 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 506 Grad_q[1][k][i] * dXdx[1][j]; 507 dx_i[j] = 1.; 508 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 509 } 510 511 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 512 KMStrainRate(grad_s, strain_rate); 513 NewtonianStress(context, strain_rate, kmstress); 514 KMUnpack(kmstress, stress); 515 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 516 517 StateConservative F_inviscid[3]; 518 FluxInviscid(context, s, F_inviscid); 519 520 CeedScalar Flux[5]; 521 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 522 523 for (CeedInt j=0; j<5; j++) v[j][i] = -wdetJb * Flux[j]; 524 525 for (int j=0; j<5; j++) jac_data_sur[j][i] = qi[j]; 526 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 527 } 528 return 0; 529 } 530 531 CEED_QFUNCTION(BoundaryIntegral_Conserv)(void *ctx, CeedInt Q, 532 const CeedScalar *const *in, CeedScalar *const *out) { 533 return BoundaryIntegral(ctx, Q, in, out, StateFromU, StateFromU_fwd); 534 } 535 536 CEED_QFUNCTION(BoundaryIntegral_Prim)(void *ctx, CeedInt Q, 537 const CeedScalar *const *in, CeedScalar *const *out) { 538 return BoundaryIntegral(ctx, Q, in, out, StateFromY, StateFromY_fwd); 539 } 540 541 // ***************************************************************************** 542 // Jacobian for "set nothing" boundary integral 543 // ***************************************************************************** 544 CEED_QFUNCTION_HELPER int BoundaryIntegral_Jacobian(void *ctx, CeedInt Q, 545 const CeedScalar *const *in, CeedScalar *const *out, 546 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 547 // *INDENT-OFF* 548 // Inputs 549 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 550 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 551 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 552 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 553 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 554 // Outputs 555 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 556 // *INDENT-ON* 557 558 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 559 const bool implicit = context->is_implicit; 560 561 CeedPragmaSIMD 562 // Quadrature Point Loop 563 for (CeedInt i=0; i<Q; i++) { 564 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 565 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 566 const CeedScalar norm[3] = {q_data_sur[1][i], 567 q_data_sur[2][i], 568 q_data_sur[3][i] 569 }; 570 const CeedScalar dXdx[2][3] = { 571 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 572 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 573 }; 574 575 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 576 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 577 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 578 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 579 580 State s = StateFromQi(context, qi, x_i); 581 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 582 583 State grad_ds[3]; 584 for (CeedInt j=0; j<3; j++) { 585 CeedScalar dx_i[3] = {0}, dqi_j[5]; 586 for (CeedInt k=0; k<5; k++) 587 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 588 Grad_dq[1][k][i] * dXdx[1][j]; 589 dx_i[j] = 1.; 590 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 591 } 592 593 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 594 KMStrainRate(grad_ds, dstrain_rate); 595 NewtonianStress(context, dstrain_rate, dkmstress); 596 KMUnpack(dkmstress, dstress); 597 KMUnpack(kmstress, stress); 598 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 599 600 StateConservative dF_inviscid[3]; 601 FluxInviscid_fwd(context, s, ds, dF_inviscid); 602 603 CeedScalar dFlux[5]; 604 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 605 606 for (int j=0; j<5; j++) v[j][i] = -wdetJb * dFlux[j]; 607 } // End Quadrature Point Loop 608 return 0; 609 } 610 611 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Conserv)(void *ctx, CeedInt Q, 612 const CeedScalar *const *in, CeedScalar *const *out) { 613 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 614 } 615 616 CEED_QFUNCTION(BoundaryIntegral_Jacobian_Prim)(void *ctx, CeedInt Q, 617 const CeedScalar *const *in, CeedScalar *const *out) { 618 return BoundaryIntegral_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 619 } 620 621 // ***************************************************************************** 622 // Outflow boundary condition, weakly setting a constant pressure 623 // ***************************************************************************** 624 CEED_QFUNCTION_HELPER int PressureOutflow(void *ctx, CeedInt Q, 625 const CeedScalar *const *in, CeedScalar *const *out, 626 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 627 // *INDENT-OFF* 628 // Inputs 629 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 630 (*Grad_q)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 631 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 632 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 633 // Outputs 634 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 635 (*jac_data_sur)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 636 // *INDENT-ON* 637 638 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 639 const bool implicit = context->is_implicit; 640 const CeedScalar P0 = context->P0; 641 642 CeedPragmaSIMD 643 // Quadrature Point Loop 644 for (CeedInt i=0; i<Q; i++) { 645 // Setup 646 // -- Interp in 647 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 648 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 649 State s = StateFromQi(context, qi, x_i); 650 s.Y.pressure = P0; 651 652 // -- Interp-to-Interp q_data 653 // For explicit mode, the surface integral is on the RHS of ODE q_dot = f(q). 654 // For implicit mode, it gets pulled to the LHS of implicit ODE/DAE g(q_dot, q). 655 // We can effect this by swapping the sign on this weight 656 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 657 658 // ---- Normal vector 659 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 660 661 const CeedScalar dXdx[2][3] = { 662 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 663 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 664 }; 665 666 State grad_s[3]; 667 for (CeedInt j=0; j<3; j++) { 668 CeedScalar dx_i[3] = {0}, dqi[5]; 669 for (CeedInt k=0; k<5; k++) 670 dqi[k] = Grad_q[0][k][i] * dXdx[0][j] + 671 Grad_q[1][k][i] * dXdx[1][j]; 672 dx_i[j] = 1.; 673 grad_s[j] = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 674 } 675 676 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 677 KMStrainRate(grad_s, strain_rate); 678 NewtonianStress(context, strain_rate, kmstress); 679 KMUnpack(kmstress, stress); 680 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 681 682 StateConservative F_inviscid[3]; 683 FluxInviscid(context, s, F_inviscid); 684 685 CeedScalar Flux[5]; 686 FluxTotal_Boundary(F_inviscid, stress, Fe, norm, Flux); 687 688 for (CeedInt j=0; j<5; j++) v[j][i] = -wdetJb * Flux[j]; 689 690 // Save values for Jacobian 691 for (int j=0; j<5; j++) jac_data_sur[j][i] = qi[j]; 692 for (int j=0; j<6; j++) jac_data_sur[5+j][i] = kmstress[j]; 693 } // End Quadrature Point Loop 694 return 0; 695 } 696 697 CEED_QFUNCTION(PressureOutflow_Conserv)(void *ctx, CeedInt Q, 698 const CeedScalar *const *in, CeedScalar *const *out) { 699 return PressureOutflow(ctx, Q, in, out, StateFromU, StateFromU_fwd); 700 } 701 702 CEED_QFUNCTION(PressureOutflow_Prim)(void *ctx, CeedInt Q, 703 const CeedScalar *const *in, CeedScalar *const *out) { 704 return PressureOutflow(ctx, Q, in, out, StateFromY, StateFromY_fwd); 705 } 706 707 // ***************************************************************************** 708 // Jacobian for weak-pressure outflow boundary condition 709 // ***************************************************************************** 710 CEED_QFUNCTION_HELPER int PressureOutflow_Jacobian(void *ctx, CeedInt Q, 711 const CeedScalar *const *in, CeedScalar *const *out, 712 StateFromQi_t StateFromQi, StateFromQi_fwd_t StateFromQi_fwd) { 713 // *INDENT-OFF* 714 // Inputs 715 const CeedScalar (*dq)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 716 (*Grad_dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 717 (*q_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 718 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 719 (*jac_data_sur)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 720 // Outputs 721 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 722 // *INDENT-ON* 723 724 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 725 const bool implicit = context->is_implicit; 726 727 CeedPragmaSIMD 728 // Quadrature Point Loop 729 for (CeedInt i=0; i<Q; i++) { 730 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 731 const CeedScalar wdetJb = (implicit ? -1. : 1.) * q_data_sur[0][i]; 732 const CeedScalar norm[3] = {q_data_sur[1][i], q_data_sur[2][i], q_data_sur[3][i]}; 733 const CeedScalar dXdx[2][3] = { 734 {q_data_sur[4][i], q_data_sur[5][i], q_data_sur[6][i]}, 735 {q_data_sur[7][i], q_data_sur[8][i], q_data_sur[9][i]} 736 }; 737 738 CeedScalar qi[5], kmstress[6], dqi[5], dx_i[3] = {0.}; 739 for (int j=0; j<5; j++) qi[j] = jac_data_sur[j][i]; 740 for (int j=0; j<6; j++) kmstress[j] = jac_data_sur[5+j][i]; 741 for (int j=0; j<5; j++) dqi[j] = dq[j][i]; 742 743 State s = StateFromQi(context, qi, x_i); 744 State ds = StateFromQi_fwd(context, s, dqi, x_i, dx_i); 745 s.Y.pressure = context->P0; 746 ds.Y.pressure = 0.; 747 748 State grad_ds[3]; 749 for (CeedInt j=0; j<3; j++) { 750 CeedScalar dx_i[3] = {0}, dqi_j[5]; 751 for (CeedInt k=0; k<5; k++) 752 dqi_j[k] = Grad_dq[0][k][i] * dXdx[0][j] + 753 Grad_dq[1][k][i] * dXdx[1][j]; 754 dx_i[j] = 1.; 755 grad_ds[j] = StateFromQi_fwd(context, s, dqi_j, x_i, dx_i); 756 } 757 758 CeedScalar dstrain_rate[6], dkmstress[6], stress[3][3], dstress[3][3], dFe[3]; 759 KMStrainRate(grad_ds, dstrain_rate); 760 NewtonianStress(context, dstrain_rate, dkmstress); 761 KMUnpack(dkmstress, dstress); 762 KMUnpack(kmstress, stress); 763 ViscousEnergyFlux_fwd(context, s.Y, ds.Y, grad_ds, stress, dstress, dFe); 764 765 StateConservative dF_inviscid[3]; 766 FluxInviscid_fwd(context, s, ds, dF_inviscid); 767 768 CeedScalar dFlux[5]; 769 FluxTotal_Boundary(dF_inviscid, dstress, dFe, norm, dFlux); 770 771 for (int j=0; j<5; j++) v[j][i] = -wdetJb * dFlux[j]; 772 } // End Quadrature Point Loop 773 return 0; 774 } 775 776 CEED_QFUNCTION(PressureOutflow_Jacobian_Conserv)(void *ctx, CeedInt Q, 777 const CeedScalar *const *in, CeedScalar *const *out) { 778 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromU, StateFromU_fwd); 779 } 780 781 CEED_QFUNCTION(PressureOutflow_Jacobian_Prim)(void *ctx, CeedInt Q, 782 const CeedScalar *const *in, CeedScalar *const *out) { 783 return PressureOutflow_Jacobian(ctx, Q, in, out, StateFromY, StateFromY_fwd); 784 } 785 786 #endif // newtonian_h 787