xref: /honee/qfunctions/monitor_cfl.h (revision c9f376050d3861e44a3522dda28dab81f9a3076c)
1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3 
4 #include <ceed/types.h>
5 #include "newtonian_state.h"
6 #include "numerics.h"
7 
8 CEED_QFUNCTION_HELPER int MonitorCFL(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var,
9                                      CeedInt dim) {
10   const NewtonianIdealGasContext gas = (const NewtonianIdealGasContext)ctx;
11   const CeedScalar(*q)[CEED_Q_VLA]   = (const CeedScalar(*)[CEED_Q_VLA])in[0];
12   const CeedScalar(*q_data)          = in[1];
13   CeedScalar(*v)                     = out[0];
14 
15   CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
16     const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]};
17     const State      s     = StateFromQ(gas, qi, state_var);
18 
19     switch (dim) {
20       case 2: {
21         CeedScalar wdetJ, dXdx[2][2], gijd_mat[2][2] = {{0.}};
22 
23         QdataUnpack_2D(Q, i, q_data, &wdetJ, dXdx);
24         MatMat2(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, gijd_mat);
25         // (1/2)^2 to account for reference element size; for length 1 square/cube element, gij should be identity matrix
26         ScaleN((CeedScalar *)gijd_mat, 0.25, Square(dim));
27         // Multiplied by timestep outside of QFunction so that it can be used for both Advection and Newtonian
28         v[i] = CalculateCFL_2D(s.Y.velocity, 1, gijd_mat);
29       } break;
30       case 3: {
31         CeedScalar wdetJ, dXdx[3][3], gijd_mat[3][3] = {{0.}};
32 
33         QdataUnpack_3D(Q, i, q_data, &wdetJ, dXdx);
34         MatMat3(dXdx, dXdx, CEED_TRANSPOSE, CEED_NOTRANSPOSE, gijd_mat);
35         // (1/2)^2 to account for reference element size; for length 1 square/cube element, gij should be identity matrix
36         ScaleN((CeedScalar *)gijd_mat, 0.25, Square(dim));
37         // Multiplied by timestep outside of QFunction so that it can be used for both Advection and Newtonian
38         v[i] = CalculateCFL_3D(s.Y.velocity, 1, gijd_mat);
39       } break;
40     }
41   }
42   return 0;
43 }
44 
45 CEED_QFUNCTION(MonitorCFL_3D_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
46   return MonitorCFL(ctx, Q, in, out, STATEVAR_CONSERVATIVE, 3);
47 }
48 
49 CEED_QFUNCTION(MonitorCFL_3D_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
50   return MonitorCFL(ctx, Q, in, out, STATEVAR_PRIMITIVE, 3);
51 }
52 
53 CEED_QFUNCTION(MonitorCFL_3D_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
54   return MonitorCFL(ctx, Q, in, out, STATEVAR_ENTROPY, 3);
55 }
56 
57 CEED_QFUNCTION(MonitorCFL_2D_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
58   return MonitorCFL(ctx, Q, in, out, STATEVAR_CONSERVATIVE, 2);
59 }
60 
61 CEED_QFUNCTION(MonitorCFL_2D_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
62   return MonitorCFL(ctx, Q, in, out, STATEVAR_PRIMITIVE, 2);
63 }
64 
65 CEED_QFUNCTION(MonitorCFL_2D_Entropy)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
66   return MonitorCFL(ctx, Q, in, out, STATEVAR_ENTROPY, 2);
67 }
68