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