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/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 104 switch (context->bubble_continuity_type) { 105 // original continuous, smooth shape 106 case ADVDIF_BUBBLE_CONTINUITY_SMOOTH: 107 q[4] = r <= rc ? (1. - r / rc) : 0.; 108 break; 109 // discontinuous, sharp back half shape 110 case ADVDIF_BUBBLE_CONTINUITY_BACK_SHARP: 111 q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) : 0.; 112 break; 113 // attempt to define a finite thickness that will get resolved under grid refinement 114 case ADVDIF_BUBBLE_CONTINUITY_THICK: 115 q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) * fmin(1.0, (center[1] - y) / 1.25) : 0.; 116 break; 117 case ADVDIF_BUBBLE_CONTINUITY_COSINE: 118 q[4] = r <= rc ? .5 + .5 * cos(r * M_PI / rc) : 0; 119 break; 120 } 121 break; 122 } 123 124 case ADVDIF_IC_COSINE_HILL: { 125 CeedScalar r = sqrt(Square(x - center[0]) + Square(y - center[1])); 126 CeedScalar half_width = context->lx / 2; 127 q[4] = r > half_width ? 0. : cos(2 * M_PI * r / half_width + M_PI) + 1.; 128 } break; 129 130 case ADVDIF_IC_SKEW: { 131 CeedScalar skewed_barrier[3] = {wind[0], wind[1], 0}; 132 CeedScalar inflow_to_point[3] = {x - context->lx / 2, y, 0}; 133 CeedScalar cross_product[3] = {0}; 134 const CeedScalar boundary_threshold = 20 * CEED_EPSILON; 135 Cross3(skewed_barrier, inflow_to_point, cross_product); 136 137 q[4] = cross_product[2] > boundary_threshold ? 0 : 1; 138 if ((x < boundary_threshold && wind[0] < boundary_threshold) || // outflow at -x boundary 139 (y < boundary_threshold && wind[1] < boundary_threshold) || // outflow at -y boundary 140 (x > context->lx - boundary_threshold && wind[0] > boundary_threshold) || // outflow at +x boundary 141 (y > context->ly - boundary_threshold && wind[1] > boundary_threshold) // outflow at +y boundary 142 ) { 143 q[4] = 0; 144 } 145 } break; 146 147 case ADVDIF_IC_WAVE: { 148 CeedScalar theta = context->wave_frequency * DotN(X, wind, dim) + context->wave_phase; 149 switch (context->wave_type) { 150 case ADVDIF_WAVE_SINE: 151 q[4] = sin(theta); 152 break; 153 case ADVDIF_WAVE_SQUARE: 154 q[4] = sin(theta) > 100 * CEED_EPSILON ? 1 : -1; 155 break; 156 } 157 } break; 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, NewtonianIdealGasContext 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 NewtonianIdealGasContext gas; 279 struct NewtonianIdealGasContext_ gas_struct = {0}; 280 gas = &gas_struct; 281 282 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 283 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 284 const State s = StateFromU(gas, qi); 285 286 CeedScalar wdetJ, dXdx[9]; 287 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 288 State grad_s[3]; 289 StatePhysicalGradientFromReference_ND(dim, Q, i, gas, s, STATEVAR_CONSERVATIVE, grad_q, dXdx, grad_s); 290 291 const CeedScalar Grad_E[3] = {grad_s[0].U.E_total, grad_s[1].U.E_total, grad_s[2].U.E_total}; 292 293 for (CeedInt f = 0; f < 4; f++) { 294 for (CeedInt j = 0; j < dim; j++) grad_v[j][f][i] = 0; // No Change in density or momentum 295 v[f][i] = wdetJ * q_dot[f][i]; // K Mass/transient term 296 } 297 298 CeedScalar div_u = 0; 299 for (CeedInt j = 0; j < dim; j++) { 300 for (CeedInt k = 0; k < dim; k++) { 301 div_u += grad_s[k].Y.velocity[j]; 302 } 303 } 304 CeedScalar uX[3] = {0.}; 305 MatVecNM(dXdx, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, uX); 306 CeedScalar strong_conv = s.U.E_total * div_u + DotN(s.Y.velocity, Grad_E, dim); 307 308 v[4][i] = wdetJ * q_dot[4][i]; // transient part (ALWAYS) 309 if (context->strong_form) { // Strong Galerkin convection term: v div(E u) 310 v[4][i] += wdetJ * strong_conv; 311 } else { // Weak Galerkin convection term: -dv \cdot (E u) 312 for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = -wdetJ * s.U.E_total * uX[j]; 313 } 314 315 { // Diffusion 316 CeedScalar Fe[3], Fe_dXdx[3] = {0.}; 317 318 for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total; 319 MatVecNM(dXdx, Fe, dim, dim, CEED_NOTRANSPOSE, Fe_dXdx); 320 for (CeedInt k = 0; k < dim; k++) grad_v[k][4][i] -= wdetJ * Fe_dXdx[k]; 321 } 322 323 const CeedScalar TauS = Tau(context, s, dXdx, dim); 324 for (CeedInt j = 0; j < dim; j++) { 325 switch (context->stabilization) { 326 case STAB_NONE: 327 break; 328 case STAB_SU: 329 grad_v[j][4][i] += wdetJ * TauS * uX[j] * strong_conv; 330 break; 331 case STAB_SUPG: { 332 CeedScalar divFdiff_i = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE ? divFdiff[i] : 0.; 333 grad_v[j][4][i] += wdetJ * TauS * uX[j] * (q_dot[4][i] + strong_conv + divFdiff_i); 334 } break; 335 } 336 } 337 } 338 return 0; 339 } 340 341 CEED_QFUNCTION(IFunction_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 342 return IFunction_AdvectionGeneric(ctx, Q, in, out, 3); 343 } 344 345 CEED_QFUNCTION(IFunction_Advection2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 346 return IFunction_AdvectionGeneric(ctx, Q, in, out, 2); 347 } 348 349 CEED_QFUNCTION_HELPER int MassFunction_AdvectionGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, CeedInt dim) { 350 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 351 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 352 const CeedScalar(*q_data) = in[2]; 353 354 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 355 CeedScalar(*grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 356 357 AdvectionContext context = (AdvectionContext)ctx; 358 struct NewtonianIdealGasContext_ gas_struct = {0}; 359 NewtonianIdealGasContext gas = &gas_struct; 360 361 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 362 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 363 const State s = StateFromU(gas, qi); 364 CeedScalar wdetJ, dXdx[9]; 365 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 366 367 for (CeedInt f = 0; f < 4; f++) { 368 for (CeedInt j = 0; j < dim; j++) grad_v[j][f][i] = 0; // No Change in density or momentum 369 v[f][i] = wdetJ * q_dot[f][i]; // K Mass/transient term 370 } 371 372 // Unstabilized mass term 373 v[4][i] = wdetJ * q_dot[4][i]; 374 375 // Stabilized mass term 376 CeedScalar uX[3] = {0.}; 377 MatVecNM(dXdx, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, uX); 378 const CeedScalar TauS = Tau(context, s, dXdx, dim); 379 for (CeedInt j = 0; j < dim; j++) { 380 switch (context->stabilization) { 381 case STAB_NONE: 382 case STAB_SU: 383 grad_v[j][4][i] = 0; 384 break; // These should be run with the unstabilized mass matrix anyways 385 case STAB_SUPG: 386 grad_v[j][4][i] = wdetJ * TauS * q_dot[4][i] * uX[j]; 387 break; 388 } 389 } 390 } 391 return 0; 392 } 393 394 CEED_QFUNCTION(MassFunction_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 395 return MassFunction_AdvectionGeneric(ctx, Q, in, out, 3); 396 } 397 398 CEED_QFUNCTION(MassFunction_Advection2D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 399 return MassFunction_AdvectionGeneric(ctx, Q, in, out, 2); 400 } 401 402 // ***************************************************************************** 403 // This QFunction implements Advection for explicit time stepping method 404 // ***************************************************************************** 405 CEED_QFUNCTION_HELPER int RHSFunction_AdvectionGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, CeedInt dim) { 406 AdvectionContext context = (AdvectionContext)ctx; 407 408 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 409 const CeedScalar(*grad_q) = in[1]; 410 const CeedScalar(*q_data) = in[2]; 411 const CeedScalar(*divFdiff) = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE ? in[4] : NULL; 412 413 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 414 CeedScalar(*grad_v)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 415 416 struct NewtonianIdealGasContext_ gas_struct = {0}; 417 NewtonianIdealGasContext gas = &gas_struct; 418 419 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 420 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 421 const State s = StateFromU(gas, qi); 422 423 CeedScalar wdetJ, dXdx[9]; 424 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 425 State grad_s[3]; 426 StatePhysicalGradientFromReference_ND(dim, Q, i, gas, s, STATEVAR_CONSERVATIVE, grad_q, dXdx, grad_s); 427 428 const CeedScalar Grad_E[3] = {grad_s[0].U.E_total, grad_s[1].U.E_total, grad_s[2].U.E_total}; 429 430 for (CeedInt f = 0; f < 4; f++) { 431 for (CeedInt j = 0; j < dim; j++) grad_v[j][f][i] = 0; // No Change in density or momentum 432 v[f][i] = 0.; 433 } 434 435 CeedScalar div_u = 0; 436 for (CeedInt j = 0; j < dim; j++) { 437 for (CeedInt k = 0; k < dim; k++) { 438 div_u += grad_s[k].Y.velocity[j]; 439 } 440 } 441 CeedScalar strong_conv = s.U.E_total * div_u + DotN(s.Y.velocity, Grad_E, dim); 442 443 CeedScalar uX[3] = {0.}; 444 MatVecNM(dXdx, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, uX); 445 446 if (context->strong_form) { // Strong Galerkin convection term: v div(E u) 447 v[4][i] = -wdetJ * strong_conv; 448 for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = 0; 449 } else { // Weak Galerkin convection term: -dv \cdot (E u) 450 for (CeedInt j = 0; j < dim; j++) grad_v[j][4][i] = wdetJ * s.U.E_total * uX[j]; 451 v[4][i] = 0.; 452 } 453 454 { // Diffusion 455 CeedScalar Fe[3], Fe_dXdx[3] = {0.}; 456 457 for (CeedInt i = 0; i < dim; i++) Fe[i] = -context->diffusion_coeff * grad_s[i].U.E_total; 458 MatVecNM(dXdx, Fe, dim, dim, CEED_NOTRANSPOSE, Fe_dXdx); 459 for (CeedInt k = 0; k < dim; k++) grad_v[k][4][i] += wdetJ * Fe_dXdx[k]; 460 } 461 462 const CeedScalar TauS = Tau(context, s, dXdx, dim); 463 for (CeedInt j = 0; j < dim; j++) { 464 switch (context->stabilization) { 465 case STAB_NONE: 466 break; 467 case STAB_SU: 468 case STAB_SUPG: { 469 CeedScalar divFdiff_i = context->divFdiff_method != DIV_DIFF_FLUX_PROJ_NONE ? divFdiff[i] : 0.; 470 grad_v[j][4][i] -= wdetJ * TauS * (strong_conv + divFdiff_i) * uX[j]; 471 } break; 472 } 473 } 474 } 475 return 0; 476 } 477 478 CEED_QFUNCTION(RHS_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 479 return RHSFunction_AdvectionGeneric(ctx, Q, in, out, 3); 480 } 481 482 CEED_QFUNCTION(RHS_Advection2d)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 483 return RHSFunction_AdvectionGeneric(ctx, Q, in, out, 2); 484 } 485 486 // ***************************************************************************** 487 // This QFunction implements consistent outflow and inflow BCs 488 // for advection 489 // 490 // Inflow and outflow faces are determined based on sign(dot(wind, normal)): 491 // sign(dot(wind, normal)) > 0 : outflow BCs 492 // sign(dot(wind, normal)) < 0 : inflow BCs 493 // 494 // Outflow BCs: 495 // The validity of the weak form of the governing equations is extended to the outflow and the current values of E are applied. 496 // 497 // Inflow BCs: 498 // A prescribed Total Energy (E_wind) is applied weakly. 499 // ***************************************************************************** 500 CEED_QFUNCTION_HELPER int Advection_InOutFlowGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, CeedInt dim) { 501 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 502 const CeedScalar(*q_data_sur) = in[2]; 503 504 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 505 AdvectionContext context = (AdvectionContext)ctx; 506 const CeedScalar E_wind = context->E_wind; 507 const CeedScalar strong_form = context->strong_form; 508 const bool is_implicit = context->implicit; 509 510 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 511 const CeedScalar rho = q[0][i]; 512 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 513 const CeedScalar E = q[4][i]; 514 515 CeedScalar wdetJb, normal[3]; 516 QdataBoundaryUnpack_ND(dim, Q, i, q_data_sur, &wdetJb, NULL, normal); 517 wdetJb *= is_implicit ? -1. : 1.; 518 519 const CeedScalar u_normal = DotN(normal, u, dim); 520 521 // No Change in density or momentum 522 for (CeedInt j = 0; j < 4; j++) { 523 v[j][i] = 0; 524 } 525 // Implementing in/outflow BCs 526 if (u_normal > 0) { // outflow 527 v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal; 528 } else { // inflow 529 v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal; 530 } 531 } 532 return 0; 533 } 534 535 CEED_QFUNCTION(Advection_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 536 return Advection_InOutFlowGeneric(ctx, Q, in, out, 3); 537 } 538 539 CEED_QFUNCTION(Advection2d_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 540 return Advection_InOutFlowGeneric(ctx, Q, in, out, 2); 541 } 542 543 // @brief Volume integral for RHS of divergence of diffusive flux direct projection 544 CEED_QFUNCTION_HELPER int DivDiffusiveFluxVolumeRHS_AdvDif_Generic(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 545 const CeedInt dim) { 546 const CeedScalar(*Grad_q) = in[0]; 547 const CeedScalar(*q_data) = in[1]; 548 CeedScalar(*Grad_v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 549 550 AdvectionContext context = (AdvectionContext)ctx; 551 552 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 553 CeedScalar wdetJ, dXdx[9], F_diff[3] = {0.}; 554 555 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 556 { // Get physical diffusive flux 557 CeedScalar Grad_qn[15], grad_E_ref[3]; 558 559 GradUnpackN(Q, i, 5, dim, Grad_q, Grad_qn); 560 CopyN(&Grad_qn[4 * dim], grad_E_ref, dim); 561 MatVecNM(dXdx, grad_E_ref, dim, dim, CEED_NOTRANSPOSE, F_diff); 562 ScaleN(F_diff, -context->diffusion_coeff, dim); 563 } 564 565 CeedScalar F_diff_dXdx[3] = {0.}; 566 MatVecNM(dXdx, F_diff, dim, dim, CEED_NOTRANSPOSE, F_diff_dXdx); 567 for (CeedInt k = 0; k < dim; k++) Grad_v[k][i] = -wdetJ * F_diff_dXdx[k]; 568 } 569 return 0; 570 } 571 572 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_AdvDif_2D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 573 return DivDiffusiveFluxVolumeRHS_AdvDif_Generic(ctx, Q, in, out, 2); 574 } 575 576 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_AdvDif_3D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 577 return DivDiffusiveFluxVolumeRHS_AdvDif_Generic(ctx, Q, in, out, 3); 578 } 579 580 // @brief Boundary integral for RHS of divergence of diffusive flux direct projection 581 CEED_QFUNCTION_HELPER int DivDiffusiveFluxBoundaryRHS_AdvDif_Generic(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 582 const CeedInt dim) { 583 const CeedScalar(*Grad_q) = in[0]; 584 const CeedScalar(*q_data) = in[1]; 585 CeedScalar(*v) = out[0]; 586 587 AdvectionContext context = (AdvectionContext)ctx; 588 589 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 590 CeedScalar wdetJ, normal[3], dXdx[9], F_diff[3] = {0.}; 591 592 QdataBoundaryGradientUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx, normal); 593 { // Get physical diffusive flux 594 CeedScalar Grad_qn[15], grad_E_ref[3]; 595 596 GradUnpackN(Q, i, 5, dim, Grad_q, Grad_qn); 597 CopyN(&Grad_qn[4 * dim], grad_E_ref, dim); 598 MatVecNM(dXdx, grad_E_ref, dim, dim, CEED_NOTRANSPOSE, F_diff); 599 ScaleN(F_diff, -context->diffusion_coeff, dim); 600 } 601 602 v[i] = wdetJ * DotN(F_diff, normal, dim); 603 } 604 return 0; 605 } 606 607 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_AdvDif_2D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 608 return DivDiffusiveFluxBoundaryRHS_AdvDif_Generic(ctx, Q, in, out, 2); 609 } 610 611 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_AdvDif_3D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 612 return DivDiffusiveFluxBoundaryRHS_AdvDif_Generic(ctx, Q, in, out, 3); 613 } 614 615 // @brief Volume integral for RHS of diffusive flux indirect projection 616 CEED_QFUNCTION_HELPER int DiffusiveFluxRHS_AdvDif_Generic(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 617 const CeedInt dim) { 618 const CeedScalar(*Grad_q) = in[0]; 619 const CeedScalar(*q_data) = in[1]; 620 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 621 622 AdvectionContext context = (AdvectionContext)ctx; 623 624 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 625 CeedScalar wdetJ, dXdx[9], F_diff[3] = {0.}; 626 627 QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx); 628 { // Get physical diffusive flux 629 CeedScalar Grad_qn[15], grad_E_ref[3]; 630 631 GradUnpackN(Q, i, 5, dim, Grad_q, Grad_qn); 632 CopyN(&Grad_qn[4 * dim], grad_E_ref, dim); 633 MatVecNM(dXdx, grad_E_ref, dim, dim, CEED_NOTRANSPOSE, F_diff); 634 ScaleN(F_diff, -context->diffusion_coeff, dim); 635 } 636 for (CeedInt k = 0; k < dim; k++) v[k][i] = wdetJ * F_diff[k]; 637 } 638 return 0; 639 } 640 641 CEED_QFUNCTION(DiffusiveFluxRHS_AdvDif_2D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 642 return DiffusiveFluxRHS_AdvDif_Generic(ctx, Q, in, out, 2); 643 } 644 645 CEED_QFUNCTION(DiffusiveFluxRHS_AdvDif_3D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 646 return DiffusiveFluxRHS_AdvDif_Generic(ctx, Q, in, out, 3); 647 } 648