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
6 #include "utils.h"
7
8 // @brief QFunction to calculate the divergence of the diffusive flux
ComputeDivDiffusiveFluxGeneric(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out,const CeedInt dim,const CeedInt num_comps)9 CEED_QFUNCTION_HELPER int ComputeDivDiffusiveFluxGeneric(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, const CeedInt dim,
10 const CeedInt num_comps) {
11 const CeedScalar *grad_q = in[0];
12 const CeedScalar(*q_data) = in[1];
13 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
14
15 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
16 CeedScalar dXdx[9];
17
18 QdataUnpack_ND(dim, Q, i, q_data, NULL, dXdx);
19 CeedPragmaSIMD for (CeedInt n = 0; n < num_comps; n++) {
20 CeedScalar grad_qn[9];
21
22 // Get gradient into dim x dim matrix form, with orientation [flux_direction][gradient_direction]
23 // Equivalent of GradUnpackN
24 const CeedInt offset = Q * n * dim; // offset to reach nth component flux gradients
25 for (CeedInt g = 0; g < dim; g++) {
26 for (CeedInt f = 0; f < dim; f++) {
27 grad_qn[f * dim + g] = grad_q[offset + (Q * num_comps * dim) * g + Q * f + i];
28 }
29 }
30 v[n][i] = 0;
31 DivergenceND(grad_qn, dXdx, dim, &v[n][i]);
32 }
33 }
34 return 0;
35 }
36
ComputeDivDiffusiveFlux3D_4(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)37 CEED_QFUNCTION(ComputeDivDiffusiveFlux3D_4)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
38 return ComputeDivDiffusiveFluxGeneric(ctx, Q, in, out, 3, 4);
39 }
40
ComputeDivDiffusiveFlux3D_1(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)41 CEED_QFUNCTION(ComputeDivDiffusiveFlux3D_1)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
42 return ComputeDivDiffusiveFluxGeneric(ctx, Q, in, out, 3, 1);
43 }
44
ComputeDivDiffusiveFlux2D_1(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)45 CEED_QFUNCTION(ComputeDivDiffusiveFlux2D_1)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
46 return ComputeDivDiffusiveFluxGeneric(ctx, Q, in, out, 2, 1);
47 }
48