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