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