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