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