1*25125139SJames Wright // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2*25125139SJames Wright // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3*25125139SJames Wright 4*25125139SJames Wright #include <ceed/types.h> 5*25125139SJames Wright #include "newtonian_state.h" 6*25125139SJames Wright 7*25125139SJames Wright CEED_QFUNCTION_HELPER int MonitorTotalKineticEnergy(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, 8*25125139SJames Wright StateVariable state_var) { 9*25125139SJames Wright const NewtonianIdealGasContext gas = (const NewtonianIdealGasContext)ctx; 10*25125139SJames Wright const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 11*25125139SJames Wright const CeedScalar(*Grad_q) = in[1]; 12*25125139SJames Wright const CeedScalar(*q_data) = in[2]; 13*25125139SJames Wright CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 14*25125139SJames Wright 15*25125139SJames Wright CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 16*25125139SJames Wright const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 17*25125139SJames Wright const State s = StateFromQ(gas, qi, state_var); 18*25125139SJames Wright CeedScalar wdetJ, dXdx[3][3], vorticity[3], kmstrain_rate[6], strain_rate[3][3]; 19*25125139SJames Wright State grad_s[3]; 20*25125139SJames Wright 21*25125139SJames Wright QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx); 22*25125139SJames Wright StatePhysicalGradientFromReference(Q, i, gas, s, state_var, Grad_q, dXdx, grad_s); 23*25125139SJames Wright 24*25125139SJames Wright v[0][i] = wdetJ * 0.5 * s.U.density * Dot3(s.Y.velocity, s.Y.velocity); 25*25125139SJames Wright KMStrainRate_State(grad_s, kmstrain_rate); 26*25125139SJames Wright { // See Kundu eq. 4.60 27*25125139SJames Wright CeedScalar div_u = kmstrain_rate[0] + kmstrain_rate[1] + kmstrain_rate[2]; 28*25125139SJames Wright KMUnpack(kmstrain_rate, strain_rate); 29*25125139SJames Wright v[1][i] = wdetJ * -2 * gas->mu * DotN((CeedScalar *)strain_rate, (CeedScalar *)strain_rate, 9); 30*25125139SJames Wright v[2][i] = wdetJ * -gas->lambda * gas->mu * Square(div_u); 31*25125139SJames Wright v[3][i] = wdetJ * s.Y.pressure * div_u; 32*25125139SJames Wright } 33*25125139SJames Wright Vorticity(grad_s, vorticity); 34*25125139SJames Wright v[4][i] = wdetJ * gas->mu * Dot3(vorticity, vorticity); 35*25125139SJames Wright } 36*25125139SJames Wright return 0; 37*25125139SJames Wright } 38*25125139SJames Wright 39*25125139SJames Wright CEED_QFUNCTION(MonitorTotalKineticEnergy_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 40*25125139SJames Wright return MonitorTotalKineticEnergy(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 41*25125139SJames Wright } 42*25125139SJames Wright 43*25125139SJames Wright CEED_QFUNCTION(MonitorTotalKineticEnergy_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 44*25125139SJames Wright return MonitorTotalKineticEnergy(ctx, Q, in, out, STATEVAR_PRIMITIVE); 45*25125139SJames Wright } 46*25125139SJames Wright 47*25125139SJames Wright CEED_QFUNCTION(MonitorTotalKineticEnergy_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 48*25125139SJames Wright return MonitorTotalKineticEnergy(ctx, Q, in, out, STATEVAR_ENTROPY); 49*25125139SJames Wright } 50