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