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 <math.h> 16 #include <ceed.h> 17 #include "newtonian_types.h" 18 19 #ifndef M_PI 20 #define M_PI 3.14159265358979323846 21 #endif 22 23 typedef struct { 24 CeedScalar pressure; 25 CeedScalar velocity[3]; 26 CeedScalar temperature; 27 } StatePrimitive; 28 29 typedef struct { 30 CeedScalar density; 31 CeedScalar momentum[3]; 32 CeedScalar E_total; 33 } StateConservative; 34 35 typedef struct { 36 StateConservative U; 37 StatePrimitive Y; 38 } State; 39 40 CEED_QFUNCTION_HELPER CeedScalar Dot3(const CeedScalar u[3], 41 const CeedScalar v[3]) { 42 return u[0]*v[0] + u[1]*v[1] + u[2]*v[2]; 43 } 44 45 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative( 46 NewtonianIdealGasContext gas, StateConservative U, const CeedScalar x[3]) { 47 StatePrimitive Y; 48 for (int i=0; i<3; i++) Y.velocity[i] = U.momentum[i] / U.density; 49 CeedScalar e_kinetic = .5 * Dot3(Y.velocity, Y.velocity); 50 CeedScalar e_potential = -Dot3(gas->g, x); 51 CeedScalar e_total = U.E_total / U.density; 52 CeedScalar e_internal = e_total - e_kinetic - e_potential; 53 Y.temperature = e_internal / gas->cv; 54 Y.pressure = (gas->cp / gas->cv - 1) * U.density * e_internal; 55 return Y; 56 } 57 58 CEED_QFUNCTION_HELPER StatePrimitive StatePrimitiveFromConservative_fwd( 59 NewtonianIdealGasContext gas, State s, StateConservative dU, 60 const CeedScalar x[3], const CeedScalar dx[3]) { 61 StatePrimitive dY; 62 for (int i=0; i<3; i++) { 63 dY.velocity[i] = (dU.momentum[i] - s.Y.velocity[i] * dU.density) / s.U.density; 64 } 65 CeedScalar e_kinetic = .5 * Dot3(s.Y.velocity, s.Y.velocity); 66 CeedScalar de_kinetic = Dot3(dY.velocity, s.Y.velocity); 67 CeedScalar e_potential = -Dot3(gas->g, x); 68 CeedScalar de_potential = -Dot3(gas->g, dx); 69 CeedScalar e_total = s.U.E_total / s.U.density; 70 CeedScalar de_total = (dU.E_total - e_total * dU.density) / s.U.density; 71 CeedScalar e_internal = e_total - e_kinetic - e_potential; 72 CeedScalar de_internal = de_total - de_kinetic - de_potential; 73 dY.temperature = de_internal / gas->cv; 74 dY.pressure = (gas->cp / gas->cv - 1) 75 * (dU.density * e_internal + s.U.density * de_internal); 76 return dY; 77 } 78 79 CEED_QFUNCTION_HELPER State StateFromU(NewtonianIdealGasContext gas, 80 const CeedScalar U[5], const CeedScalar x[3]) { 81 State s; 82 s.U.density = U[0]; 83 s.U.momentum[0] = U[1]; 84 s.U.momentum[1] = U[2]; 85 s.U.momentum[2] = U[3]; 86 s.U.E_total = U[4]; 87 s.Y = StatePrimitiveFromConservative(gas, s.U, x); 88 return s; 89 } 90 91 CEED_QFUNCTION_HELPER void FluxInviscid(NewtonianIdealGasContext gas, State s, 92 StateConservative Flux[3]) { 93 for (int i=0; i<3; i++) { 94 Flux[i].density = s.U.momentum[i]; 95 for (int j=0; j<3; j++) 96 Flux[i].momentum[j] = s.U.momentum[i] * s.Y.velocity[j] 97 + s.Y.pressure * (i == j); 98 Flux[i].E_total = (s.U.E_total + s.Y.pressure) * s.Y.velocity[i]; 99 } 100 } 101 102 CEED_QFUNCTION_HELPER void FluxInviscid_fwd(NewtonianIdealGasContext gas, 103 State s, State ds, StateConservative dFlux[3]) { 104 for (int i=0; i<3; i++) { 105 dFlux[i].density = ds.U.momentum[i]; 106 for (int j=0; j<3; j++) 107 dFlux[i].momentum[j] = ds.U.momentum[i] * s.Y.velocity[j] + 108 s.U.momentum[i] * ds.Y.velocity[j] + ds.Y.pressure * (i == j); 109 dFlux[i].E_total = (ds.U.E_total + ds.Y.pressure) * s.Y.velocity[i] + 110 (s.U.E_total + s.Y.pressure) * ds.Y.velocity[i]; 111 } 112 } 113 114 // Kelvin-Mandel notation 115 CEED_QFUNCTION_HELPER void KMStrainRate(const State grad_s[3], 116 CeedScalar strain_rate[6]) { 117 const CeedScalar weight = 1 / sqrt(2.); 118 strain_rate[0] = grad_s[0].Y.velocity[0]; 119 strain_rate[1] = grad_s[1].Y.velocity[1]; 120 strain_rate[2] = grad_s[2].Y.velocity[2]; 121 strain_rate[3] = weight * (grad_s[2].Y.velocity[1] + grad_s[1].Y.velocity[2]); 122 strain_rate[4] = weight * (grad_s[2].Y.velocity[0] + grad_s[0].Y.velocity[2]); 123 strain_rate[5] = weight * (grad_s[1].Y.velocity[0] + grad_s[0].Y.velocity[1]); 124 } 125 126 CEED_QFUNCTION_HELPER void KMUnpack(const CeedScalar v[6], CeedScalar A[3][3]) { 127 const CeedScalar weight = 1 / sqrt(2.); 128 A[0][0] = v[0]; 129 A[1][1] = v[1]; 130 A[2][2] = v[2]; 131 A[2][1] = A[1][2] = weight * v[3]; 132 A[2][0] = A[0][2] = weight * v[4]; 133 A[1][0] = A[0][1] = weight * v[5]; 134 } 135 136 CEED_QFUNCTION_HELPER void NewtonianStress(NewtonianIdealGasContext gas, 137 const CeedScalar strain_rate[6], CeedScalar stress[6]) { 138 CeedScalar div_u = strain_rate[0] + strain_rate[1] + strain_rate[2]; 139 for (int i=0; i<6; i++) { 140 stress[i] = gas->mu * (2 * strain_rate[i] + gas->lambda * div_u * (i < 3)); 141 } 142 } 143 144 CEED_QFUNCTION_HELPER void ViscousEnergyFlux(NewtonianIdealGasContext gas, 145 StatePrimitive Y, const State grad_s[3], const CeedScalar stress[3][3], 146 CeedScalar Fe[3]) { 147 for (int i=0; i<3; i++) { 148 Fe[i] = - Y.velocity[0] * stress[0][i] 149 - Y.velocity[1] * stress[1][i] 150 - Y.velocity[2] * stress[2][i] 151 - gas->k * grad_s[i].Y.temperature; 152 } 153 } 154 155 // ***************************************************************************** 156 // Helper function for computing flux Jacobian 157 // ***************************************************************************** 158 CEED_QFUNCTION_HELPER void computeFluxJacobian_NS(CeedScalar dF[3][5][5], 159 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 160 const CeedScalar gamma, const CeedScalar g[3], const CeedScalar x[3]) { 161 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 162 CeedScalar e_potential = -(g[0]*x[0] + g[1]*x[1] + g[2]*x[2]); 163 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 164 for (CeedInt j=0; j<3; j++) { // Rows of each Jacobian matrix 165 dF[i][j+1][0] = ((i==j) ? ((gamma-1.)*(u_sq/2. - e_potential)) : 0.) - 166 u[i]*u[j]; 167 for (CeedInt k=0; k<3; k++) { // Columns of each Jacobian matrix 168 dF[i][0][k+1] = ((i==k) ? 1. : 0.); 169 dF[i][j+1][k+1] = ((j==k) ? u[i] : 0.) + 170 ((i==k) ? u[j] : 0.) - 171 ((i==j) ? u[k] : 0.) * (gamma-1.); 172 dF[i][4][k+1] = ((i==k) ? (E*gamma/rho - (gamma-1.)*u_sq/2.) : 0.) - 173 (gamma-1.)*u[i]*u[k]; 174 } 175 dF[i][j+1][4] = ((i==j) ? (gamma-1.) : 0.); 176 } 177 dF[i][4][0] = u[i] * ((gamma-1.)*u_sq - E*gamma/rho); 178 dF[i][4][4] = u[i] * gamma; 179 } 180 } 181 182 // ***************************************************************************** 183 // Helper function for computing flux Jacobian of Primitive variables 184 // ***************************************************************************** 185 CEED_QFUNCTION_HELPER void computeFluxJacobian_NSp(CeedScalar dF[3][5][5], 186 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 187 const CeedScalar Rd, const CeedScalar cv) { 188 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 189 // TODO Add in gravity's contribution 190 191 CeedScalar T = ( E / rho - u_sq / 2. ) / cv; 192 CeedScalar drdT = -rho / T; 193 CeedScalar drdP = 1. / ( Rd * T); 194 CeedScalar etot = E / rho ; 195 CeedScalar e2p = drdP * etot + 1. ; 196 CeedScalar e3p = ( E + rho * Rd * T ); 197 CeedScalar e4p = drdT * etot + rho * cv ; 198 199 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 200 for (CeedInt j=0; j<3; j++) { // j counts F^{m_j} 201 // [row][col] of A_i 202 dF[i][j+1][0] = drdP * u[i] * u[j] + ((i==j) ? 1. : 0.); // F^{{m_j} wrt p 203 for (CeedInt k=0; k<3; k++) { // k counts the wrt vel_k 204 dF[i][0][k+1] = ((i==k) ? rho : 0.); // F^c wrt u_k 205 dF[i][j+1][k+1] = (((j==k) ? u[i] : 0.) + // F^m_j wrt u_k 206 ((i==k) ? u[j] : 0.) ) * rho; 207 dF[i][4][k+1] = rho * u[i] * u[k] 208 + ((i==k) ? e3p : 0.) ; // F^e wrt u_k 209 } 210 dF[i][j+1][4] = drdT * u[i] * u[j]; // F^{m_j} wrt T 211 } 212 dF[i][4][0] = u[i] * e2p; // F^e wrt p 213 dF[i][4][4] = u[i] * e4p; // F^e wrt T 214 dF[i][0][0] = u[i] * drdP; // F^c wrt p 215 dF[i][0][4] = u[i] * drdT; // F^c wrt T 216 } 217 } 218 219 CEED_QFUNCTION_HELPER void PrimitiveToConservative_fwd(const CeedScalar rho, 220 const CeedScalar u[3], const CeedScalar E, const CeedScalar Rd, 221 const CeedScalar cv, const CeedScalar dY[5], CeedScalar dU[5]) { 222 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; 223 CeedScalar T = ( E / rho - u_sq / 2. ) / cv; 224 CeedScalar drdT = -rho / T; 225 CeedScalar drdP = 1. / ( Rd * T); 226 dU[0] = drdP * dY[0] + drdT * dY[4]; 227 CeedScalar de_kinetic = 0; 228 for (int i=0; i<3; i++) { 229 dU[1+i] = dU[0] * u[i] + rho * dY[1+i]; 230 de_kinetic += u[i] * dY[1+i]; 231 } 232 dU[4] = rho * cv * dY[4] + dU[0] * cv * T // internal energy: rho * e 233 + rho * de_kinetic + .5 * dU[0] * u_sq; // kinetic energy: .5 * rho * |u|^2 234 } 235 236 // ***************************************************************************** 237 // Helper function for computing Tau elements (stabilization constant) 238 // Model from: 239 // PHASTA 240 // 241 // Tau[i] = itau=0 which is diagonal-Shakib (3 values still but not spatial) 242 // 243 // Where NOT UPDATED YET 244 // ***************************************************************************** 245 CEED_QFUNCTION_HELPER void Tau_diagPrim(CeedScalar Tau_d[3], 246 const CeedScalar dXdx[3][3], const CeedScalar u[3], 247 const CeedScalar cv, const NewtonianIdealGasContext newt_ctx, 248 const CeedScalar mu, const CeedScalar dt, 249 const CeedScalar rho) { 250 // Context 251 const CeedScalar Ctau_t = newt_ctx->Ctau_t; 252 const CeedScalar Ctau_v = newt_ctx->Ctau_v; 253 const CeedScalar Ctau_C = newt_ctx->Ctau_C; 254 const CeedScalar Ctau_M = newt_ctx->Ctau_M; 255 const CeedScalar Ctau_E = newt_ctx->Ctau_E; 256 CeedScalar gijd[6]; 257 CeedScalar tau; 258 CeedScalar dts; 259 CeedScalar fact; 260 261 //*INDENT-OFF* 262 gijd[0] = dXdx[0][0] * dXdx[0][0] 263 + dXdx[1][0] * dXdx[1][0] 264 + dXdx[2][0] * dXdx[2][0]; 265 266 gijd[1] = dXdx[0][0] * dXdx[0][1] 267 + dXdx[1][0] * dXdx[1][1] 268 + dXdx[2][0] * dXdx[2][1]; 269 270 gijd[2] = dXdx[0][1] * dXdx[0][1] 271 + dXdx[1][1] * dXdx[1][1] 272 + dXdx[2][1] * dXdx[2][1]; 273 274 gijd[3] = dXdx[0][0] * dXdx[0][2] 275 + dXdx[1][0] * dXdx[1][2] 276 + dXdx[2][0] * dXdx[2][2]; 277 278 gijd[4] = dXdx[0][1] * dXdx[0][2] 279 + dXdx[1][1] * dXdx[1][2] 280 + dXdx[2][1] * dXdx[2][2]; 281 282 gijd[5] = dXdx[0][2] * dXdx[0][2] 283 + dXdx[1][2] * dXdx[1][2] 284 + dXdx[2][2] * dXdx[2][2]; 285 //*INDENT-ON* 286 287 dts = Ctau_t / dt ; 288 289 tau = rho*rho*((4. * dts * dts) 290 + u[0] * ( u[0] * gijd[0] + 2. * ( u[1] * gijd[1] + u[2] * gijd[3])) 291 + u[1] * ( u[1] * gijd[2] + 2. * u[2] * gijd[4]) 292 + u[2] * u[2] * gijd[5]) 293 + Ctau_v* mu * mu * 294 (gijd[0]*gijd[0] + gijd[2]*gijd[2] + gijd[5]*gijd[5] + 295 + 2. * (gijd[1]*gijd[1] + gijd[3]*gijd[3] + gijd[4]*gijd[4])); 296 297 fact=sqrt(tau); 298 299 Tau_d[0] = Ctau_C * fact / (rho*(gijd[0] + gijd[2] + gijd[5]))*0.125; 300 301 Tau_d[1] = Ctau_M / fact; 302 Tau_d[2] = Ctau_E / ( fact * cv ); 303 304 // consider putting back the way I initially had it Ctau_E * Tau_d[1] /cv 305 // to avoid a division if the compiler is smart enough to see that cv IS 306 // a constant that it could invert once for all elements 307 // but in that case energy tau is scaled by the product of Ctau_E * Ctau_M 308 // OR we could absorb cv into Ctau_E but this puts more burden on user to 309 // know how to change constants with a change of fluid or units. Same for 310 // Ctau_v * mu * mu IF AND ONLY IF we don't add viscosity law =f(T) 311 } 312 313 // ***************************************************************************** 314 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 315 // ***************************************************************************** 316 CEED_QFUNCTION(ICsNewtonianIG)(void *ctx, CeedInt Q, 317 const CeedScalar *const *in, CeedScalar *const *out) { 318 // Inputs 319 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 320 321 // Outputs 322 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 323 324 // Context 325 const SetupContext context = (SetupContext)ctx; 326 const CeedScalar theta0 = context->theta0; 327 const CeedScalar P0 = context->P0; 328 const CeedScalar cv = context->cv; 329 const CeedScalar cp = context->cp; 330 const CeedScalar *g = context->g; 331 const CeedScalar Rd = cp - cv; 332 333 // Quadrature Point Loop 334 CeedPragmaSIMD 335 for (CeedInt i=0; i<Q; i++) { 336 CeedScalar q[5] = {0.}; 337 338 // Setup 339 // -- Coordinates 340 const CeedScalar x[3] = {X[0][i], X[1][i], X[2][i]}; 341 const CeedScalar e_potential = -(g[0]*x[0] + g[1]*x[1] + g[2]*x[2]); 342 343 // -- Density 344 const CeedScalar rho = P0 / (Rd*theta0); 345 346 // Initial Conditions 347 q[0] = rho; 348 q[1] = 0.0; 349 q[2] = 0.0; 350 q[3] = 0.0; 351 q[4] = rho * (cv*theta0 + e_potential); 352 353 for (CeedInt j=0; j<5; j++) 354 q0[j][i] = q[j]; 355 } // End of Quadrature Point Loop 356 return 0; 357 } 358 359 // ***************************************************************************** 360 // This QFunction implements the following formulation of Navier-Stokes with 361 // explicit time stepping method 362 // 363 // This is 3D compressible Navier-Stokes in conservation form with state 364 // variables of density, momentum density, and total energy density. 365 // 366 // State Variables: q = ( rho, U1, U2, U3, E ) 367 // rho - Mass Density 368 // Ui - Momentum Density, Ui = rho ui 369 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 370 // 371 // Navier-Stokes Equations: 372 // drho/dt + div( U ) = 0 373 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 374 // dE/dt + div( (E + P) u ) = div( Fe ) 375 // 376 // Viscous Stress: 377 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 378 // 379 // Thermal Stress: 380 // Fe = u Fu + k grad( T ) 381 // Equation of State 382 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 383 // 384 // Stabilization: 385 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 386 // f1 = rho sqrt(ui uj gij) 387 // gij = dXi/dX * dXi/dX 388 // TauC = Cc f1 / (8 gii) 389 // TauM = min( 1 , 1 / f1 ) 390 // TauE = TauM / (Ce cv) 391 // 392 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 393 // 394 // Constants: 395 // lambda = - 2 / 3, From Stokes hypothesis 396 // mu , Dynamic viscosity 397 // k , Thermal conductivity 398 // cv , Specific heat, constant volume 399 // cp , Specific heat, constant pressure 400 // g , Gravity 401 // gamma = cp / cv, Specific heat ratio 402 // 403 // We require the product of the inverse of the Jacobian (dXdx_j,k) and 404 // its transpose (dXdx_k,j) to properly compute integrals of the form: 405 // int( gradv gradu ) 406 // 407 // ***************************************************************************** 408 CEED_QFUNCTION(RHSFunction_Newtonian)(void *ctx, CeedInt Q, 409 const CeedScalar *const *in, CeedScalar *const *out) { 410 // *INDENT-OFF* 411 // Inputs 412 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 413 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 414 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 415 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 416 // Outputs 417 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 418 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 419 // *INDENT-ON* 420 421 // Context 422 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 423 const CeedScalar mu = context->mu; 424 const CeedScalar cv = context->cv; 425 const CeedScalar cp = context->cp; 426 const CeedScalar *g = context->g; 427 const CeedScalar dt = context->dt; 428 const CeedScalar gamma = cp / cv; 429 const CeedScalar Rd = cp - cv; 430 431 CeedPragmaSIMD 432 // Quadrature Point Loop 433 for (CeedInt i=0; i<Q; i++) { 434 CeedScalar U[5]; 435 for (int j=0; j<5; j++) U[j] = q[j][i]; 436 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 437 State s = StateFromU(context, U, x_i); 438 439 // -- Interp-to-Interp q_data 440 const CeedScalar wdetJ = q_data[0][i]; 441 // -- Interp-to-Grad q_data 442 // ---- Inverse of change of coordinate matrix: X_i,j 443 // *INDENT-OFF* 444 const CeedScalar dXdx[3][3] = {{q_data[1][i], 445 q_data[2][i], 446 q_data[3][i]}, 447 {q_data[4][i], 448 q_data[5][i], 449 q_data[6][i]}, 450 {q_data[7][i], 451 q_data[8][i], 452 q_data[9][i]} 453 }; 454 // *INDENT-ON* 455 456 State grad_s[3]; 457 for (int j=0; j<3; j++) { 458 CeedScalar dx_i[3] = {0}; 459 grad_s[j].U.density = dq[0][0][i] * dXdx[0][j] 460 + dq[1][0][i] * dXdx[1][j] + dq[2][0][i] * dXdx[2][j]; 461 for (int k=0; k<3; k++) grad_s[j].U.momentum[k] = dq[0][k+1][i] * dXdx[0][j] 462 + dq[1][k+1][i] * dXdx[1][j] + dq[2][k+1][i] * dXdx[2][j]; 463 grad_s[j].U.E_total = dq[0][4][i] * dXdx[0][j] + dq[1][4][i] * dXdx[1][j] + 464 dq[2][4][i] * dXdx[2][j]; 465 dx_i[j] = 1.; 466 grad_s[j].Y = StatePrimitiveFromConservative_fwd(context, s, grad_s[j].U, 467 x_i, dx_i); 468 } 469 470 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 471 KMStrainRate(grad_s, strain_rate); 472 NewtonianStress(context, strain_rate, kmstress); 473 KMUnpack(kmstress, stress); 474 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 475 476 StateConservative F_inviscid[3]; 477 FluxInviscid(context, s, F_inviscid); 478 479 // Total flux 480 CeedScalar Flux[5][3]; 481 for (int j=0; j<3; j++) { 482 Flux[0][j] = F_inviscid[j].density; 483 for (int k=0; k<3; k++) 484 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 485 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 486 } 487 488 for (int j=0; j<3; j++) { 489 for (int k=0; k<5; k++) { 490 dv[j][k][i] = wdetJ * (dXdx[j][0] * Flux[k][0] + 491 dXdx[j][1] * Flux[k][1] + 492 dXdx[j][2] * Flux[k][2]); 493 } 494 } 495 496 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 497 for (int j=0; j<5; j++) 498 v[j][i] = wdetJ * body_force[j]; 499 500 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 501 CeedScalar jacob_F_conv[3][5][5] = {0}; 502 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 503 gamma, g, x_i); 504 CeedScalar grad_U[5][3]; 505 for (int j=0; j<3; j++) { 506 grad_U[0][j] = grad_s[j].U.density; 507 for (int k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 508 grad_U[4][j] = grad_s[j].U.E_total; 509 } 510 511 // strong_conv = dF/dq * dq/dx (Strong convection) 512 CeedScalar strong_conv[5] = {0}; 513 for (int j=0; j<3; j++) 514 for (int k=0; k<5; k++) 515 for (int l=0; l<5; l++) 516 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 517 518 // -- Stabilization method: none, SU, or SUPG 519 CeedScalar stab[5][3] = {{0.}}; 520 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 521 CeedScalar Tau_d[3] = {0.}; 522 switch (context->stabilization) { 523 case STAB_NONE: // Galerkin 524 break; 525 case STAB_SU: // SU 526 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 527 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 528 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 529 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 530 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 531 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 532 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 533 tau_strong_conv, 534 tau_strong_conv_conservative); 535 for (int j=0; j<3; j++) 536 for (int k=0; k<5; k++) 537 for (int l=0; l<5; l++) 538 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 539 540 for (int j=0; j<5; j++) 541 for (int k=0; k<3; k++) 542 dv[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 543 stab[j][1] * dXdx[k][1] + 544 stab[j][2] * dXdx[k][2]); 545 break; 546 case STAB_SUPG: // SUPG is not implemented for explicit scheme 547 break; 548 } 549 550 } // End Quadrature Point Loop 551 552 // Return 553 return 0; 554 } 555 556 // ***************************************************************************** 557 // This QFunction implements the Navier-Stokes equations (mentioned above) with 558 // implicit time stepping method 559 // 560 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 561 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 562 // (diffussive terms will be added later) 563 // 564 // ***************************************************************************** 565 CEED_QFUNCTION(IFunction_Newtonian)(void *ctx, CeedInt Q, 566 const CeedScalar *const *in, 567 CeedScalar *const *out) { 568 // *INDENT-OFF* 569 // Inputs 570 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 571 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 572 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 573 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 574 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 575 // Outputs 576 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 577 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 578 // *INDENT-ON* 579 // Context 580 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 581 const CeedScalar mu = context->mu; 582 const CeedScalar cv = context->cv; 583 const CeedScalar cp = context->cp; 584 const CeedScalar *g = context->g; 585 const CeedScalar dt = context->dt; 586 const CeedScalar gamma = cp / cv; 587 const CeedScalar Rd = cp-cv; 588 589 CeedPragmaSIMD 590 // Quadrature Point Loop 591 for (CeedInt i=0; i<Q; i++) { 592 CeedScalar U[5]; 593 for (int j=0; j<5; j++) U[j] = q[j][i]; 594 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 595 State s = StateFromU(context, U, x_i); 596 597 // -- Interp-to-Interp q_data 598 const CeedScalar wdetJ = q_data[0][i]; 599 // -- Interp-to-Grad q_data 600 // ---- Inverse of change of coordinate matrix: X_i,j 601 // *INDENT-OFF* 602 const CeedScalar dXdx[3][3] = {{q_data[1][i], 603 q_data[2][i], 604 q_data[3][i]}, 605 {q_data[4][i], 606 q_data[5][i], 607 q_data[6][i]}, 608 {q_data[7][i], 609 q_data[8][i], 610 q_data[9][i]} 611 }; 612 // *INDENT-ON* 613 State grad_s[3]; 614 for (int j=0; j<3; j++) { 615 CeedScalar dx_i[3]; 616 grad_s[j].U.density = dq[0][0][i] * dXdx[0][j] 617 + dq[1][0][i] * dXdx[1][j] + dq[2][0][i] * dXdx[2][j]; 618 for (int k=0; k<3; k++) grad_s[j].U.momentum[k] = dq[0][k+1][i] * dXdx[0][j] 619 + dq[1][k+1][i] * dXdx[1][j] + dq[2][k+1][i] * dXdx[2][j]; 620 grad_s[j].U.E_total = dq[0][4][i] * dXdx[0][j] + dq[1][4][i] * dXdx[1][j] + 621 dq[2][4][i] * dXdx[2][j]; 622 dx_i[j] = 1.; 623 grad_s[j].Y = StatePrimitiveFromConservative_fwd(context, s, grad_s[j].U, 624 x_i, dx_i); 625 } 626 627 CeedScalar strain_rate[6], kmstress[6], stress[3][3], Fe[3]; 628 KMStrainRate(grad_s, strain_rate); 629 NewtonianStress(context, strain_rate, kmstress); 630 KMUnpack(kmstress, stress); 631 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 632 633 StateConservative F_inviscid[3]; 634 FluxInviscid(context, s, F_inviscid); 635 636 637 // Total flux 638 CeedScalar Flux[5][3]; 639 for (int j=0; j<3; j++) { 640 Flux[0][j] = F_inviscid[j].density; 641 for (int k=0; k<3; k++) 642 Flux[k+1][j] = F_inviscid[j].momentum[k] - stress[k][j]; 643 Flux[4][j] = F_inviscid[j].E_total + Fe[j]; 644 } 645 646 for (int j=0; j<3; j++) { 647 for (int k=0; k<5; k++) { 648 dv[j][k][i] = -wdetJ * (dXdx[j][0] * Flux[k][0] + 649 dXdx[j][1] * Flux[k][1] + 650 dXdx[j][2] * Flux[k][2]); 651 } 652 } 653 654 const CeedScalar body_force[5] = {0, s.U.density *g[0], s.U.density *g[1], s.U.density *g[2], 0}; 655 for (int j=0; j<5; j++) 656 v[j][i] = wdetJ * (q_dot[j][i] - body_force[j]); 657 658 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 659 CeedScalar jacob_F_conv[3][5][5] = {0}; 660 computeFluxJacobian_NS(jacob_F_conv, s.U.density, s.Y.velocity, s.U.E_total, 661 gamma, g, x_i); 662 CeedScalar grad_U[5][3]; 663 for (int j=0; j<3; j++) { 664 grad_U[0][j] = grad_s[j].U.density; 665 for (int k=0; k<3; k++) grad_U[k+1][j] = grad_s[j].U.momentum[k]; 666 grad_U[4][j] = grad_s[j].U.E_total; 667 } 668 669 // strong_conv = dF/dq * dq/dx (Strong convection) 670 CeedScalar strong_conv[5] = {0}; 671 for (int j=0; j<3; j++) 672 for (int k=0; k<5; k++) 673 for (int l=0; l<5; l++) 674 strong_conv[k] += jacob_F_conv[j][k][l] * grad_U[l][j]; 675 676 // Strong residual 677 CeedScalar strong_res[5]; 678 for (int j=0; j<5; j++) 679 strong_res[j] = q_dot[j][i] + strong_conv[j] - body_force[j]; 680 681 // -- Stabilization method: none, SU, or SUPG 682 CeedScalar stab[5][3] = {{0.}}; 683 CeedScalar tau_strong_res[5] = {0.}, tau_strong_res_conservative[5] = {0}; 684 CeedScalar tau_strong_conv[5] = {0.}, tau_strong_conv_conservative[5] = {0}; 685 CeedScalar Tau_d[3] = {0.}; 686 switch (context->stabilization) { 687 case STAB_NONE: // Galerkin 688 break; 689 case STAB_SU: // SU 690 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 691 tau_strong_conv[0] = Tau_d[0] * strong_conv[0]; 692 tau_strong_conv[1] = Tau_d[1] * strong_conv[1]; 693 tau_strong_conv[2] = Tau_d[1] * strong_conv[2]; 694 tau_strong_conv[3] = Tau_d[1] * strong_conv[3]; 695 tau_strong_conv[4] = Tau_d[2] * strong_conv[4]; 696 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 697 tau_strong_conv, tau_strong_conv_conservative); 698 for (int j=0; j<3; j++) 699 for (int k=0; k<5; k++) 700 for (int l=0; l<5; l++) 701 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_conv_conservative[l]; 702 703 for (int j=0; j<5; j++) 704 for (int k=0; k<3; k++) 705 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 706 stab[j][1] * dXdx[k][1] + 707 stab[j][2] * dXdx[k][2]); 708 break; 709 case STAB_SUPG: // SUPG 710 Tau_diagPrim(Tau_d, dXdx, s.Y.velocity, cv, context, mu, dt, s.U.density); 711 tau_strong_res[0] = Tau_d[0] * strong_res[0]; 712 tau_strong_res[1] = Tau_d[1] * strong_res[1]; 713 tau_strong_res[2] = Tau_d[1] * strong_res[2]; 714 tau_strong_res[3] = Tau_d[1] * strong_res[3]; 715 tau_strong_res[4] = Tau_d[2] * strong_res[4]; 716 // Alternate route (useful later with primitive variable code) 717 // this function was verified against PHASTA for as IC that was as close as possible 718 // computeFluxJacobian_NSp(jacob_F_conv_p, rho, u, E, Rd, cv); 719 // it has also been verified to compute a correct through the following 720 // stab[k][j] += jacob_F_conv_p[j][k][l] * tau_strong_res[l] // flux Jacobian wrt primitive 721 // applied in the triple loop below 722 // However, it is more flops than using the existing Jacobian wrt q after q_{,Y} viz 723 PrimitiveToConservative_fwd(s.U.density, s.Y.velocity, s.U.E_total, Rd, cv, 724 tau_strong_res, tau_strong_res_conservative); 725 for (int j=0; j<3; j++) 726 for (int k=0; k<5; k++) 727 for (int l=0; l<5; l++) 728 stab[k][j] += jacob_F_conv[j][k][l] * tau_strong_res_conservative[l]; 729 730 for (int j=0; j<5; j++) 731 for (int k=0; k<3; k++) 732 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 733 stab[j][1] * dXdx[k][1] + 734 stab[j][2] * dXdx[k][2]); 735 break; 736 } 737 738 } // End Quadrature Point Loop 739 740 // Return 741 return 0; 742 } 743 // ***************************************************************************** 744 #endif // newtonian_h 745