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