xref: /honee/qfunctions/monitor_cfl.h (revision 87fd7f33e5ff541d346273b55e9eb051d4de0b43)
1*87fd7f33SJames Wright // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2*87fd7f33SJames Wright // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3*87fd7f33SJames Wright 
4*87fd7f33SJames Wright #include <ceed/types.h>
5*87fd7f33SJames Wright #include "newtonian_state.h"
6*87fd7f33SJames Wright 
7*87fd7f33SJames Wright CEED_QFUNCTION_HELPER int MonitorCFL(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var,
8*87fd7f33SJames Wright                                      CeedInt dim) {
9*87fd7f33SJames Wright   const NewtonianIdealGasContext gas = (const NewtonianIdealGasContext)ctx;
10*87fd7f33SJames Wright   const CeedScalar(*q)[CEED_Q_VLA]   = (const CeedScalar(*)[CEED_Q_VLA])in[0];
11*87fd7f33SJames Wright   const CeedScalar(*q_data)          = in[1];
12*87fd7f33SJames Wright   CeedScalar(*v)                     = out[0];
13*87fd7f33SJames Wright 
14*87fd7f33SJames Wright   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
15*87fd7f33SJames Wright     const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
16*87fd7f33SJames Wright     const State      s     = StateFromQ(gas, qi, state_var);
17*87fd7f33SJames Wright     CeedScalar       wdetJ, dXdx[9], gij_uj[3] = {0.}, gijd_mat[9] = {0.};
18*87fd7f33SJames Wright 
19*87fd7f33SJames Wright     QdataUnpack_ND(dim, Q, i, q_data, &wdetJ, dXdx);
20*87fd7f33SJames Wright     MatMatN(dXdx, dXdx, dim, CEED_TRANSPOSE, CEED_NOTRANSPOSE, gijd_mat);
21*87fd7f33SJames Wright     // (1/2)^2 to account for reference element size; for length 1 square/cube element, gij should be identity matrix
22*87fd7f33SJames Wright     ScaleN((CeedScalar *)gijd_mat, 0.25, Square(dim));
23*87fd7f33SJames Wright     // u_i g_ij u_j
24*87fd7f33SJames Wright     MatVecNM(gijd_mat, s.Y.velocity, dim, dim, CEED_NOTRANSPOSE, gij_uj);
25*87fd7f33SJames Wright     v[i] = sqrt(Dot3(s.Y.velocity, gij_uj));
26*87fd7f33SJames Wright   }
27*87fd7f33SJames Wright   return 0;
28*87fd7f33SJames Wright }
29*87fd7f33SJames Wright 
30*87fd7f33SJames Wright CEED_QFUNCTION(MonitorCFL_3D_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
31*87fd7f33SJames Wright   return MonitorCFL(ctx, Q, in, out, STATEVAR_CONSERVATIVE, 3);
32*87fd7f33SJames Wright }
33*87fd7f33SJames Wright 
34*87fd7f33SJames Wright CEED_QFUNCTION(MonitorCFL_3D_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
35*87fd7f33SJames Wright   return MonitorCFL(ctx, Q, in, out, STATEVAR_PRIMITIVE, 3);
36*87fd7f33SJames Wright }
37*87fd7f33SJames Wright 
38*87fd7f33SJames Wright CEED_QFUNCTION(MonitorCFL_3D_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
39*87fd7f33SJames Wright   return MonitorCFL(ctx, Q, in, out, STATEVAR_ENTROPY, 3);
40*87fd7f33SJames Wright }
41*87fd7f33SJames Wright 
42*87fd7f33SJames Wright CEED_QFUNCTION(MonitorCFL_2D_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
43*87fd7f33SJames Wright   return MonitorCFL(ctx, Q, in, out, STATEVAR_CONSERVATIVE, 2);
44*87fd7f33SJames Wright }
45*87fd7f33SJames Wright 
46*87fd7f33SJames Wright CEED_QFUNCTION(MonitorCFL_2D_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
47*87fd7f33SJames Wright   return MonitorCFL(ctx, Q, in, out, STATEVAR_PRIMITIVE, 2);
48*87fd7f33SJames Wright }
49*87fd7f33SJames Wright 
50*87fd7f33SJames Wright CEED_QFUNCTION(MonitorCFL_2D_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
51*87fd7f33SJames Wright   return MonitorCFL(ctx, Q, in, out, STATEVAR_ENTROPY, 2);
52*87fd7f33SJames Wright }
53