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