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