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 18 #ifndef M_PI 19 #define M_PI 3.14159265358979323846 20 #endif 21 22 #ifndef setup_context_struct 23 #define setup_context_struct 24 typedef struct SetupContext_ *SetupContext; 25 struct SetupContext_ { 26 CeedScalar theta0; 27 CeedScalar thetaC; 28 CeedScalar P0; 29 CeedScalar N; 30 CeedScalar cv; 31 CeedScalar cp; 32 CeedScalar g; 33 CeedScalar rc; 34 CeedScalar lx; 35 CeedScalar ly; 36 CeedScalar lz; 37 CeedScalar center[3]; 38 CeedScalar dc_axis[3]; 39 CeedScalar wind[3]; 40 CeedScalar time; 41 int wind_type; // See WindType: 0=ROTATION, 1=TRANSLATION 42 int bubble_type; // See BubbleType: 0=SPHERE, 1=CYLINDER 43 int bubble_continuity_type; // See BubbleContinuityType: 0=SMOOTH, 1=BACK_SHARP 2=THICK 44 }; 45 #endif 46 47 #ifndef newtonian_context_struct 48 #define newtonian_context_struct 49 typedef enum { 50 STAB_NONE = 0, 51 STAB_SU = 1, // Streamline Upwind 52 STAB_SUPG = 2, // Streamline Upwind Petrov-Galerkin 53 } StabilizationType; 54 55 typedef struct NewtonianIdealGasContext_ *NewtonianIdealGasContext; 56 struct NewtonianIdealGasContext_ { 57 CeedScalar lambda; 58 CeedScalar mu; 59 CeedScalar k; 60 CeedScalar cv; 61 CeedScalar cp; 62 CeedScalar g; 63 CeedScalar c_tau; 64 StabilizationType stabilization; 65 }; 66 #endif 67 68 // ***************************************************************************** 69 // Helper function for computing flux Jacobian 70 // ***************************************************************************** 71 CEED_QFUNCTION_HELPER void computeFluxJacobian_NS(CeedScalar dF[3][5][5], 72 const CeedScalar rho, const CeedScalar u[3], const CeedScalar E, 73 const CeedScalar gamma, const CeedScalar g, CeedScalar z) { 74 CeedScalar u_sq = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; // Velocity square 75 for (CeedInt i=0; i<3; i++) { // Jacobian matrices for 3 directions 76 for (CeedInt j=0; j<3; j++) { // Rows of each Jacobian matrix 77 dF[i][j+1][0] = ((i==j) ? ((gamma-1.)*(u_sq/2. - g*z)) : 0.) - u[i]*u[j]; 78 for (CeedInt k=0; k<3; k++) { // Columns of each Jacobian matrix 79 dF[i][0][k+1] = ((i==k) ? 1. : 0.); 80 dF[i][j+1][k+1] = ((j==k) ? u[i] : 0.) + 81 ((i==k) ? u[j] : 0.) - 82 ((i==j) ? u[k] : 0.) * (gamma-1.); 83 dF[i][4][k+1] = ((i==k) ? (E*gamma/rho - (gamma-1.)*u_sq/2.) : 0.) - 84 (gamma-1.)*u[i]*u[k]; 85 } 86 dF[i][j+1][4] = ((i==j) ? (gamma-1.) : 0.); 87 } 88 dF[i][4][0] = u[i] * ((gamma-1.)*u_sq - E*gamma/rho); 89 dF[i][4][4] = u[i] * gamma; 90 } 91 } 92 93 // ***************************************************************************** 94 // Helper function for computing Tau elements (stabilization constant) 95 // Model from: 96 // Stabilized Methods for Compressible Flows, Hughes et al 2010 97 // 98 // Spatial criterion #2 - Tau is a 3x3 diagonal matrix 99 // Tau[i] = c_tau h[i] Xi(Pe) / rho(A[i]) (no sum) 100 // 101 // Where 102 // c_tau = stabilization constant (0.5 is reported as "optimal") 103 // h[i] = 2 length(dxdX[i]) 104 // Pe = Peclet number ( Pe = sqrt(u u) / dot(dXdx,u) diffusivity ) 105 // Xi(Pe) = coth Pe - 1. / Pe (1. at large local Peclet number ) 106 // rho(A[i]) = spectral radius of the convective flux Jacobian i, 107 // wave speed in direction i 108 // ***************************************************************************** 109 CEED_QFUNCTION_HELPER void Tau_spatial(CeedScalar Tau_x[3], 110 const CeedScalar dXdx[3][3], const CeedScalar u[3], 111 const CeedScalar sound_speed, const CeedScalar c_tau) { 112 for (int i=0; i<3; i++) { 113 // length of element in direction i 114 CeedScalar h = 2 / sqrt(dXdx[0][i]*dXdx[0][i] + dXdx[1][i]*dXdx[1][i] + 115 dXdx[2][i]*dXdx[2][i]); 116 // fastest wave in direction i 117 CeedScalar fastest_wave = fabs(u[i]) + sound_speed; 118 Tau_x[i] = c_tau * h / fastest_wave; 119 } 120 } 121 122 // ***************************************************************************** 123 // This QFunction sets a "still" initial condition for generic Newtonian IG problems 124 // ***************************************************************************** 125 CEED_QFUNCTION(ICsNewtonianIG)(void *ctx, CeedInt Q, 126 const CeedScalar *const *in, CeedScalar *const *out) { 127 // Inputs 128 const CeedScalar (*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 129 130 // Outputs 131 CeedScalar (*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 132 133 // Quadrature Point Loop 134 CeedPragmaSIMD 135 for (CeedInt i=0; i<Q; i++) { 136 CeedScalar q[5] = {0.}; 137 138 // Context 139 const SetupContext context = (SetupContext)ctx; 140 const CeedScalar theta0 = context->theta0; 141 const CeedScalar P0 = context->P0; 142 const CeedScalar N = context->N; 143 const CeedScalar cv = context->cv; 144 const CeedScalar cp = context->cp; 145 const CeedScalar g = context->g; 146 const CeedScalar Rd = cp - cv; 147 148 // Setup 149 // -- Coordinates 150 const CeedScalar z = X[2][i]; 151 152 // -- Exner pressure, hydrostatic balance 153 const CeedScalar Pi = 1. + g*g*(exp(-N*N*z/g) - 1.) / (cp*theta0*N*N); 154 155 // -- Density 156 const CeedScalar rho = P0 * pow(Pi, cv/Rd) / (Rd*theta0); 157 158 // Initial Conditions 159 q[0] = rho; 160 q[1] = 0.0; 161 q[2] = 0.0; 162 q[3] = 0.0; 163 q[4] = rho * (cv*theta0*Pi + g*z); 164 165 for (CeedInt j=0; j<5; j++) 166 q0[j][i] = q[j]; 167 } // End of Quadrature Point Loop 168 return 0; 169 } 170 171 // ***************************************************************************** 172 // This QFunction implements the following formulation of Navier-Stokes with 173 // explicit time stepping method 174 // 175 // This is 3D compressible Navier-Stokes in conservation form with state 176 // variables of density, momentum density, and total energy density. 177 // 178 // State Variables: q = ( rho, U1, U2, U3, E ) 179 // rho - Mass Density 180 // Ui - Momentum Density, Ui = rho ui 181 // E - Total Energy Density, E = rho (cv T + (u u)/2 + g z) 182 // 183 // Navier-Stokes Equations: 184 // drho/dt + div( U ) = 0 185 // dU/dt + div( rho (u x u) + P I3 ) + rho g khat = div( Fu ) 186 // dE/dt + div( (E + P) u ) = div( Fe ) 187 // 188 // Viscous Stress: 189 // Fu = mu (grad( u ) + grad( u )^T + lambda div ( u ) I3) 190 // 191 // Thermal Stress: 192 // Fe = u Fu + k grad( T ) 193 // 194 // Equation of State: 195 // P = (gamma - 1) (E - rho (u u) / 2 - rho g z) 196 // 197 // Stabilization: 198 // Tau = diag(TauC, TauM, TauM, TauM, TauE) 199 // f1 = rho sqrt(ui uj gij) 200 // gij = dXi/dX * dXi/dX 201 // TauC = Cc f1 / (8 gii) 202 // TauM = min( 1 , 1 / f1 ) 203 // TauE = TauM / (Ce cv) 204 // 205 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 206 // 207 // Constants: 208 // lambda = - 2 / 3, From Stokes hypothesis 209 // mu , Dynamic viscosity 210 // k , Thermal conductivity 211 // cv , Specific heat, constant volume 212 // cp , Specific heat, constant pressure 213 // g , Gravity 214 // gamma = cp / cv, Specific heat ratio 215 // 216 // We require the product of the inverse of the Jacobian (dXdx_j,k) and 217 // its transpose (dXdx_k,j) to properly compute integrals of the form: 218 // int( gradv gradu ) 219 // 220 // ***************************************************************************** 221 CEED_QFUNCTION(Newtonian)(void *ctx, CeedInt Q, 222 const CeedScalar *const *in, CeedScalar *const *out) { 223 // *INDENT-OFF* 224 // Inputs 225 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 226 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 227 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 228 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 229 // Outputs 230 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 231 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 232 // *INDENT-ON* 233 234 // Context 235 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 236 const CeedScalar lambda = context->lambda; 237 const CeedScalar mu = context->mu; 238 const CeedScalar k = context->k; 239 const CeedScalar cv = context->cv; 240 const CeedScalar cp = context->cp; 241 const CeedScalar g = context->g; 242 const CeedScalar c_tau = context->c_tau; 243 const CeedScalar gamma = cp / cv; 244 245 CeedPragmaSIMD 246 // Quadrature Point Loop 247 for (CeedInt i=0; i<Q; i++) { 248 // *INDENT-OFF* 249 // Setup 250 // -- Interp in 251 const CeedScalar rho = q[0][i]; 252 const CeedScalar u[3] = {q[1][i] / rho, 253 q[2][i] / rho, 254 q[3][i] / rho 255 }; 256 const CeedScalar E = q[4][i]; 257 // -- Grad in 258 const CeedScalar drho[3] = {dq[0][0][i], 259 dq[1][0][i], 260 dq[2][0][i] 261 }; 262 const CeedScalar dU[3][3] = {{dq[0][1][i], 263 dq[1][1][i], 264 dq[2][1][i]}, 265 {dq[0][2][i], 266 dq[1][2][i], 267 dq[2][2][i]}, 268 {dq[0][3][i], 269 dq[1][3][i], 270 dq[2][3][i]} 271 }; 272 const CeedScalar dE[3] = {dq[0][4][i], 273 dq[1][4][i], 274 dq[2][4][i] 275 }; 276 // -- Interp-to-Interp q_data 277 const CeedScalar wdetJ = q_data[0][i]; 278 // -- Interp-to-Grad q_data 279 // ---- Inverse of change of coordinate matrix: X_i,j 280 // *INDENT-OFF* 281 const CeedScalar dXdx[3][3] = {{q_data[1][i], 282 q_data[2][i], 283 q_data[3][i]}, 284 {q_data[4][i], 285 q_data[5][i], 286 q_data[6][i]}, 287 {q_data[7][i], 288 q_data[8][i], 289 q_data[9][i]} 290 }; 291 // *INDENT-ON* 292 // -- Grad-to-Grad q_data 293 // dU/dx 294 CeedScalar du[3][3] = {{0}}; 295 CeedScalar drhodx[3] = {0}; 296 CeedScalar dEdx[3] = {0}; 297 CeedScalar dUdx[3][3] = {{0}}; 298 CeedScalar dXdxdXdxT[3][3] = {{0}}; 299 for (int j=0; j<3; j++) { 300 for (int k=0; k<3; k++) { 301 du[j][k] = (dU[j][k] - drho[k]*u[j]) / rho; 302 drhodx[j] += drho[k] * dXdx[k][j]; 303 dEdx[j] += dE[k] * dXdx[k][j]; 304 for (int l=0; l<3; l++) { 305 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 306 dXdxdXdxT[j][k] += dXdx[j][l]*dXdx[k][l]; //dXdx_j,k * dXdx_k,j 307 } 308 } 309 } 310 CeedScalar dudx[3][3] = {{0}}; 311 for (int j=0; j<3; j++) 312 for (int k=0; k<3; k++) 313 for (int l=0; l<3; l++) 314 dudx[j][k] += du[j][l] * dXdx[l][k]; 315 // -- grad_T 316 const CeedScalar grad_T[3] = {(dEdx[0]/rho - E*drhodx[0]/(rho*rho) - /* *NOPAD* */ 317 (u[0]*dudx[0][0] + u[1]*dudx[1][0] + u[2]*dudx[2][0]))/cv, 318 (dEdx[1]/rho - E*drhodx[1]/(rho*rho) - /* *NOPAD* */ 319 (u[0]*dudx[0][1] + u[1]*dudx[1][1] + u[2]*dudx[2][1]))/cv, 320 (dEdx[2]/rho - E*drhodx[2]/(rho*rho) - /* *NOPAD* */ 321 (u[0]*dudx[0][2] + u[1]*dudx[1][2] + u[2]*dudx[2][2]) - g)/cv 322 }; 323 324 // -- Fuvisc 325 // ---- Symmetric 3x3 matrix 326 const CeedScalar Fu[6] = {mu*(dudx[0][0] * (2 + lambda) + /* *NOPAD* */ 327 lambda * (dudx[1][1] + dudx[2][2])), 328 mu*(dudx[0][1] + dudx[1][0]), /* *NOPAD* */ 329 mu*(dudx[0][2] + dudx[2][0]), /* *NOPAD* */ 330 mu*(dudx[1][1] * (2 + lambda) + /* *NOPAD* */ 331 lambda * (dudx[0][0] + dudx[2][2])), 332 mu*(dudx[1][2] + dudx[2][1]), /* *NOPAD* */ 333 mu*(dudx[2][2] * (2 + lambda) + /* *NOPAD* */ 334 lambda * (dudx[0][0] + dudx[1][1])) 335 }; 336 // -- Fevisc 337 const CeedScalar Fe[3] = {u[0]*Fu[0] + u[1]*Fu[1] + u[2]*Fu[2] + /* *NOPAD* */ 338 k*grad_T[0], /* *NOPAD* */ 339 u[0]*Fu[1] + u[1]*Fu[3] + u[2]*Fu[4] + /* *NOPAD* */ 340 k*grad_T[1], /* *NOPAD* */ 341 u[0]*Fu[2] + u[1]*Fu[4] + u[2]*Fu[5] + /* *NOPAD* */ 342 k*grad_T[2] /* *NOPAD* */ 343 }; 344 // Pressure 345 const CeedScalar 346 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 347 E_potential = rho*g*x[2][i], 348 E_internal = E - E_kinetic - E_potential, 349 P = E_internal * (gamma - 1.); // P = pressure 350 351 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 352 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 353 computeFluxJacobian_NS(jacob_F_conv, rho, u, E, gamma, g, x[2][i]); 354 355 // jacob_F_conv_T = jacob_F_conv^T 356 CeedScalar jacob_F_conv_T[3][5][5]; 357 for (int j=0; j<3; j++) 358 for (int k=0; k<5; k++) 359 for (int l=0; l<5; l++) 360 jacob_F_conv_T[j][k][l] = jacob_F_conv[j][l][k]; 361 362 // dqdx collects drhodx, dUdx and dEdx in one vector 363 CeedScalar dqdx[5][3]; 364 for (int j=0; j<3; j++) { 365 dqdx[0][j] = drhodx[j]; 366 dqdx[4][j] = dEdx[j]; 367 for (int k=0; k<3; k++) 368 dqdx[k+1][j] = dUdx[k][j]; 369 } 370 371 // strong_conv = dF/dq * dq/dx (Strong convection) 372 CeedScalar strong_conv[5] = {0}; 373 for (int j=0; j<3; j++) 374 for (int k=0; k<5; k++) 375 for (int l=0; l<5; l++) 376 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 377 378 // Body force 379 const CeedScalar body_force[5] = {0, 0, 0, -rho*g, 0}; 380 381 // The Physics 382 // Zero dv so all future terms can safely sum into it 383 for (int j=0; j<5; j++) 384 for (int k=0; k<3; k++) 385 dv[k][j][i] = 0; 386 387 // -- Density 388 // ---- u rho 389 for (int j=0; j<3; j++) 390 dv[j][0][i] += wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 391 rho*u[2]*dXdx[j][2]); 392 // -- Momentum 393 // ---- rho (u x u) + P I3 394 for (int j=0; j<3; j++) 395 for (int k=0; k<3; k++) 396 dv[k][j+1][i] += wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 397 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 398 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 399 // ---- Fuvisc 400 const CeedInt Fuviscidx[3][3] = {{0, 1, 2}, {1, 3, 4}, {2, 4, 5}}; // symmetric matrix indices 401 for (int j=0; j<3; j++) 402 for (int k=0; k<3; k++) 403 dv[k][j+1][i] -= wdetJ*(Fu[Fuviscidx[j][0]]*dXdx[k][0] + 404 Fu[Fuviscidx[j][1]]*dXdx[k][1] + 405 Fu[Fuviscidx[j][2]]*dXdx[k][2]); 406 // -- Total Energy Density 407 // ---- (E + P) u 408 for (int j=0; j<3; j++) 409 dv[j][4][i] += wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 410 u[2]*dXdx[j][2]); 411 // ---- Fevisc 412 for (int j=0; j<3; j++) 413 dv[j][4][i] -= wdetJ * (Fe[0]*dXdx[j][0] + Fe[1]*dXdx[j][1] + 414 Fe[2]*dXdx[j][2]); 415 // Body Force 416 for (int j=0; j<5; j++) 417 v[j][i] = wdetJ * body_force[j]; 418 419 // Stabilization 420 // -- Tau elements 421 const CeedScalar sound_speed = sqrt(gamma * P / rho); 422 CeedScalar Tau_x[3] = {0.}; 423 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 424 425 // -- Stabilization method: none or SU 426 CeedScalar stab[5][3]; 427 switch (context->stabilization) { 428 case STAB_NONE: // Galerkin 429 break; 430 case STAB_SU: // SU 431 for (int j=0; j<3; j++) 432 for (int k=0; k<5; k++) 433 for (int l=0; l<5; l++) 434 stab[k][j] = jacob_F_conv_T[j][k][l] * Tau_x[j] * strong_conv[l]; 435 436 for (int j=0; j<5; j++) 437 for (int k=0; k<3; k++) 438 dv[k][j][i] -= wdetJ*(stab[j][0] * dXdx[k][0] + 439 stab[j][1] * dXdx[k][1] + 440 stab[j][2] * dXdx[k][2]); 441 break; 442 case STAB_SUPG: // SUPG is not implemented for explicit scheme 443 break; 444 } 445 446 } // End Quadrature Point Loop 447 448 // Return 449 return 0; 450 } 451 452 // ***************************************************************************** 453 // This QFunction implements the Navier-Stokes equations (mentioned above) with 454 // implicit time stepping method 455 // 456 // SU = Galerkin + grad(v) . ( Ai^T * Tau * (Aj q,j) ) 457 // SUPG = Galerkin + grad(v) . ( Ai^T * Tau * (q_dot + Aj q,j - body force) ) 458 // (diffussive terms will be added later) 459 // 460 // ***************************************************************************** 461 CEED_QFUNCTION(IFunction_Newtonian)(void *ctx, CeedInt Q, 462 const CeedScalar *const *in, 463 CeedScalar *const *out) { 464 // *INDENT-OFF* 465 // Inputs 466 const CeedScalar (*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0], 467 (*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1], 468 (*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], 469 (*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3], 470 (*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 471 // Outputs 472 CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], 473 (*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 474 // *INDENT-ON* 475 // Context 476 NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 477 const CeedScalar lambda = context->lambda; 478 const CeedScalar mu = context->mu; 479 const CeedScalar k = context->k; 480 const CeedScalar cv = context->cv; 481 const CeedScalar cp = context->cp; 482 const CeedScalar g = context->g; 483 const CeedScalar c_tau = context->c_tau; 484 const CeedScalar gamma = cp / cv; 485 486 CeedPragmaSIMD 487 // Quadrature Point Loop 488 for (CeedInt i=0; i<Q; i++) { 489 // Setup 490 // -- Interp in 491 const CeedScalar rho = q[0][i]; 492 const CeedScalar u[3] = {q[1][i] / rho, 493 q[2][i] / rho, 494 q[3][i] / rho 495 }; 496 const CeedScalar E = q[4][i]; 497 // -- Grad in 498 const CeedScalar drho[3] = {dq[0][0][i], 499 dq[1][0][i], 500 dq[2][0][i] 501 }; 502 // *INDENT-OFF* 503 const CeedScalar dU[3][3] = {{dq[0][1][i], 504 dq[1][1][i], 505 dq[2][1][i]}, 506 {dq[0][2][i], 507 dq[1][2][i], 508 dq[2][2][i]}, 509 {dq[0][3][i], 510 dq[1][3][i], 511 dq[2][3][i]} 512 }; 513 // *INDENT-ON* 514 const CeedScalar dE[3] = {dq[0][4][i], 515 dq[1][4][i], 516 dq[2][4][i] 517 }; 518 // -- Interp-to-Interp q_data 519 const CeedScalar wdetJ = q_data[0][i]; 520 // -- Interp-to-Grad q_data 521 // ---- Inverse of change of coordinate matrix: X_i,j 522 // *INDENT-OFF* 523 const CeedScalar dXdx[3][3] = {{q_data[1][i], 524 q_data[2][i], 525 q_data[3][i]}, 526 {q_data[4][i], 527 q_data[5][i], 528 q_data[6][i]}, 529 {q_data[7][i], 530 q_data[8][i], 531 q_data[9][i]} 532 }; 533 // *INDENT-ON* 534 // -- Grad-to-Grad q_data 535 // dU/dx 536 CeedScalar du[3][3] = {{0}}; 537 CeedScalar drhodx[3] = {0}; 538 CeedScalar dEdx[3] = {0}; 539 CeedScalar dUdx[3][3] = {{0}}; 540 CeedScalar dXdxdXdxT[3][3] = {{0}}; 541 for (int j=0; j<3; j++) { 542 for (int k=0; k<3; k++) { 543 du[j][k] = (dU[j][k] - drho[k]*u[j]) / rho; 544 drhodx[j] += drho[k] * dXdx[k][j]; 545 dEdx[j] += dE[k] * dXdx[k][j]; 546 for (int l=0; l<3; l++) { 547 dUdx[j][k] += dU[j][l] * dXdx[l][k]; 548 dXdxdXdxT[j][k] += dXdx[j][l]*dXdx[k][l]; //dXdx_j,k * dXdx_k,j 549 } 550 } 551 } 552 CeedScalar dudx[3][3] = {{0}}; 553 for (int j=0; j<3; j++) 554 for (int k=0; k<3; k++) 555 for (int l=0; l<3; l++) 556 dudx[j][k] += du[j][l] * dXdx[l][k]; 557 // -- grad_T 558 const CeedScalar grad_T[3] = {(dEdx[0]/rho - E*drhodx[0]/(rho*rho) - /* *NOPAD* */ 559 (u[0]*dudx[0][0] + u[1]*dudx[1][0] + u[2]*dudx[2][0]))/cv, 560 (dEdx[1]/rho - E*drhodx[1]/(rho*rho) - /* *NOPAD* */ 561 (u[0]*dudx[0][1] + u[1]*dudx[1][1] + u[2]*dudx[2][1]))/cv, 562 (dEdx[2]/rho - E*drhodx[2]/(rho*rho) - /* *NOPAD* */ 563 (u[0]*dudx[0][2] + u[1]*dudx[1][2] + u[2]*dudx[2][2]) - g)/cv 564 }; 565 // -- Fuvisc 566 // ---- Symmetric 3x3 matrix 567 const CeedScalar Fu[6] = {mu*(dudx[0][0] * (2 + lambda) + /* *NOPAD* */ 568 lambda * (dudx[1][1] + dudx[2][2])), 569 mu*(dudx[0][1] + dudx[1][0]), /* *NOPAD* */ 570 mu*(dudx[0][2] + dudx[2][0]), /* *NOPAD* */ 571 mu*(dudx[1][1] * (2 + lambda) + /* *NOPAD* */ 572 lambda * (dudx[0][0] + dudx[2][2])), 573 mu*(dudx[1][2] + dudx[2][1]), /* *NOPAD* */ 574 mu*(dudx[2][2] * (2 + lambda) + /* *NOPAD* */ 575 lambda * (dudx[0][0] + dudx[1][1])) 576 }; 577 // -- Fevisc 578 const CeedScalar Fe[3] = {u[0]*Fu[0] + u[1]*Fu[1] + u[2]*Fu[2] + /* *NOPAD* */ 579 k*grad_T[0], /* *NOPAD* */ 580 u[0]*Fu[1] + u[1]*Fu[3] + u[2]*Fu[4] + /* *NOPAD* */ 581 k*grad_T[1], /* *NOPAD* */ 582 u[0]*Fu[2] + u[1]*Fu[4] + u[2]*Fu[5] + /* *NOPAD* */ 583 k*grad_T[2] /* *NOPAD* */ 584 }; 585 // Pressure 586 const CeedScalar 587 E_kinetic = 0.5 * rho * (u[0]*u[0] + u[1]*u[1] + u[2]*u[2]), 588 E_potential = rho*g*x[2][i], 589 E_internal = E - E_kinetic - E_potential, 590 P = E_internal * (gamma - 1.); // P = pressure 591 592 // jacob_F_conv[3][5][5] = dF(convective)/dq at each direction 593 CeedScalar jacob_F_conv[3][5][5] = {{{0.}}}; 594 computeFluxJacobian_NS(jacob_F_conv, rho, u, E, gamma, g, x[2][i]); 595 596 // jacob_F_conv_T = jacob_F_conv^T 597 CeedScalar jacob_F_conv_T[3][5][5]; 598 for (int j=0; j<3; j++) 599 for (int k=0; k<5; k++) 600 for (int l=0; l<5; l++) 601 jacob_F_conv_T[j][k][l] = jacob_F_conv[j][l][k]; 602 // dqdx collects drhodx, dUdx and dEdx in one vector 603 CeedScalar dqdx[5][3]; 604 for (int j=0; j<3; j++) { 605 dqdx[0][j] = drhodx[j]; 606 dqdx[4][j] = dEdx[j]; 607 for (int k=0; k<3; k++) 608 dqdx[k+1][j] = dUdx[k][j]; 609 } 610 // strong_conv = dF/dq * dq/dx (Strong convection) 611 CeedScalar strong_conv[5] = {0}; 612 for (int j=0; j<3; j++) 613 for (int k=0; k<5; k++) 614 for (int l=0; l<5; l++) 615 strong_conv[k] += jacob_F_conv[j][k][l] * dqdx[l][j]; 616 617 // Body force 618 const CeedScalar body_force[5] = {0, 0, 0, -rho*g, 0}; 619 620 // Strong residual 621 CeedScalar strong_res[5]; 622 for (int j=0; j<5; j++) 623 strong_res[j] = q_dot[j][i] + strong_conv[j] - body_force[j]; 624 625 // The Physics 626 //-----mass matrix 627 for (int j=0; j<5; j++) 628 v[j][i] = wdetJ*q_dot[j][i]; 629 630 // Zero dv so all future terms can safely sum into it 631 for (int j=0; j<5; j++) 632 for (int k=0; k<3; k++) 633 dv[k][j][i] = 0; 634 635 // -- Density 636 // ---- u rho 637 for (int j=0; j<3; j++) 638 dv[j][0][i] -= wdetJ*(rho*u[0]*dXdx[j][0] + rho*u[1]*dXdx[j][1] + 639 rho*u[2]*dXdx[j][2]); 640 // -- Momentum 641 // ---- rho (u x u) + P I3 642 for (int j=0; j<3; j++) 643 for (int k=0; k<3; k++) 644 dv[k][j+1][i] -= wdetJ*((rho*u[j]*u[0] + (j==0?P:0))*dXdx[k][0] + 645 (rho*u[j]*u[1] + (j==1?P:0))*dXdx[k][1] + 646 (rho*u[j]*u[2] + (j==2?P:0))*dXdx[k][2]); 647 // ---- Fuvisc 648 const CeedInt Fuviscidx[3][3] = {{0, 1, 2}, {1, 3, 4}, {2, 4, 5}}; // symmetric matrix indices 649 for (int j=0; j<3; j++) 650 for (int k=0; k<3; k++) 651 dv[k][j+1][i] += wdetJ*(Fu[Fuviscidx[j][0]]*dXdx[k][0] + 652 Fu[Fuviscidx[j][1]]*dXdx[k][1] + 653 Fu[Fuviscidx[j][2]]*dXdx[k][2]); 654 // -- Total Energy Density 655 // ---- (E + P) u 656 for (int j=0; j<3; j++) 657 dv[j][4][i] -= wdetJ * (E + P) * (u[0]*dXdx[j][0] + u[1]*dXdx[j][1] + 658 u[2]*dXdx[j][2]); 659 // ---- Fevisc 660 for (int j=0; j<3; j++) 661 dv[j][4][i] += wdetJ * (Fe[0]*dXdx[j][0] + Fe[1]*dXdx[j][1] + 662 Fe[2]*dXdx[j][2]); 663 // Body Force 664 for (int j=0; j<5; j++) 665 v[j][i] -= wdetJ*body_force[j]; 666 667 // Stabilization 668 // -- Tau elements 669 const CeedScalar sound_speed = sqrt(gamma * P / rho); 670 CeedScalar Tau_x[3] = {0.}; 671 Tau_spatial(Tau_x, dXdx, u, sound_speed, c_tau); 672 673 // -- Stabilization method: none, SU, or SUPG 674 CeedScalar stab[5][3]; 675 switch (context->stabilization) { 676 case STAB_NONE: // Galerkin 677 break; 678 case STAB_SU: // SU 679 for (int j=0; j<3; j++) 680 for (int k=0; k<5; k++) 681 for (int l=0; l<5; l++) 682 stab[k][j] = jacob_F_conv_T[j][k][l] * Tau_x[j] * strong_conv[l]; 683 684 for (int j=0; j<5; j++) 685 for (int k=0; k<3; k++) 686 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 687 stab[j][1] * dXdx[k][1] + 688 stab[j][2] * dXdx[k][2]); 689 break; 690 case STAB_SUPG: // SUPG 691 for (int j=0; j<3; j++) 692 for (int k=0; k<5; k++) 693 for (int l=0; l<5; l++) 694 stab[k][j] = jacob_F_conv_T[j][k][l] * Tau_x[j] * strong_res[l]; 695 696 for (int j=0; j<5; j++) 697 for (int k=0; k<3; k++) 698 dv[k][j][i] += wdetJ*(stab[j][0] * dXdx[k][0] + 699 stab[j][1] * dXdx[k][1] + 700 stab[j][2] * dXdx[k][2]); 701 break; 702 } 703 704 } // End Quadrature Point Loop 705 706 // Return 707 return 0; 708 } 709 // ***************************************************************************** 710 #endif // newtonian_h 711