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