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