1 // Copyright (c) 2017-2022, 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 "utils.h" 18 19 typedef struct SetupContextAdv_ *SetupContextAdv; 20 struct SetupContextAdv_ { 21 CeedScalar rc; 22 CeedScalar lx; 23 CeedScalar ly; 24 CeedScalar lz; 25 CeedScalar wind[3]; 26 CeedScalar time; 27 int wind_type; // See WindType: 0=ROTATION, 1=TRANSLATION 28 int bubble_type; // See BubbleType: 0=SPHERE, 1=CYLINDER 29 int bubble_continuity_type; // See BubbleContinuityType: 0=SMOOTH, 1=BACK_SHARP 2=THICK 30 }; 31 32 typedef struct AdvectionContext_ *AdvectionContext; 33 struct AdvectionContext_ { 34 CeedScalar CtauS; 35 CeedScalar strong_form; 36 CeedScalar E_wind; 37 bool implicit; 38 int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG 39 }; 40 41 // ***************************************************************************** 42 // This QFunction sets the initial conditions and the boundary conditions 43 // for two test cases: ROTATION and TRANSLATION 44 // 45 // -- ROTATION (default) 46 // Initial Conditions: 47 // Mass Density: 48 // Constant mass density of 1.0 49 // Momentum Density: 50 // Rotational 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 // 0.0 flux 62 // 63 // -- TRANSLATION 64 // Initial Conditions: 65 // Mass Density: 66 // Constant mass density of 1.0 67 // Momentum Density: 68 // Constant rectilinear field in x,y 69 // Energy Density: 70 // Maximum of 1. x0 decreasing linearly to 0. as radial distance 71 // increases to (1.-r/rc), then 0. everywhere else 72 // 73 // Boundary Conditions: 74 // Mass Density: 75 // 0.0 flux 76 // Momentum Density: 77 // 0.0 78 // Energy Density: 79 // Inflow BCs: 80 // E = E_wind 81 // Outflow BCs: 82 // E = E(boundary) 83 // Both In/Outflow BCs for E are applied weakly in the 84 // QFunction "Advection_Sur" 85 // 86 // ***************************************************************************** 87 88 // ***************************************************************************** 89 // This helper function provides support for the exact, time-dependent solution (currently not implemented) and IC formulation for 3D advection 90 // ***************************************************************************** 91 CEED_QFUNCTION_HELPER CeedInt Exact_Advection(CeedInt dim, CeedScalar time, const CeedScalar X[], CeedInt Nf, CeedScalar q[], void *ctx) { 92 const SetupContextAdv context = (SetupContextAdv)ctx; 93 const CeedScalar rc = context->rc; 94 const CeedScalar lx = context->lx; 95 const CeedScalar ly = context->ly; 96 const CeedScalar lz = context->lz; 97 const CeedScalar *wind = context->wind; 98 99 // Setup 100 const CeedScalar x0[3] = {0.25 * lx, 0.5 * ly, 0.5 * lz}; 101 const CeedScalar center[3] = {0.5 * lx, 0.5 * ly, 0.5 * lz}; 102 103 // -- Coordinates 104 const CeedScalar x = X[0]; 105 const CeedScalar y = X[1]; 106 const CeedScalar z = X[2]; 107 108 // -- Energy 109 CeedScalar r = 0.; 110 switch (context->bubble_type) { 111 // original sphere 112 case 0: { // (dim=3) 113 r = sqrt(Square(x - x0[0]) + Square(y - x0[1]) + Square(z - x0[2])); 114 } break; 115 // cylinder (needs periodicity to work properly) 116 case 1: { // (dim=2) 117 r = sqrt(Square(x - x0[0]) + Square(y - x0[1])); 118 } break; 119 } 120 121 // Initial Conditions 122 switch (context->wind_type) { 123 case 0: // Rotation 124 q[0] = 1.; 125 q[1] = -(y - center[1]); 126 q[2] = (x - center[0]); 127 q[3] = 0; 128 break; 129 case 1: // Translation 130 q[0] = 1.; 131 q[1] = wind[0]; 132 q[2] = wind[1]; 133 q[3] = wind[2]; 134 break; 135 } 136 137 switch (context->bubble_continuity_type) { 138 // original continuous, smooth shape 139 case 0: { 140 q[4] = r <= rc ? (1. - r / rc) : 0.; 141 } break; 142 // discontinuous, sharp back half shape 143 case 1: { 144 q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) : 0.; 145 } break; 146 // attempt to define a finite thickness that will get resolved under grid refinement 147 case 2: { 148 q[4] = ((r <= rc) && (y < center[1])) ? (1. - r / rc) * fmin(1.0, (center[1] - y) / 1.25) : 0.; 149 } break; 150 } 151 return 0; 152 } 153 154 // ***************************************************************************** 155 // This QFunction sets the initial conditions for 3D advection 156 // ***************************************************************************** 157 CEED_QFUNCTION(ICsAdvection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 158 // Inputs 159 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 160 // Outputs 161 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 162 163 // Quadrature Point Loop 164 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 165 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 166 CeedScalar q[5] = {0.}; 167 168 Exact_Advection(3, 0., x, 5, q, ctx); 169 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 170 } // End of Quadrature Point Loop 171 172 // Return 173 return 0; 174 } 175 176 // ***************************************************************************** 177 // This QFunction implements the following formulation of the advection equation 178 // 179 // This is 3D advection given in two formulations based upon the weak form. 180 // 181 // State Variables: q = ( rho, U1, U2, U3, E ) 182 // rho - Mass Density 183 // Ui - Momentum Density , Ui = rho ui 184 // E - Total Energy Density 185 // 186 // Advection Equation: 187 // dE/dt + div( E u ) = 0 188 // ***************************************************************************** 189 CEED_QFUNCTION(Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 190 // Inputs 191 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 192 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 193 const CeedScalar(*q_data) = in[2]; 194 195 // Outputs 196 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 197 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 198 199 // Context 200 AdvectionContext context = (AdvectionContext)ctx; 201 const CeedScalar CtauS = context->CtauS; 202 const CeedScalar strong_form = context->strong_form; 203 204 // Quadrature Point Loop 205 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 206 // Setup 207 // -- Interp in 208 const CeedScalar rho = q[0][i]; 209 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 210 const CeedScalar E = q[4][i]; 211 // -- Grad in 212 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]}; 213 const CeedScalar du[3][3] = { 214 {(dq[0][1][i] - drho[0] * u[0]) / rho, (dq[1][1][i] - drho[1] * u[0]) / rho, (dq[2][1][i] - drho[2] * u[0]) / rho}, 215 {(dq[0][2][i] - drho[0] * u[1]) / rho, (dq[1][2][i] - drho[1] * u[1]) / rho, (dq[2][2][i] - drho[2] * u[1]) / rho}, 216 {(dq[0][3][i] - drho[0] * u[2]) / rho, (dq[1][3][i] - drho[1] * u[2]) / rho, (dq[2][3][i] - drho[2] * u[2]) / rho} 217 }; 218 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]}; 219 CeedScalar wdetJ, dXdx[3][3]; 220 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 221 // The Physics 222 // Note with the order that du was filled and the order that dXdx was filled 223 // du[j][k]= du_j / dX_K (note cap K to be clear this is u_{j,xi_k}) 224 // dXdx[k][j] = dX_K / dx_j 225 // X_K=Kth reference element coordinate (note cap X and K instead of xi_k} 226 // x_j and u_j are jth physical position and velocity components 227 228 // No Change in density or momentum 229 for (CeedInt f = 0; f < 4; f++) { 230 for (CeedInt j = 0; j < 3; j++) dv[j][f][i] = 0; 231 v[f][i] = 0; 232 } 233 234 // -- Total Energy 235 // Evaluate the strong form using div(E u) = u . grad(E) + E div(u) 236 // or in index notation: (u_j E)_{,j} = u_j E_j + E u_{j,j} 237 CeedScalar div_u = 0, u_dot_grad_E = 0; 238 for (CeedInt j = 0; j < 3; j++) { 239 CeedScalar dEdx_j = 0; 240 for (CeedInt k = 0; k < 3; k++) { 241 div_u += du[j][k] * dXdx[k][j]; // u_{j,j} = u_{j,K} X_{K,j} 242 dEdx_j += dE[k] * dXdx[k][j]; 243 } 244 u_dot_grad_E += u[j] * dEdx_j; 245 } 246 CeedScalar strong_conv = E * div_u + u_dot_grad_E; 247 248 // Weak Galerkin convection term: dv \cdot (E u) 249 for (CeedInt j = 0; j < 3; j++) dv[j][4][i] = (1 - strong_form) * wdetJ * E * (u[0] * dXdx[j][0] + u[1] * dXdx[j][1] + u[2] * dXdx[j][2]); 250 v[4][i] = 0; 251 252 // Strong Galerkin convection term: - v div(E u) 253 v[4][i] = -strong_form * wdetJ * strong_conv; 254 255 // Stabilization requires a measure of element transit time in the velocity 256 // field u. 257 CeedScalar uX[3]; 258 for (CeedInt j = 0; j < 3; j++) uX[j] = dXdx[j][0] * u[0] + dXdx[j][1] * u[1] + dXdx[j][2] * u[2]; 259 const CeedScalar TauS = CtauS / sqrt(uX[0] * uX[0] + uX[1] * uX[1] + uX[2] * uX[2]); 260 for (CeedInt j = 0; j < 3; j++) dv[j][4][i] -= wdetJ * TauS * strong_conv * uX[j]; 261 } // End Quadrature Point Loop 262 263 return 0; 264 } 265 266 // ***************************************************************************** 267 // This QFunction implements 3D (mentioned above) with implicit time stepping method 268 // ***************************************************************************** 269 CEED_QFUNCTION(IFunction_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 270 // Inputs 271 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 272 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 273 const CeedScalar(*q_dot)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 274 const CeedScalar(*q_data) = in[3]; 275 276 // Outputs 277 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 278 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 279 280 AdvectionContext context = (AdvectionContext)ctx; 281 const CeedScalar CtauS = context->CtauS; 282 const CeedScalar strong_form = context->strong_form; 283 284 // Quadrature Point Loop 285 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 286 // Setup 287 // -- Interp in 288 const CeedScalar rho = q[0][i]; 289 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 290 const CeedScalar E = q[4][i]; 291 // -- Grad in 292 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]}; 293 const CeedScalar du[3][3] = { 294 {(dq[0][1][i] - drho[0] * u[0]) / rho, (dq[1][1][i] - drho[1] * u[0]) / rho, (dq[2][1][i] - drho[2] * u[0]) / rho}, 295 {(dq[0][2][i] - drho[0] * u[1]) / rho, (dq[1][2][i] - drho[1] * u[1]) / rho, (dq[2][2][i] - drho[2] * u[1]) / rho}, 296 {(dq[0][3][i] - drho[0] * u[2]) / rho, (dq[1][3][i] - drho[1] * u[2]) / rho, (dq[2][3][i] - drho[2] * u[2]) / rho} 297 }; 298 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]}; 299 CeedScalar wdetJ, dXdx[3][3]; 300 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 301 // The Physics 302 // Note with the order that du was filled and the order that dXdx was filled 303 // du[j][k]= du_j / dX_K (note cap K to be clear this is u_{j,xi_k} ) 304 // dXdx[k][j] = dX_K / dx_j 305 // X_K=Kth reference element coordinate (note cap X and K instead of xi_k} 306 // x_j and u_j are jth physical position and velocity components 307 308 // No Change in density or momentum 309 for (CeedInt f = 0; f < 4; f++) { 310 for (CeedInt j = 0; j < 3; j++) dv[j][f][i] = 0; 311 v[f][i] = wdetJ * q_dot[f][i]; // K Mass/transient term 312 } 313 314 // -- Total Energy 315 // Evaluate the strong form using div(E u) = u . grad(E) + E div(u) 316 // or in index notation: (u_j E)_{,j} = u_j E_j + E u_{j,j} 317 CeedScalar div_u = 0, u_dot_grad_E = 0; 318 for (CeedInt j = 0; j < 3; j++) { 319 CeedScalar dEdx_j = 0; 320 for (CeedInt k = 0; k < 3; k++) { 321 div_u += du[j][k] * dXdx[k][j]; // u_{j,j} = u_{j,K} X_{K,j} 322 dEdx_j += dE[k] * dXdx[k][j]; 323 } 324 u_dot_grad_E += u[j] * dEdx_j; 325 } 326 CeedScalar strong_conv = E * div_u + u_dot_grad_E; 327 CeedScalar strong_res = q_dot[4][i] + strong_conv; 328 329 v[4][i] = wdetJ * q_dot[4][i]; // transient part (ALWAYS) 330 331 // Weak Galerkin convection term: -dv \cdot (E u) 332 for (CeedInt j = 0; j < 3; j++) dv[j][4][i] = -wdetJ * (1 - strong_form) * E * (u[0] * dXdx[j][0] + u[1] * dXdx[j][1] + u[2] * dXdx[j][2]); 333 334 // Strong Galerkin convection term: v div(E u) 335 v[4][i] += wdetJ * strong_form * strong_conv; 336 337 // Stabilization requires a measure of element transit time in the velocity 338 // field u. 339 CeedScalar uX[3]; 340 for (CeedInt j = 0; j < 3; j++) uX[j] = dXdx[j][0] * u[0] + dXdx[j][1] * u[1] + dXdx[j][2] * u[2]; 341 const CeedScalar TauS = CtauS / sqrt(uX[0] * uX[0] + uX[1] * uX[1] + uX[2] * uX[2]); 342 343 for (CeedInt j = 0; j < 3; j++) switch (context->stabilization) { 344 case 0: 345 break; 346 case 1: 347 dv[j][4][i] += wdetJ * TauS * strong_conv * uX[j]; // SU 348 break; 349 case 2: 350 dv[j][4][i] += wdetJ * TauS * strong_res * uX[j]; // SUPG 351 break; 352 } 353 } // End Quadrature Point Loop 354 355 return 0; 356 } 357 358 // ***************************************************************************** 359 // This QFunction implements consistent outflow and inflow BCs 360 // for 3D advection 361 // 362 // Inflow and outflow faces are determined based on sign(dot(wind, normal)): 363 // sign(dot(wind, normal)) > 0 : outflow BCs 364 // sign(dot(wind, normal)) < 0 : inflow BCs 365 // 366 // Outflow BCs: 367 // The validity of the weak form of the governing equations is extended to the outflow and the current values of E are applied. 368 // 369 // Inflow BCs: 370 // A prescribed Total Energy (E_wind) is applied weakly. 371 // ***************************************************************************** 372 CEED_QFUNCTION(Advection_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 373 // Inputs 374 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 375 const CeedScalar(*q_data_sur) = in[2]; 376 377 // Outputs 378 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 379 AdvectionContext context = (AdvectionContext)ctx; 380 const CeedScalar E_wind = context->E_wind; 381 const CeedScalar strong_form = context->strong_form; 382 const bool is_implicit = context->implicit; 383 384 // Quadrature Point Loop 385 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 386 // Setup 387 // -- Interp in 388 const CeedScalar rho = q[0][i]; 389 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 390 const CeedScalar E = q[4][i]; 391 392 CeedScalar wdetJb, norm[3]; 393 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 394 wdetJb *= is_implicit ? -1. : 1.; 395 396 // Normal velocity 397 const CeedScalar u_normal = norm[0] * u[0] + norm[1] * u[1] + norm[2] * u[2]; 398 399 // No Change in density or momentum 400 for (CeedInt j = 0; j < 4; j++) { 401 v[j][i] = 0; 402 } 403 // Implementing in/outflow BCs 404 if (u_normal > 0) { // outflow 405 v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal; 406 } else { // inflow 407 v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal; 408 } 409 } // End Quadrature Point Loop 410 return 0; 411 } 412 // ***************************************************************************** 413 414 #endif // advection_h 415