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