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 CeedScalar *jac_data = out[2]; 271 272 AdvectionContext context = (AdvectionContext)ctx; 273 const CeedScalar CtauS = context->CtauS; 274 const CeedScalar strong_form = context->strong_form; 275 const CeedScalar zeros[14] = {0.}; 276 277 // Quadrature Point Loop 278 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 279 // Setup 280 // -- Interp in 281 const CeedScalar rho = q[0][i]; 282 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 283 const CeedScalar E = q[4][i]; 284 // -- Grad in 285 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]}; 286 const CeedScalar du[3][3] = { 287 {(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}, 288 {(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}, 289 {(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} 290 }; 291 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]}; 292 CeedScalar wdetJ, dXdx[3][3]; 293 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 294 // The Physics 295 // Note with the order that du was filled and the order that dXdx was filled 296 // du[j][k]= du_j / dX_K (note cap K to be clear this is u_{j,xi_k} ) 297 // dXdx[k][j] = dX_K / dx_j 298 // X_K=Kth reference element coordinate (note cap X and K instead of xi_k} 299 // x_j and u_j are jth physical position and velocity components 300 301 // No Change in density or momentum 302 for (CeedInt f = 0; f < 4; f++) { 303 for (CeedInt j = 0; j < 3; j++) dv[j][f][i] = 0; 304 v[f][i] = wdetJ * q_dot[f][i]; // K Mass/transient term 305 } 306 307 // -- Total Energy 308 // Evaluate the strong form using div(E u) = u . grad(E) + E div(u) 309 // or in index notation: (u_j E)_{,j} = u_j E_j + E u_{j,j} 310 CeedScalar div_u = 0, u_dot_grad_E = 0; 311 for (CeedInt j = 0; j < 3; j++) { 312 CeedScalar dEdx_j = 0; 313 for (CeedInt k = 0; k < 3; k++) { 314 div_u += du[j][k] * dXdx[k][j]; // u_{j,j} = u_{j,K} X_{K,j} 315 dEdx_j += dE[k] * dXdx[k][j]; 316 } 317 u_dot_grad_E += u[j] * dEdx_j; 318 } 319 CeedScalar strong_conv = E * div_u + u_dot_grad_E; 320 CeedScalar strong_res = q_dot[4][i] + strong_conv; 321 322 v[4][i] = wdetJ * q_dot[4][i]; // transient part (ALWAYS) 323 324 // Weak Galerkin convection term: -dv \cdot (E u) 325 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]); 326 327 // Strong Galerkin convection term: v div(E u) 328 v[4][i] += wdetJ * strong_form * strong_conv; 329 330 // Stabilization requires a measure of element transit time in the velocity 331 // field u. 332 CeedScalar uX[3]; 333 for (CeedInt j = 0; j < 3; j++) uX[j] = dXdx[j][0] * u[0] + dXdx[j][1] * u[1] + dXdx[j][2] * u[2]; 334 const CeedScalar TauS = CtauS / sqrt(uX[0] * uX[0] + uX[1] * uX[1] + uX[2] * uX[2]); 335 336 for (CeedInt j = 0; j < 3; j++) switch (context->stabilization) { 337 case STAB_NONE: 338 break; 339 case STAB_SU: 340 dv[j][4][i] += wdetJ * TauS * strong_conv * uX[j]; 341 break; 342 case STAB_SUPG: 343 dv[j][4][i] += wdetJ * TauS * strong_res * uX[j]; 344 break; 345 } 346 StoredValuesPack(Q, i, 0, 14, zeros, jac_data); 347 } // End Quadrature Point Loop 348 349 return 0; 350 } 351 352 // ***************************************************************************** 353 // This QFunction implements consistent outflow and inflow BCs 354 // for 3D advection 355 // 356 // Inflow and outflow faces are determined based on sign(dot(wind, normal)): 357 // sign(dot(wind, normal)) > 0 : outflow BCs 358 // sign(dot(wind, normal)) < 0 : inflow BCs 359 // 360 // Outflow BCs: 361 // The validity of the weak form of the governing equations is extended to the outflow and the current values of E are applied. 362 // 363 // Inflow BCs: 364 // A prescribed Total Energy (E_wind) is applied weakly. 365 // ***************************************************************************** 366 CEED_QFUNCTION(Advection_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 367 // Inputs 368 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 369 const CeedScalar(*q_data_sur) = in[2]; 370 371 // Outputs 372 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 373 AdvectionContext context = (AdvectionContext)ctx; 374 const CeedScalar E_wind = context->E_wind; 375 const CeedScalar strong_form = context->strong_form; 376 const bool is_implicit = context->implicit; 377 378 // Quadrature Point Loop 379 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 380 // Setup 381 // -- Interp in 382 const CeedScalar rho = q[0][i]; 383 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 384 const CeedScalar E = q[4][i]; 385 386 CeedScalar wdetJb, norm[3]; 387 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 388 wdetJb *= is_implicit ? -1. : 1.; 389 390 // Normal velocity 391 const CeedScalar u_normal = norm[0] * u[0] + norm[1] * u[1] + norm[2] * u[2]; 392 393 // No Change in density or momentum 394 for (CeedInt j = 0; j < 4; j++) { 395 v[j][i] = 0; 396 } 397 // Implementing in/outflow BCs 398 if (u_normal > 0) { // outflow 399 v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal; 400 } else { // inflow 401 v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal; 402 } 403 } // End Quadrature Point Loop 404 return 0; 405 } 406 // ***************************************************************************** 407 408 #endif // advection_h 409