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 "advection_generic.h" 18 #include "advection_types.h" 19 #include "newtonian_state.h" 20 #include "newtonian_types.h" 21 #include "stabilization_types.h" 22 #include "utils.h" 23 24 // ***************************************************************************** 25 // This QFunction sets the initial conditions for 3D advection 26 // ***************************************************************************** 27 CEED_QFUNCTION(ICsAdvection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 28 const CeedScalar(*X)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 29 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 30 31 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 32 const CeedScalar x[] = {X[0][i], X[1][i], X[2][i]}; 33 CeedScalar q[5] = {0.}; 34 35 Exact_AdvectionGeneric(3, 0., x, 5, q, ctx); 36 for (CeedInt j = 0; j < 5; j++) q0[j][i] = q[j]; 37 } 38 return 0; 39 } 40 41 // ***************************************************************************** 42 // This QFunction implements the following formulation of the advection equation 43 // 44 // This is 3D advection given in two formulations based upon the weak form. 45 // 46 // State Variables: q = ( rho, U1, U2, U3, E ) 47 // rho - Mass Density 48 // Ui - Momentum Density , Ui = rho ui 49 // E - Total Energy Density 50 // 51 // Advection Equation: 52 // dE/dt + div( E u ) = 0 53 // ***************************************************************************** 54 CEED_QFUNCTION(Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 55 // Inputs 56 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 57 const CeedScalar(*dq)[5][CEED_Q_VLA] = (const CeedScalar(*)[5][CEED_Q_VLA])in[1]; 58 const CeedScalar(*q_data) = in[2]; 59 60 // Outputs 61 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 62 CeedScalar(*dv)[5][CEED_Q_VLA] = (CeedScalar(*)[5][CEED_Q_VLA])out[1]; 63 64 // Context 65 AdvectionContext context = (AdvectionContext)ctx; 66 const CeedScalar CtauS = context->CtauS; 67 const CeedScalar strong_form = context->strong_form; 68 69 // Quadrature Point Loop 70 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 71 // Setup 72 // -- Interp in 73 const CeedScalar rho = q[0][i]; 74 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 75 const CeedScalar E = q[4][i]; 76 // -- Grad in 77 const CeedScalar drho[3] = {dq[0][0][i], dq[1][0][i], dq[2][0][i]}; 78 const CeedScalar du[3][3] = { 79 {(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}, 80 {(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}, 81 {(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} 82 }; 83 const CeedScalar dE[3] = {dq[0][4][i], dq[1][4][i], dq[2][4][i]}; 84 CeedScalar wdetJ, dXdx[3][3]; 85 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 86 // The Physics 87 // Note with the order that du was filled and the order that dXdx was filled 88 // du[j][k]= du_j / dX_K (note cap K to be clear this is u_{j,xi_k}) 89 // dXdx[k][j] = dX_K / dx_j 90 // X_K=Kth reference element coordinate (note cap X and K instead of xi_k} 91 // x_j and u_j are jth physical position and velocity components 92 93 // No Change in density or momentum 94 for (CeedInt f = 0; f < 4; f++) { 95 for (CeedInt j = 0; j < 3; j++) dv[j][f][i] = 0; 96 v[f][i] = 0; 97 } 98 99 // -- Total Energy 100 // Evaluate the strong form using div(E u) = u . grad(E) + E div(u) 101 // or in index notation: (u_j E)_{,j} = u_j E_j + E u_{j,j} 102 CeedScalar div_u = 0, u_dot_grad_E = 0; 103 for (CeedInt j = 0; j < 3; j++) { 104 CeedScalar dEdx_j = 0; 105 for (CeedInt k = 0; k < 3; k++) { 106 div_u += du[j][k] * dXdx[k][j]; // u_{j,j} = u_{j,K} X_{K,j} 107 dEdx_j += dE[k] * dXdx[k][j]; 108 } 109 u_dot_grad_E += u[j] * dEdx_j; 110 } 111 CeedScalar strong_conv = E * div_u + u_dot_grad_E; 112 113 // Weak Galerkin convection term: dv \cdot (E u) 114 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]); 115 v[4][i] = 0; 116 117 // Strong Galerkin convection term: - v div(E u) 118 v[4][i] = -strong_form * wdetJ * strong_conv; 119 120 // Stabilization requires a measure of element transit time in the velocity 121 // field u. 122 CeedScalar uX[3]; 123 for (CeedInt j = 0; j < 3; j++) uX[j] = dXdx[j][0] * u[0] + dXdx[j][1] * u[1] + dXdx[j][2] * u[2]; 124 const CeedScalar TauS = CtauS / sqrt(Dot3(uX, uX)); 125 for (CeedInt j = 0; j < 3; j++) dv[j][4][i] -= wdetJ * TauS * strong_conv * uX[j]; 126 } // End Quadrature Point Loop 127 128 return 0; 129 } 130 131 CEED_QFUNCTION(IFunction_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 132 IFunction_AdvectionGeneric(ctx, Q, in, out, 3); 133 return 0; 134 } 135 136 // ***************************************************************************** 137 // This QFunction implements consistent outflow and inflow BCs 138 // for 3D advection 139 // 140 // Inflow and outflow faces are determined based on sign(dot(wind, normal)): 141 // sign(dot(wind, normal)) > 0 : outflow BCs 142 // sign(dot(wind, normal)) < 0 : inflow BCs 143 // 144 // Outflow BCs: 145 // The validity of the weak form of the governing equations is extended to the outflow and the current values of E are applied. 146 // 147 // Inflow BCs: 148 // A prescribed Total Energy (E_wind) is applied weakly. 149 // ***************************************************************************** 150 CEED_QFUNCTION(Advection_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 151 // Inputs 152 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 153 const CeedScalar(*q_data_sur) = in[2]; 154 155 // Outputs 156 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 157 AdvectionContext context = (AdvectionContext)ctx; 158 const CeedScalar E_wind = context->E_wind; 159 const CeedScalar strong_form = context->strong_form; 160 const bool is_implicit = context->implicit; 161 162 // Quadrature Point Loop 163 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 164 // Setup 165 // -- Interp in 166 const CeedScalar rho = q[0][i]; 167 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 168 const CeedScalar E = q[4][i]; 169 170 CeedScalar wdetJb, norm[3]; 171 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 172 wdetJb *= is_implicit ? -1. : 1.; 173 174 // Normal velocity 175 const CeedScalar u_normal = norm[0] * u[0] + norm[1] * u[1] + norm[2] * u[2]; 176 177 // No Change in density or momentum 178 for (CeedInt j = 0; j < 4; j++) { 179 v[j][i] = 0; 180 } 181 // Implementing in/outflow BCs 182 if (u_normal > 0) { // outflow 183 v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal; 184 } else { // inflow 185 v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal; 186 } 187 } // End Quadrature Point Loop 188 return 0; 189 } 190 // ***************************************************************************** 191 192 #endif // advection_h 193