1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors.
2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause
3
4 /// @file
5 /// Element anisotropy tensor, as defined in 'Invariant data-driven subgrid stress modeling in the strain-rate eigenframe for large eddy simulation'
6 /// Prakash et al. 2022
7 #include <ceed/types.h>
8
9 #include "utils.h"
10 #include "utils_eigensolver_jacobi.h"
11
12 // @brief Get Anisotropy tensor from xi_{i,j}
13 // @details A_ij = \Delta_{ij} / ||\Delta_ij||, \Delta_ij = (xi_{i,j})^(-1/2)
AnisotropyTensor(const CeedScalar km_g_ij[6],CeedScalar A_ij[3][3],CeedScalar * delta,const CeedInt n_sweeps)14 CEED_QFUNCTION_HELPER void AnisotropyTensor(const CeedScalar km_g_ij[6], CeedScalar A_ij[3][3], CeedScalar *delta, const CeedInt n_sweeps) {
15 CeedScalar evals[3], evecs[3][3], evals_evecs[3][3] = {{0.}}, g_ij[3][3];
16 CeedInt work_vector[3];
17
18 // Invert square root of metric tensor to get \Delta_ij
19 KMUnpack(km_g_ij, g_ij);
20 Diagonalize3(g_ij, evals, evecs, work_vector, SORT_DECREASING_EVALS, true, n_sweeps);
21 for (int i = 0; i < 3; i++) evals[i] = 1 / sqrt(evals[i]);
22 MatDiag3(evecs, evals, CEED_NOTRANSPOSE, evals_evecs);
23 MatMat3(evecs, evals_evecs, CEED_TRANSPOSE, CEED_NOTRANSPOSE, A_ij); // A_ij = E^T D E
24
25 // Scale by delta to get anisotropy tensor
26 *delta = Norm3(evals);
27 ScaleN((CeedScalar *)A_ij, 1 / *delta, 9);
28 // NOTE Need 2 factor to get physical element size (rather than projected onto [-1,1]^dim)
29 // Should attempt to auto-determine this from the quadrature point coordinates in reference space
30 *delta *= 2;
31 }
32
33 // @brief RHS for L^2 projection of anisotropic tensor and it's Frobenius norm
AnisotropyTensorProjection(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)34 CEED_QFUNCTION(AnisotropyTensorProjection)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
35 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
36 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
37
38 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
39 const CeedScalar wdetJ = q_data[0][i];
40 const CeedScalar dXdx[3][3] = {
41 {q_data[1][i], q_data[2][i], q_data[3][i]},
42 {q_data[4][i], q_data[5][i], q_data[6][i]},
43 {q_data[7][i], q_data[8][i], q_data[9][i]}
44 };
45
46 CeedScalar km_g_ij[6] = {0.}, A_ij[3][3] = {{0.}}, km_A_ij[6], delta;
47 KMMetricTensor(dXdx, km_g_ij);
48 AnisotropyTensor(km_g_ij, A_ij, &delta, 15);
49 KMPack(A_ij, km_A_ij);
50
51 for (CeedInt j = 0; j < 6; j++) v[j][i] = wdetJ * km_A_ij[j];
52 v[6][i] = wdetJ * delta;
53 }
54 return 0;
55 }
56
57 // @brief Get anisotropic tensor and it's Frobenius norm at quadrature points
AnisotropyTensorCollocate(void * ctx,CeedInt Q,const CeedScalar * const * in,CeedScalar * const * out)58 CEED_QFUNCTION(AnisotropyTensorCollocate)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
59 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
60 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
61
62 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) {
63 const CeedScalar dXdx[3][3] = {
64 {q_data[1][i], q_data[2][i], q_data[3][i]},
65 {q_data[4][i], q_data[5][i], q_data[6][i]},
66 {q_data[7][i], q_data[8][i], q_data[9][i]}
67 };
68
69 CeedScalar km_g_ij[6] = {0.}, A_ij[3][3] = {{0.}}, km_A_ij[6], delta;
70 KMMetricTensor(dXdx, km_g_ij);
71 AnisotropyTensor(km_g_ij, A_ij, &delta, 15);
72 KMPack(A_ij, km_A_ij);
73
74 for (CeedInt j = 0; j < 6; j++) v[j][i] = km_A_ij[j];
75 v[6][i] = delta;
76 }
77 return 0;
78 }
79