1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Advection initial condition and operator for Navier-Stokes example using PETSc 10 #include <ceed.h> 11 #include <math.h> 12 13 #include "advection_types.h" 14 #include "newtonian_state.h" 15 #include "newtonian_types.h" 16 #include "stabilization_types.h" 17 #include "utils.h" 18 19 // ***************************************************************************** 20 // This QFunction sets the initial conditions and the boundary conditions 21 // for two test cases: ROTATION and TRANSLATION 22 // 23 // -- ROTATION (default) 24 // Initial Conditions: 25 // Mass Density: 26 // Constant mass density of 1.0 27 // Momentum Density: 28 // Rotational field in x,y 29 // Energy Density: 30 // Maximum of 1. x0 decreasing linearly to 0. as radial distance 31 // increases to (1.-r/rc), then 0. everywhere else 32 // 33 // Boundary Conditions: 34 // Mass Density: 35 // 0.0 flux 36 // Momentum Density: 37 // 0.0 38 // Energy Density: 39 // 0.0 flux 40 // 41 // -- TRANSLATION 42 // Initial Conditions: 43 // Mass Density: 44 // Constant mass density of 1.0 45 // Momentum Density: 46 // Constant rectilinear field in x,y 47 // Energy Density: 48 // Maximum of 1. x0 decreasing linearly to 0. as radial distance 49 // increases to (1.-r/rc), then 0. everywhere else 50 // 51 // Boundary Conditions: 52 // Mass Density: 53 // 0.0 flux 54 // Momentum Density: 55 // 0.0 56 // Energy Density: 57 // Inflow BCs: 58 // E = E_wind 59 // Outflow BCs: 60 // E = E(boundary) 61 // Both In/Outflow BCs for E are applied weakly in the 62 // QFunction "Advection2d_Sur" 63 // 64 // ***************************************************************************** 65 66 // ***************************************************************************** 67 // This helper function provides the exact, time-dependent solution and IC formulation for 2D advection 68 // ***************************************************************************** 69 CEED_QFUNCTION_HELPER CeedInt Exact_AdvectionGeneric(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) { 70 const SetupContextAdv context = (SetupContextAdv)ctx; 71 const CeedScalar rc = context->rc; 72 const CeedScalar lx = context->lx; 73 const CeedScalar ly = context->ly; 74 const CeedScalar lz = dim == 2 ? 0. : context->lz; 75 const CeedScalar *wind = context->wind; 76 77 const CeedScalar center[3] = {0.5 * lx, 0.5 * ly, 0.5 * lz}; 78 const CeedScalar theta = dim == 2 ? M_PI / 3 : M_PI; 79 const CeedScalar x0[3] = {center[0] + .25 * lx * cos(theta + time), center[1] + .25 * ly * sin(theta + time), 0.5 * lz}; 80 81 const CeedScalar x = X[0], y = X[1], z = dim == 2 ? 0. : X[2]; 82 83 CeedScalar r = 0.; 84 switch (context->initial_condition_type) { 85 case ADVECTIONIC_BUBBLE_SPHERE: 86 case ADVECTIONIC_BUBBLE_CYLINDER: 87 r = sqrt(Square(x - x0[0]) + Square(y - x0[1]) + Square(z - x0[2])); 88 break; 89 case ADVECTIONIC_COSINE_HILL: 90 r = sqrt(Square(x - center[0]) + Square(y - center[1])); 91 break; 92 case ADVECTIONIC_SKEW: 93 break; 94 } 95 96 switch (context->wind_type) { 97 case WIND_ROTATION: 98 q[0] = 1.; 99 q[1] = -(y - center[1]); 100 q[2] = (x - center[0]); 101 q[3] = 0; 102 break; 103 case WIND_TRANSLATION: 104 q[0] = 1.; 105 q[1] = wind[0]; 106 q[2] = wind[1]; 107 q[3] = dim == 2 ? 0. : wind[2]; 108 break; 109 default: 110 return 1; 111 } 112 113 switch (context->initial_condition_type) { 114 case ADVECTIONIC_BUBBLE_SPHERE: 115 case ADVECTIONIC_BUBBLE_CYLINDER: 116 switch (context->bubble_continuity_type) { 117 // original continuous, smooth shape 118 case BUBBLE_CONTINUITY_SMOOTH: 119 q[4] = r <= rc ? (1. - r / rc) : 0.; 120 break; 121 // discontinuous, sharp back half shape 122 case BUBBLE_CONTINUITY_BACK_SHARP: 123 q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) : 0.; 124 break; 125 // attempt to define a finite thickness that will get resolved under grid refinement 126 case BUBBLE_CONTINUITY_THICK: 127 q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) * fmin(1.0, (center[1] - y) / 1.25) : 0.; 128 break; 129 case BUBBLE_CONTINUITY_COSINE: 130 q[4] = r <= rc ? .5 + .5 * cos(r * M_PI / rc) : 0; 131 break; 132 } 133 break; 134 case ADVECTIONIC_COSINE_HILL: { 135 CeedScalar half_width = context->lx / 2; 136 q[4] = r > half_width ? 0. : cos(2 * M_PI * r / half_width + M_PI) + 1.; 137 } break; 138 case ADVECTIONIC_SKEW: { 139 CeedScalar skewed_barrier[3] = {wind[0], wind[1], 0}; 140 CeedScalar inflow_to_point[3] = {x - context->lx / 2, y, 0}; 141 CeedScalar cross_product[3] = {0}; 142 const CeedScalar boundary_threshold = 20 * CEED_EPSILON; 143 Cross3(skewed_barrier, inflow_to_point, cross_product); 144 145 q[4] = cross_product[2] > boundary_threshold ? 0 : 1; 146 if ((x < boundary_threshold && wind[0] < boundary_threshold) || // outflow at -x boundary 147 (y < boundary_threshold && wind[1] < boundary_threshold) || // outflow at -y boundary 148 (x > context->lx - boundary_threshold && wind[0] > boundary_threshold) || // outflow at +x boundary 149 (y > context->ly - boundary_threshold && wind[1] > boundary_threshold) // outflow at +y boundary 150 ) { 151 q[4] = 0; 152 } 153 } break; 154 } 155 return 0; 156 } 157 158 // ***************************************************************************** 159 // This QFunction sets the initial conditions for 3D advection 160 // ***************************************************************************** 161 CEED_QFUNCTION(ICsAdvection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 162 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 163 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 164 165 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 166 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 167 CeedScalar q[5] = {0.}; 168 169 Exact_AdvectionGeneric(3, 0., x, 5, q, ctx); 170 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 171 } 172 return 0; 173 } 174 175 // ***************************************************************************** 176 // This QFunction sets the initial conditions for 2D advection 177 // ***************************************************************************** 178 CEED_QFUNCTION(ICsAdvection2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 179 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 180 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 181 const SetupContextAdv context = (SetupContextAdv)ctx; 182 183 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 184 const CeedScalar x[] = {X[0][i], X[1][i]}; 185 CeedScalar q[5] = {0.}; 186 187 Exact_AdvectionGeneric(2, context->time, x, 5, q, ctx); 188 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 189 } 190 return 0; 191 } 192 193 CEED_QFUNCTION_HELPER void QdataUnpack_ND(CeedInt N, CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar *dXdx) { 194 // Cannot directly use QdataUnpack* helper functions due to SYCL online compiler incompatabilities 195 switch (N) { 196 case 2: 197 StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 198 StoredValuesUnpack(Q, i, 1, 4, q_data, dXdx); 199 break; 200 case 3: 201 StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 202 StoredValuesUnpack(Q, i, 1, 9, q_data, dXdx); 203 break; 204 } 205 } 206 207 CEED_QFUNCTION_HELPER int QdataBoundaryUnpack_ND(CeedInt N, CeedInt Q, CeedInt i, const CeedScalar *q_data, CeedScalar *wdetJ, CeedScalar *dXdx, 208 CeedScalar *normal) { 209 // Cannot directly use QdataBoundaryUnpack* helper functions due to SYCL online compiler incompatabilities 210 switch (N) { 211 case 2: 212 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 213 if (normal) StoredValuesUnpack(Q, i, 1, 2, q_data, normal); 214 break; 215 case 3: 216 if (wdetJ) StoredValuesUnpack(Q, i, 0, 1, q_data, wdetJ); 217 if (normal) StoredValuesUnpack(Q, i, 1, 3, q_data, normal); 218 if (dXdx) StoredValuesUnpack(Q, i, 4, 6, q_data, (CeedScalar *)dXdx); 219 break; 220 } 221 return CEED_ERROR_SUCCESS; 222 } 223 224 CEED_QFUNCTION_HELPER void StatePhysicalGradientFromReference_ND(CeedInt N, CeedInt Q, CeedInt i, NewtonianIdealGasContext gas, State s, 225 StateVariable state_var, const CeedScalar *grad_q, const CeedScalar *dXdx, 226 State *grad_s) { 227 switch (N) { 228 case 2: { 229 for (CeedInt k = 0; k < 2; k++) { 230 CeedScalar dqi[5]; 231 for (CeedInt j = 0; j < 5; j++) { 232 dqi[j] = grad_q[(Q * 5) * 0 + Q * j + i] * dXdx[0 * N + k] + grad_q[(Q * 5) * 1 + Q * j + i] * dXdx[1 * N + k]; 233 } 234 grad_s[k] = StateFromQ_fwd(gas, s, dqi, state_var); 235 } 236 CeedScalar U[5] = {0.}; 237 grad_s[2] = StateFromU(gas, U); 238 } break; 239 case 3: 240 // Cannot directly use StatePhysicalGradientFromReference helper functions due to SYCL online compiler incompatabilities 241 for (CeedInt k = 0; k < 3; k++) { 242 CeedScalar dqi[5]; 243 for (CeedInt j = 0; j < 5; j++) { 244 dqi[j] = grad_q[(Q * 5) * 0 + Q * j + i] * dXdx[0 * N + k] + grad_q[(Q * 5) * 1 + Q * j + i] * dXdx[1 * N + k] + 245 grad_q[(Q * 5) * 2 + Q * j + i] * dXdx[2 * N + k]; 246 } 247 grad_s[k] = StateFromQ_fwd(gas, s, dqi, state_var); 248 } 249 break; 250 } 251 } 252 253 // @brief Calculate the stabilization constant \tau 254 CEED_QFUNCTION_HELPER CeedScalar Tau(AdvectionContext context, const State s, const CeedScalar *dXdx, CeedInt dim) { 255 switch (context->stabilization_tau) { 256 case STAB_TAU_CTAU: { 257 CeedScalar uX[3] = {0.}; 258 259 MatVecNM(dXdx, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, uX); 260 return context->CtauS / sqrt(DotN(uX, uX, dim)); 261 } break; 262 case STAB_TAU_ADVDIFF_SHAKIB: { 263 CeedScalar gijd_mat[9] = {0.}, gij_uj[3] = {0.}; 264 265 MatMatN(dXdx, dXdx, dim, CEED_TRANSPOSE, CEED_NOTRANSPOSE, gijd_mat); 266 MatVecNM(gijd_mat, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, gij_uj); 267 return 1 / sqrt(Square(2 * context->Ctau_t / context->dt) + DotN(s.Y.velocity, gij_uj, dim) * context->Ctau_a); 268 } break; 269 default: 270 return 0.; 271 } 272 } 273 274 // ***************************************************************************** 275 // This QFunction implements Advection for implicit time stepping method 276 // ***************************************************************************** 277 CEED_QFUNCTION_HELPER void IFunction_AdvectionGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, CeedInt dim) { 278 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 279 const CeedScalar(*grad_q) = in[1]; 280 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 281 const CeedScalar(*q_data) = in[3]; 282 283 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 284 CeedScalar(*grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 285 CeedScalar *jac_data = out[2]; 286 287 AdvectionContext context = (AdvectionContext)ctx; 288 const CeedScalar zeros[14] = {0.}; 289 NewtonianIdealGasContext gas; 290 struct NewtonianIdealGasContext_ gas_struct = {0}; 291 gas = &gas_struct; 292 293 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 294 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 295 const State s = StateFromU(gas, qi); 296 297 CeedScalar wdetJ, dXdx[9]; 298 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 299 State grad_s[3]; 300 StatePhysicalGradientFromReference_ND(dim, Q, i, gas, s, STATEVAR_CONSERVATIVE, grad_q, dXdx, grad_s); 301 302 const CeedScalar Grad_E[3] = {grad_s[0].U.E_total, grad_s[1].U.E_total, grad_s[2].U.E_total}; 303 304 for (CeedInt f = 0; f < 4; f++) { 305 for (CeedInt j = 0; j < dim; j++) grad_v[j][f][i] = 0; // No Change in density or momentum 306 v[f][i] = wdetJ * q_dot[f][i]; // K Mass/transient term 307 } 308 309 CeedScalar div_u = 0; 310 for (CeedInt j = 0; j < dim; j++) { 311 for (CeedInt k = 0; k < dim; k++) { 312 div_u += grad_s[k].Y.velocity[j]; 313 } 314 } 315 CeedScalar strong_conv = s.U.E_total * div_u + DotN(s.Y.velocity, Grad_E, dim); 316 CeedScalar strong_res = q_dot[4][i] + strong_conv; 317 318 v[4][i] = wdetJ * q_dot[4][i]; // transient part (ALWAYS) 319 320 CeedScalar uX[3] = {0.}; 321 MatVecNM(dXdx, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, uX); 322 323 if (context->strong_form) { // Strong Galerkin convection term: v div(E u) 324 v[4][i] += wdetJ * strong_conv; 325 } else { // Weak Galerkin convection term: -dv \cdot (E u) 326 for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = -wdetJ * s.U.E_total * uX[j]; 327 } 328 329 { // Diffusion 330 CeedScalar Fe[3], Fe_dXdx[3] = {0.}; 331 332 for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total; 333 MatVecNM(dXdx, Fe, dim, dim, CEED_NOTRANSPOSE, Fe_dXdx); 334 for (CeedInt k = 0; k < dim; k++) grad_v[k][4][i] -= wdetJ * Fe_dXdx[k]; 335 } 336 337 const CeedScalar TauS = Tau(context, s, dXdx, dim); 338 for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) { 339 case STAB_NONE: 340 break; 341 case STAB_SU: 342 grad_v[j][4][i] += wdetJ * TauS * strong_conv * uX[j]; 343 break; 344 case STAB_SUPG: 345 grad_v[j][4][i] += wdetJ * TauS * strong_res * uX[j]; 346 break; 347 } 348 StoredValuesPack(Q, i, 0, 14, zeros, jac_data); 349 } 350 } 351 352 CEED_QFUNCTION(IFunction_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 353 IFunction_AdvectionGeneric(ctx, Q, in, out, 3); 354 return 0; 355 } 356 357 CEED_QFUNCTION(IFunction_Advection2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 358 IFunction_AdvectionGeneric(ctx, Q, in, out, 2); 359 return 0; 360 } 361 362 CEED_QFUNCTION_HELPER void MassFunction_AdvectionGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, CeedInt dim) { 363 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 364 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 365 const CeedScalar(*q_data) = in[2]; 366 367 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 368 CeedScalar(*grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 369 370 AdvectionContext context = (AdvectionContext)ctx; 371 struct NewtonianIdealGasContext_ gas_struct = {0}; 372 NewtonianIdealGasContext gas = &gas_struct; 373 374 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 375 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 376 const State s = StateFromU(gas, qi); 377 CeedScalar wdetJ, dXdx[9]; 378 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 379 380 for (CeedInt f = 0; f < 4; f++) { 381 for (CeedInt j = 0; j < dim; j++) grad_v[j][f][i] = 0; // No Change in density or momentum 382 v[f][i] = wdetJ * q_dot[f][i]; // K Mass/transient term 383 } 384 385 // Unstabilized mass term 386 v[4][i] = wdetJ * q_dot[4][i]; 387 388 // Stabilized mass term 389 CeedScalar uX[3] = {0.}; 390 MatVecNM(dXdx, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, uX); 391 const CeedScalar TauS = Tau(context, s, dXdx, dim); 392 for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) { 393 case STAB_NONE: 394 case STAB_SU: 395 grad_v[j][4][i] = 0; 396 break; // These should be run with the unstabilized mass matrix anyways 397 case STAB_SUPG: 398 grad_v[j][4][i] = wdetJ * TauS * q_dot[4][i] * uX[j]; 399 break; 400 } 401 } 402 } 403 404 CEED_QFUNCTION(MassFunction_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 405 MassFunction_AdvectionGeneric(ctx, Q, in, out, 3); 406 return 0; 407 } 408 409 CEED_QFUNCTION(MassFunction_Advection2D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 410 MassFunction_AdvectionGeneric(ctx, Q, in, out, 2); 411 return 0; 412 } 413 414 // ***************************************************************************** 415 // This QFunction implements Advection for explicit time stepping method 416 // ***************************************************************************** 417 CEED_QFUNCTION_HELPER void RHSFunction_AdvectionGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, CeedInt dim) { 418 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 419 const CeedScalar(*grad_q) = in[1]; 420 const CeedScalar(*q_data) = in[2]; 421 422 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 423 CeedScalar(*grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 424 425 AdvectionContext context = (AdvectionContext)ctx; 426 struct NewtonianIdealGasContext_ gas_struct = {0}; 427 NewtonianIdealGasContext gas = &gas_struct; 428 429 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 430 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 431 const State s = StateFromU(gas, qi); 432 433 CeedScalar wdetJ, dXdx[9]; 434 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 435 State grad_s[3]; 436 StatePhysicalGradientFromReference_ND(dim, Q, i, gas, s, STATEVAR_CONSERVATIVE, grad_q, dXdx, grad_s); 437 438 const CeedScalar Grad_E[3] = {grad_s[0].U.E_total, grad_s[1].U.E_total, grad_s[2].U.E_total}; 439 440 for (CeedInt f = 0; f < 4; f++) { 441 for (CeedInt j = 0; j < dim; j++) grad_v[j][f][i] = 0; // No Change in density or momentum 442 v[f][i] = 0.; 443 } 444 445 CeedScalar div_u = 0; 446 for (CeedInt j = 0; j < dim; j++) { 447 for (CeedInt k = 0; k < dim; k++) { 448 div_u += grad_s[k].Y.velocity[j]; 449 } 450 } 451 CeedScalar strong_conv = s.U.E_total * div_u + DotN(s.Y.velocity, Grad_E, dim); 452 453 CeedScalar uX[3] = {0.}; 454 MatVecNM(dXdx, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, uX); 455 456 if (context->strong_form) { // Strong Galerkin convection term: v div(E u) 457 v[4][i] = -wdetJ * strong_conv; 458 for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = 0; 459 } else { // Weak Galerkin convection term: -dv \cdot (E u) 460 for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = wdetJ * s.U.E_total * uX[j]; 461 v[4][i] = 0.; 462 } 463 464 { // Diffusion 465 CeedScalar Fe[3], Fe_dXdx[3] = {0.}; 466 467 for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total; 468 MatVecNM(dXdx, Fe, dim, dim, CEED_NOTRANSPOSE, Fe_dXdx); 469 for (CeedInt k = 0; k < dim; k++) grad_v[k][4][i] += wdetJ * Fe_dXdx[k]; 470 } 471 472 const CeedScalar TauS = Tau(context, s, dXdx, dim); 473 for (CeedInt j = 0; j < dim; j++) switch (context->stabilization) { 474 case STAB_NONE: 475 break; 476 case STAB_SU: 477 case STAB_SUPG: 478 grad_v[j][4][i] -= wdetJ * TauS * strong_conv * uX[j]; 479 break; 480 } 481 } 482 } 483 484 CEED_QFUNCTION(RHS_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 485 RHSFunction_AdvectionGeneric(ctx, Q, in, out, 3); 486 return 0; 487 } 488 489 CEED_QFUNCTION(RHS_Advection2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 490 RHSFunction_AdvectionGeneric(ctx, Q, in, out, 2); 491 return 0; 492 } 493 494 // ***************************************************************************** 495 // This QFunction implements consistent outflow and inflow BCs 496 // for advection 497 // 498 // Inflow and outflow faces are determined based on sign(dot(wind, normal)): 499 // sign(dot(wind, normal)) > 0 : outflow BCs 500 // sign(dot(wind, normal)) < 0 : inflow BCs 501 // 502 // Outflow BCs: 503 // The validity of the weak form of the governing equations is extended to the outflow and the current values of E are applied. 504 // 505 // Inflow BCs: 506 // A prescribed Total Energy (E_wind) is applied weakly. 507 // ***************************************************************************** 508 CEED_QFUNCTION(Advection_InOutFlowGeneric)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, CeedInt dim) { 509 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 510 const CeedScalar(*q_data_sur) = in[2]; 511 512 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 513 AdvectionContext context = (AdvectionContext)ctx; 514 const CeedScalar E_wind = context->E_wind; 515 const CeedScalar strong_form = context->strong_form; 516 const bool is_implicit = context->implicit; 517 518 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 519 const CeedScalar rho = q[0][i]; 520 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 521 const CeedScalar E = q[4][i]; 522 523 CeedScalar wdetJb, norm[3]; 524 QdataBoundaryUnpack_ND(dim, Q, i, q_data_sur, &wdetJb, NULL, norm); 525 wdetJb *= is_implicit ? -1. : 1.; 526 527 const CeedScalar u_normal = DotN(norm, u, dim); 528 529 // No Change in density or momentum 530 for (CeedInt j = 0; j < 4; j++) { 531 v[j][i] = 0; 532 } 533 // Implementing in/outflow BCs 534 if (u_normal > 0) { // outflow 535 v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal; 536 } else { // inflow 537 v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal; 538 } 539 } 540 return 0; 541 } 542 543 CEED_QFUNCTION(Advection_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 544 Advection_InOutFlowGeneric(ctx, Q, in, out, 3); 545 return 0; 546 } 547 548 CEED_QFUNCTION(Advection2d_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 549 Advection_InOutFlowGeneric(ctx, Q, in, out, 2); 550 return 0; 551 } 552