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