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 RHSFunction_AdvectionGeneric(ctx, Q, in, out, 3); 56 return 0; 57 } 58 59 CEED_QFUNCTION(IFunction_Advection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 60 IFunction_AdvectionGeneric(ctx, Q, in, out, 3); 61 return 0; 62 } 63 64 // ***************************************************************************** 65 // This QFunction implements consistent outflow and inflow BCs 66 // for 3D advection 67 // 68 // Inflow and outflow faces are determined based on sign(dot(wind, normal)): 69 // sign(dot(wind, normal)) > 0 : outflow BCs 70 // sign(dot(wind, normal)) < 0 : inflow BCs 71 // 72 // Outflow BCs: 73 // The validity of the weak form of the governing equations is extended to the outflow and the current values of E are applied. 74 // 75 // Inflow BCs: 76 // A prescribed Total Energy (E_wind) is applied weakly. 77 // ***************************************************************************** 78 CEED_QFUNCTION(Advection_InOutFlow)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 79 // Inputs 80 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 81 const CeedScalar(*q_data_sur) = in[2]; 82 83 // Outputs 84 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 85 AdvectionContext context = (AdvectionContext)ctx; 86 const CeedScalar E_wind = context->E_wind; 87 const CeedScalar strong_form = context->strong_form; 88 const bool is_implicit = context->implicit; 89 90 // Quadrature Point Loop 91 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 92 // Setup 93 // -- Interp in 94 const CeedScalar rho = q[0][i]; 95 const CeedScalar u[3] = {q[1][i] / rho, q[2][i] / rho, q[3][i] / rho}; 96 const CeedScalar E = q[4][i]; 97 98 CeedScalar wdetJb, norm[3]; 99 QdataBoundaryUnpack_3D(Q, i, q_data_sur, &wdetJb, NULL, norm); 100 wdetJb *= is_implicit ? -1. : 1.; 101 102 // Normal velocity 103 const CeedScalar u_normal = norm[0] * u[0] + norm[1] * u[1] + norm[2] * u[2]; 104 105 // No Change in density or momentum 106 for (CeedInt j = 0; j < 4; j++) { 107 v[j][i] = 0; 108 } 109 // Implementing in/outflow BCs 110 if (u_normal > 0) { // outflow 111 v[4][i] = -(1 - strong_form) * wdetJb * E * u_normal; 112 } else { // inflow 113 v[4][i] = -(1 - strong_form) * wdetJb * E_wind * u_normal; 114 } 115 } // End Quadrature Point Loop 116 return 0; 117 } 118 // ***************************************************************************** 119 120 #endif // advection_h 121