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