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 /// Structs and helper functions for data-driven subgrid-stress modeling 10 /// See 'Invariant data-driven subgrid stress modeling in the strain-rate eigenframe for large eddy simulation' 2022 and 'S-frame discrepancy 11 /// correction models for data-informed Reynolds stress closure' 2022 12 13 #ifndef sgs_dd_utils_h 14 #define sgs_dd_utils_h 15 16 #include <ceed.h> 17 18 #include "newtonian_state.h" 19 #include "newtonian_types.h" 20 #include "utils.h" 21 #include "utils_eigensolver_jacobi.h" 22 23 // @brief Calculate the inverse of the multiplicity, reducing to a single component 24 CEED_QFUNCTION(InverseMultiplicity)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 25 const CeedScalar(*multiplicity)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 26 CeedScalar(*inv_multiplicity) = (CeedScalar(*))out[0]; 27 28 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) inv_multiplicity[i] = 1.0 / multiplicity[0][i]; 29 return 0; 30 } 31 32 // @brief Calculate Frobenius norm of velocity gradient from eigenframe quantities 33 CEED_QFUNCTION_HELPER CeedScalar VelocityGradientMagnitude(const CeedScalar strain_sframe[3], const CeedScalar vorticity_sframe[3]) { 34 return sqrt(Dot3(strain_sframe, strain_sframe) + 0.5 * Dot3(vorticity_sframe, vorticity_sframe)); 35 }; 36 37 // @brief Change the order of basis vectors so that they align with vector and obey right-hand rule 38 // @details The e_1 and e_3 basis vectors are the closest aligned to the vector. The e_2 is set via e_3 x e_1 39 // The basis vectors are assumed to form the rows of the basis matrix. 40 CEED_QFUNCTION_HELPER void OrientBasisWithVector(CeedScalar basis[3][3], const CeedScalar vector[3]) { 41 CeedScalar alignment[3] = {0.}, cross[3]; 42 43 MatVec3(basis, vector, CEED_NOTRANSPOSE, alignment); 44 45 if (alignment[0] < 0) ScaleN(basis[0], -1, 3); 46 if (alignment[2] < 0) ScaleN(basis[2], -1, 3); 47 48 Cross3(basis[2], basis[0], cross); 49 CeedScalar basis_1_orientation = Dot3(cross, basis[1]); 50 if (basis_1_orientation < 0) ScaleN(basis[1], -1, 3); 51 } 52 53 /** 54 * @brief Compute model inputs for anisotropic data-driven model 55 * 56 * @param[in] grad_velo_aniso Gradient of velocity in physical (anisotropic) coordinates 57 * @param[in] km_A_ij Anisotropy tensor, in Kelvin-Mandel notation 58 * @param[in] delta Length used to create anisotropy tensor 59 * @param[in] viscosity Kinematic viscosity 60 * @param[out] eigenvectors Eigenvectors of the (anisotropic) velocity gradient 61 * @param[out] inputs Data-driven model inputs 62 * @param[out] grad_velo_magnitude Frobenius norm of the velocity gradient 63 */ 64 CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropicInputs(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta, 65 const CeedScalar viscosity, CeedScalar eigenvectors[3][3], CeedScalar inputs[6], 66 CeedScalar *grad_velo_magnitude) { 67 CeedScalar strain_sframe[3] = {0.}, vorticity_sframe[3] = {0.}; 68 CeedScalar A_ij[3][3] = {{0.}}, grad_velo_iso[3][3] = {{0.}}; 69 70 // -- Transform physical, anisotropic velocity gradient to isotropic 71 KMUnpack(km_A_ij, A_ij); 72 MatMat3(grad_velo_aniso, A_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, grad_velo_iso); 73 74 { // -- Get Eigenframe 75 CeedScalar kmstrain_iso[6], strain_iso[3][3]; 76 CeedInt work_vector[3] = {0}; 77 KMStrainRate(grad_velo_iso, kmstrain_iso); 78 KMUnpack(kmstrain_iso, strain_iso); 79 Diagonalize3(strain_iso, strain_sframe, eigenvectors, work_vector, SORT_DECREASING_EVALS, true, 5); 80 } 81 82 { // -- Get vorticity in S-frame 83 CeedScalar rotation_iso[3][3]; 84 RotationRate(grad_velo_iso, rotation_iso); 85 CeedScalar vorticity_iso[3] = {-2 * rotation_iso[1][2], 2 * rotation_iso[0][2], -2 * rotation_iso[0][1]}; 86 OrientBasisWithVector(eigenvectors, vorticity_iso); 87 MatVec3(eigenvectors, vorticity_iso, CEED_NOTRANSPOSE, vorticity_sframe); 88 } 89 90 // -- Calculate DD model inputs 91 *grad_velo_magnitude = VelocityGradientMagnitude(strain_sframe, vorticity_sframe); 92 inputs[0] = strain_sframe[0]; 93 inputs[1] = strain_sframe[1]; 94 inputs[2] = strain_sframe[2]; 95 inputs[3] = vorticity_sframe[0]; 96 inputs[4] = vorticity_sframe[1]; 97 inputs[5] = viscosity / Square(delta); 98 ScaleN(inputs, 1 / (*grad_velo_magnitude + CEED_EPSILON), 6); 99 } 100 101 #endif // sgs_dd_utils_h 102