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