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