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