1 // Copyright (c) 2017-2024, 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 #pragma once 13 14 #include <ceed.h> 15 16 #include "newtonian_state.h" 17 #include "newtonian_types.h" 18 #include "utils.h" 19 #include "utils_eigensolver_jacobi.h" 20 21 // @brief Calculate Frobenius norm of velocity gradient from eigenframe quantities 22 CEED_QFUNCTION_HELPER CeedScalar VelocityGradientMagnitude(const CeedScalar strain_sframe[3], const CeedScalar vorticity_sframe[3]) { 23 return sqrt(Dot3(strain_sframe, strain_sframe) + 0.5 * Dot3(vorticity_sframe, vorticity_sframe)); 24 }; 25 26 // @brief Change the order of basis vectors so that they align with vector and obey right-hand rule 27 // @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 28 // The basis vectors are assumed to form the rows of the basis matrix. 29 CEED_QFUNCTION_HELPER void OrientBasisWithVector(CeedScalar basis[3][3], const CeedScalar vector[3]) { 30 CeedScalar alignment[3] = {0.}, cross[3]; 31 32 MatVec3(basis, vector, CEED_NOTRANSPOSE, alignment); 33 34 if (alignment[0] < 0) ScaleN(basis[0], -1, 3); 35 if (alignment[2] < 0) ScaleN(basis[2], -1, 3); 36 37 Cross3(basis[2], basis[0], cross); 38 CeedScalar basis_1_orientation = Dot3(cross, basis[1]); 39 if (basis_1_orientation < 0) ScaleN(basis[1], -1, 3); 40 } 41 42 // @brief Denormalize outputs using min-max (de-)normalization 43 CEED_QFUNCTION_HELPER void DenormalizeDDOutputs(CeedScalar output[6], const CeedScalar new_bounds[6][2], const CeedScalar old_bounds[6][2]) { 44 CeedScalar bounds_ratio; 45 for (int i = 0; i < 6; i++) { 46 bounds_ratio = (new_bounds[i][1] - new_bounds[i][0]) / (old_bounds[i][1] - old_bounds[i][0]); 47 output[i] = bounds_ratio * (output[i] - old_bounds[i][1]) + new_bounds[i][1]; 48 } 49 } 50 51 /** 52 * @brief Compute model inputs for anisotropic data-driven model 53 * 54 * @param[in] grad_velo_aniso Gradient of velocity in physical (anisotropic) coordinates 55 * @param[in] km_A_ij Anisotropy tensor, in Kelvin-Mandel notation 56 * @param[in] delta Length used to create anisotropy tensor 57 * @param[in] viscosity Kinematic viscosity 58 * @param[out] eigenvectors Eigenvectors of the (anisotropic) velocity gradient 59 * @param[out] inputs Data-driven model inputs 60 * @param[out] grad_velo_magnitude Frobenius norm of the velocity gradient 61 */ 62 CEED_QFUNCTION_HELPER void ComputeSgsDDInputs(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta, 63 const CeedScalar viscosity, CeedScalar eigenvectors[3][3], CeedScalar inputs[6], 64 CeedScalar *grad_velo_magnitude) { 65 CeedScalar strain_sframe[3] = {0.}, vorticity_sframe[3] = {0.}; 66 CeedScalar A_ij[3][3] = {{0.}}, grad_velo_iso[3][3] = {{0.}}; 67 68 // -- Transform physical, anisotropic velocity gradient to isotropic 69 KMUnpack(km_A_ij, A_ij); 70 MatMat3(grad_velo_aniso, A_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, grad_velo_iso); 71 72 { // -- Get Eigenframe 73 CeedScalar kmstrain_iso[6], strain_iso[3][3]; 74 CeedInt work_vector[3] = {0}; 75 KMStrainRate(grad_velo_iso, kmstrain_iso); 76 KMUnpack(kmstrain_iso, strain_iso); 77 Diagonalize3(strain_iso, strain_sframe, eigenvectors, work_vector, SORT_DECREASING_EVALS, true, 5); 78 } 79 80 { // -- Get vorticity in S-frame 81 CeedScalar rotation_iso[3][3]; 82 RotationRate(grad_velo_iso, rotation_iso); 83 CeedScalar vorticity_iso[3] = {-2 * rotation_iso[1][2], 2 * rotation_iso[0][2], -2 * rotation_iso[0][1]}; 84 OrientBasisWithVector(eigenvectors, vorticity_iso); 85 MatVec3(eigenvectors, vorticity_iso, CEED_NOTRANSPOSE, vorticity_sframe); 86 } 87 88 // -- Calculate DD model inputs 89 *grad_velo_magnitude = VelocityGradientMagnitude(strain_sframe, vorticity_sframe); 90 inputs[0] = strain_sframe[0]; 91 inputs[1] = strain_sframe[1]; 92 inputs[2] = strain_sframe[2]; 93 inputs[3] = vorticity_sframe[0]; 94 inputs[4] = vorticity_sframe[1]; 95 inputs[5] = viscosity / Square(delta); 96 ScaleN(inputs, 1 / (*grad_velo_magnitude + CEED_EPSILON), 6); 97 } 98 99 /** 100 * @brief Compute the physical SGS stresses from the neural-network output 101 * 102 * @param[in,out] outputs Outputs from the neural-network 103 * @param[in] delta Length used to create anisotropy tensor 104 * @param[in] eigenvectors Eigenvectors of the (anisotropic) velocity gradient 105 * @param[in] new_bounds Bounds used for min-max de-normalization 106 * @param[in] grad_velo_magnitude Magnitude of the velocity gradient 107 * @param[out] kmsgs_stress Physical SGS stresses in Kelvin-Mandel notation 108 */ 109 CEED_QFUNCTION_HELPER void ComputeSgsDDOutputs(CeedScalar outputs[6], const CeedScalar delta, const CeedScalar eigenvectors[3][3], 110 const CeedScalar new_bounds[6][2], const CeedScalar grad_velo_magnitude, CeedScalar kmsgs_stress[6]) { 111 CeedScalar old_bounds[6][2] = {{0}}; 112 for (int j = 0; j < 6; j++) old_bounds[j][1] = 1; 113 DenormalizeDDOutputs(outputs, new_bounds, old_bounds); 114 115 // Re-dimensionalize sgs_stress 116 ScaleN(outputs, Square(delta) * Square(grad_velo_magnitude), 6); 117 118 CeedScalar sgs_stress[3][3] = {{0.}}; 119 { // Rotate SGS Stress back to physical frame, SGS_physical = E^T SGS_sframe E 120 CeedScalar Evec_sgs[3][3] = {{0.}}; 121 const CeedScalar sgs_sframe[3][3] = { 122 {outputs[0], outputs[3], outputs[4]}, 123 {outputs[3], outputs[1], outputs[5]}, 124 {outputs[4], outputs[5], outputs[2]}, 125 }; 126 MatMat3(eigenvectors, sgs_sframe, CEED_TRANSPOSE, CEED_NOTRANSPOSE, Evec_sgs); 127 MatMat3(Evec_sgs, eigenvectors, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, sgs_stress); 128 } 129 130 KMPack(sgs_stress, kmsgs_stress); 131 } 132