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