1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 #include <ceed.h> 5 6 #include "newtonian_state.h" 7 #include "utils.h" 8 9 // @brief Volume integral for RHS of divergence of diffusive flux projection 10 CEED_QFUNCTION_HELPER int DivDiffusiveFluxVolumeRHS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 11 StateVariable state_var) { 12 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 13 const CeedScalar(*Grad_q) = in[1]; 14 const CeedScalar(*q_data) = in[2]; 15 CeedScalar(*Grad_v)[4][CEED_Q_VLA] = (CeedScalar(*)[4][CEED_Q_VLA])out[0]; 16 17 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 18 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 19 20 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 21 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 22 const State s = StateFromQ(context, qi, state_var); 23 CeedScalar wdetJ, dXdx[3][3]; 24 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 25 26 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 27 { // Get stress and Fe 28 State grad_s[3]; 29 CeedScalar strain_rate[6], kmstress[6]; 30 31 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 32 KMStrainRate_State(grad_s, strain_rate); 33 NewtonianStress(context, strain_rate, kmstress); 34 KMUnpack(kmstress, stress); 35 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 36 } 37 38 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 39 40 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 41 for (CeedInt k = 0; k < 3; k++) { 42 Grad_v[k][j - 1][i] = -wdetJ * Dot3(dXdx[k], Fdiff[j]); 43 } 44 } 45 } 46 return 0; 47 } 48 49 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 50 return DivDiffusiveFluxVolumeRHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 51 } 52 53 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 54 return DivDiffusiveFluxVolumeRHS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 55 } 56 57 CEED_QFUNCTION(DivDiffusiveFluxVolumeRHS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 58 return DivDiffusiveFluxVolumeRHS(ctx, Q, in, out, STATEVAR_ENTROPY); 59 } 60 61 // @brief Boundary integral for RHS of divergence of diffusive flux projection 62 CEED_QFUNCTION_HELPER int DivDiffusiveFluxBoundaryRHS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 63 StateVariable state_var) { 64 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 65 const CeedScalar(*Grad_q) = in[1]; 66 const CeedScalar(*q_data) = in[2]; 67 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 68 69 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 70 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 71 72 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 73 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 74 const State s = StateFromQ(context, qi, state_var); 75 CeedScalar wdetJ, dXdx[3][3], normal[3]; 76 CeedScalar stress[3][3], Fe[3], Fdiff[5]; 77 78 QdataBoundaryGradientUnpack_3D(Q, i, q_data, &wdetJ, dXdx, normal); 79 { // Get stress and Fe 80 State grad_s[3]; 81 CeedScalar strain_rate[6], kmstress[6]; 82 83 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 84 KMStrainRate_State(grad_s, strain_rate); 85 NewtonianStress(context, strain_rate, kmstress); 86 KMUnpack(kmstress, stress); 87 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 88 } 89 90 FluxTotal_Boundary(ZeroInviscidFluxes, stress, Fe, normal, Fdiff); 91 92 // Continuity has no diffusive flux, therefore skip 93 for (CeedInt j = 1; j < 5; j++) v[j - 1][i] = wdetJ * Fdiff[j]; 94 } 95 return 0; 96 } 97 98 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 99 return DivDiffusiveFluxBoundaryRHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 100 } 101 102 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 103 return DivDiffusiveFluxBoundaryRHS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 104 } 105 106 CEED_QFUNCTION(DivDiffusiveFluxBoundaryRHS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 107 return DivDiffusiveFluxBoundaryRHS(ctx, Q, in, out, STATEVAR_ENTROPY); 108 } 109 110 // @brief Integral for RHS of diffusive flux projection 111 CEED_QFUNCTION_HELPER int DiffusiveFluxRHS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 112 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 113 const CeedScalar(*Grad_q) = in[1]; 114 const CeedScalar(*q_data) = in[2]; 115 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 116 117 const NewtonianIdealGasContext context = (NewtonianIdealGasContext)ctx; 118 const StateConservative ZeroInviscidFluxes[3] = {{0}}; 119 120 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 121 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 122 const State s = StateFromQ(context, qi, state_var); 123 CeedScalar wdetJ, dXdx[3][3]; 124 CeedScalar stress[3][3], Fe[3], Fdiff[5][3]; 125 126 QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 127 { // Get stress and Fe 128 State grad_s[3]; 129 CeedScalar strain_rate[6], kmstress[6]; 130 131 StatePhysicalGradientFromReference(Q, i, context, s, state_var, Grad_q, dXdx, grad_s); 132 KMStrainRate_State(grad_s, strain_rate); 133 NewtonianStress(context, strain_rate, kmstress); 134 KMUnpack(kmstress, stress); 135 ViscousEnergyFlux(context, s.Y, grad_s, stress, Fe); 136 } 137 138 FluxTotal(ZeroInviscidFluxes, stress, Fe, Fdiff); 139 140 for (CeedInt j = 1; j < 5; j++) { // Continuity has no diffusive flux, therefore skip 141 for (CeedInt k = 0; k < 3; k++) { 142 v[(j - 1) * 3 + k][i] = wdetJ * Fdiff[j][k]; 143 } 144 } 145 } 146 return 0; 147 } 148 149 CEED_QFUNCTION(DiffusiveFluxRHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 150 return DiffusiveFluxRHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 151 } 152 153 CEED_QFUNCTION(DiffusiveFluxRHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 154 return DiffusiveFluxRHS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 155 } 156 157 CEED_QFUNCTION(DiffusiveFluxRHS_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 158 return DiffusiveFluxRHS(ctx, Q, in, out, STATEVAR_ENTROPY); 159 } 160 161 // @brief QFunction to calculate the divergence of the diffusive flux 162 CEED_QFUNCTION_HELPER int ComputeDivDiffusiveFluxGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, const CeedInt dim, 163 const CeedInt num_comps) { 164 const CeedScalar *grad_q = in[0]; 165 const CeedScalar(*q_data) = in[1]; 166 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 167 168 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 169 CeedScalar dXdx[9]; 170 171 QdataUnpack_ND(dim, Q, i, q_data, NULL, dXdx); 172 CeedPragmaSIMD for (CeedInt n = 0; n < num_comps; n++) { 173 CeedScalar grad_qn[9]; 174 175 // Get gradient into dim x dim matrix form, with orientation [flux_direction][gradient_direction] 176 // Equivalent of GradUnpackN 177 const CeedInt offset = Q * n * dim; // offset to reach nth component flux gradients 178 for (CeedInt g = 0; g < dim; g++) { 179 for (CeedInt f = 0; f < dim; f++) { 180 grad_qn[f * dim + g] = grad_q[offset + (Q * num_comps * dim) * g + Q * f + i]; 181 } 182 } 183 v[n][i] = 0; 184 DivergenceND(grad_qn, dXdx, dim, &v[n][i]); 185 } 186 } 187 return 0; 188 } 189 190 CEED_QFUNCTION(ComputeDivDiffusiveFlux3D_4)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 191 return ComputeDivDiffusiveFluxGeneric(ctx, Q, in, out, 3, 4); 192 } 193