Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fluids: Decompose SGS DD input and output processing #1223

Merged
merged 3 commits into from
Aug 24, 2023
Merged

Fluids: Decompose SGS DD input and output processing #1223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
7dfc0c3
fluids: Extract data-driven model input calculation
jrwrigh Apr 16, 2023
1abde36
fluids: Move some DD SGS functions to sgs_dd_utils.h
jrwrigh Apr 16, 2023
b32ff2a
fluids: Extract out the DD output processing
jrwrigh May 25, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
98 changes: 6 additions & 92 deletions examples/fluids/qfunctions/sgs_dd_model.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
// This file is part of CEED: http://github.com/ceed

/// @file
/// Structs and helper functions for data-driven subgrid-stress modeling
/// Structs and helper functions to evaluate data-driven subgrid-stress modeling
/// See 'Invariant data-driven subgrid stress modeling in the strain-rate eigenframe for large eddy simulation' 2022 and 'S-frame discrepancy
/// correction models for data-informed Reynolds stress closure' 2022

Expand All @@ -17,6 +17,7 @@

#include "newtonian_state.h"
#include "newtonian_types.h"
#include "sgs_dd_utils.h"
#include "utils.h"
#include "utils_eigensolver_jacobi.h"

Expand All @@ -37,45 +38,6 @@ struct SGS_DD_ModelContext_ {
CeedScalar data[1];
};

// @brief Calculate the inverse of the multiplicity, reducing to a single component
CEED_QFUNCTION(InverseMultiplicity)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
const CeedScalar(*multiplicity)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
CeedScalar(*inv_multiplicity) = (CeedScalar(*))out[0];

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) inv_multiplicity[i] = 1.0 / multiplicity[0][i];
return 0;
}

// @brief Calculate Frobenius norm of velocity gradient from eigenframe quantities
CEED_QFUNCTION_HELPER CeedScalar VelocityGradientMagnitude(const CeedScalar strain_sframe[3], const CeedScalar vorticity_sframe[3]) {
return sqrt(Dot3(strain_sframe, strain_sframe) + 0.5 * Dot3(vorticity_sframe, vorticity_sframe));
};

// @brief Denormalize outputs using min-max (de-)normalization
CEED_QFUNCTION_HELPER void DenormalizeDDOutputs(CeedScalar output[6], const CeedScalar (*new_bounds)[2], const CeedScalar old_bounds[6][2]) {
CeedScalar bounds_ratio;
for (int i = 0; i < 6; i++) {
bounds_ratio = (new_bounds[i][1] - new_bounds[i][0]) / (old_bounds[i][1] - old_bounds[i][0]);
output[i] = bounds_ratio * (output[i] - old_bounds[i][1]) + new_bounds[i][1];
}
}

// @brief Change the order of basis vectors so that they align with vector and obey right-hand rule
// @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
// The basis vectors are assumed to form the rows of the basis matrix.
CEED_QFUNCTION_HELPER void OrientBasisWithVector(CeedScalar basis[3][3], const CeedScalar vector[3]) {
CeedScalar alignment[3] = {0.}, cross[3];

MatVec3(basis, vector, CEED_NOTRANSPOSE, alignment);

if (alignment[0] < 0) ScaleN(basis[0], -1, 3);
if (alignment[2] < 0) ScaleN(basis[2], -1, 3);

Cross3(basis[2], basis[0], cross);
CeedScalar basis_1_orientation = Dot3(cross, basis[1]);
if (basis_1_orientation < 0) ScaleN(basis[1], -1, 3);
}

CEED_QFUNCTION_HELPER void LeakyReLU(CeedScalar *x, const CeedScalar alpha, const CeedInt N) {
for (CeedInt i = 0; i < N; i++) x[i] *= (x[i] < 0 ? alpha : 1.);
}
Expand All @@ -100,60 +62,12 @@ CEED_QFUNCTION_HELPER void DataDrivenInference(const CeedScalar *inputs, CeedSca

CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropic(const CeedScalar grad_velo_aniso[3][3], const CeedScalar km_A_ij[6], const CeedScalar delta,
const CeedScalar viscosity, CeedScalar kmsgs_stress[6], SGS_DDModelContext sgsdd_ctx) {
CeedScalar strain_sframe[3] = {0.}, vorticity_sframe[3] = {0.}, eigenvectors[3][3];
CeedScalar A_ij[3][3] = {{0.}}, grad_velo_iso[3][3] = {{0.}};

// -- Unpack anisotropy tensor
KMUnpack(km_A_ij, A_ij);

// -- Transform physical, anisotropic velocity gradient to isotropic
MatMat3(grad_velo_aniso, A_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, grad_velo_iso);

{ // -- Get Eigenframe
CeedScalar kmstrain_iso[6], strain_iso[3][3];
CeedInt work_vector[3] = {0};
KMStrainRate(grad_velo_iso, kmstrain_iso);
KMUnpack(kmstrain_iso, strain_iso);
Diagonalize3(strain_iso, strain_sframe, eigenvectors, work_vector, SORT_DECREASING_EVALS, true, 5);
}

{ // -- Get vorticity in S-frame
CeedScalar rotation_iso[3][3];
RotationRate(grad_velo_iso, rotation_iso);
CeedScalar vorticity_iso[3] = {-2 * rotation_iso[1][2], 2 * rotation_iso[0][2], -2 * rotation_iso[0][1]};
OrientBasisWithVector(eigenvectors, vorticity_iso);
MatVec3(eigenvectors, vorticity_iso, CEED_NOTRANSPOSE, vorticity_sframe);
}

// -- Setup DD model inputs
const CeedScalar grad_velo_magnitude = VelocityGradientMagnitude(strain_sframe, vorticity_sframe);
CeedScalar inputs[6] = {strain_sframe[0], strain_sframe[1], strain_sframe[2], vorticity_sframe[0], vorticity_sframe[1], viscosity / Square(delta)};
ScaleN(inputs, 1 / (grad_velo_magnitude + CEED_EPSILON), 6);

CeedScalar sgs_sframe_sym[6] = {0.};
DataDrivenInference(inputs, sgs_sframe_sym, sgsdd_ctx);

CeedScalar old_bounds[6][2] = {{0}};
for (int j = 0; j < 6; j++) old_bounds[j][1] = 1;
CeedScalar inputs[6], grad_velo_magnitude, eigenvectors[3][3], sgs_sframe_sym[6] = {0.};
const CeedScalar(*new_bounds)[2] = (const CeedScalar(*)[2]) & sgsdd_ctx->data[sgsdd_ctx->offsets.out_scaling];
DenormalizeDDOutputs(sgs_sframe_sym, new_bounds, old_bounds);

// Re-dimensionalize sgs_stress
ScaleN(sgs_sframe_sym, Square(delta) * Square(grad_velo_magnitude), 6);

CeedScalar sgs_stress[3][3] = {{0.}};
{ // Rotate SGS Stress back to physical frame, SGS_physical = E^T SGS_sframe E
CeedScalar Evec_sgs[3][3] = {{0.}};
const CeedScalar sgs_sframe[3][3] = {
{sgs_sframe_sym[0], sgs_sframe_sym[3], sgs_sframe_sym[4]},
{sgs_sframe_sym[3], sgs_sframe_sym[1], sgs_sframe_sym[5]},
{sgs_sframe_sym[4], sgs_sframe_sym[5], sgs_sframe_sym[2]},
};
MatMat3(eigenvectors, sgs_sframe, CEED_TRANSPOSE, CEED_NOTRANSPOSE, Evec_sgs);
MatMat3(Evec_sgs, eigenvectors, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, sgs_stress);
}

KMPack(sgs_stress, kmsgs_stress);
ComputeSGS_DDAnisotropicInputs(grad_velo_aniso, km_A_ij, delta, viscosity, eigenvectors, inputs, &grad_velo_magnitude);
DataDrivenInference(inputs, sgs_sframe_sym, sgsdd_ctx);
ComputeSGS_DDAnisotropicOutputs(sgs_sframe_sym, delta, eigenvectors, new_bounds, grad_velo_magnitude, kmsgs_stress);
}

// @brief Calculate subgrid stress at nodes using anisotropic data-driven model
Expand Down
145 changes: 145 additions & 0 deletions examples/fluids/qfunctions/sgs_dd_utils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,145 @@
// Copyright (c) 2017-2023, Lawrence Livermore National Security, LLC and other CEED contributors.
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
//
// SPDX-License-Identifier: BSD-2-Clause
//
// This file is part of CEED: http://github.com/ceed

/// @file
/// Structs and helper functions for data-driven subgrid-stress modeling
/// See 'Invariant data-driven subgrid stress modeling in the strain-rate eigenframe for large eddy simulation' 2022 and 'S-frame discrepancy
/// correction models for data-informed Reynolds stress closure' 2022

#ifndef sgs_dd_utils_h
#define sgs_dd_utils_h

#include <ceed.h>

#include "newtonian_state.h"
#include "newtonian_types.h"
#include "utils.h"
#include "utils_eigensolver_jacobi.h"

// @brief Calculate the inverse of the multiplicity, reducing to a single component
CEED_QFUNCTION(InverseMultiplicity)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) {
const CeedScalar(*multiplicity)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0];
CeedScalar(*inv_multiplicity) = (CeedScalar(*))out[0];

CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) inv_multiplicity[i] = 1.0 / multiplicity[0][i];
return 0;
}
jrwrigh marked this conversation as resolved.
Show resolved Hide resolved

// @brief Calculate Frobenius norm of velocity gradient from eigenframe quantities
CEED_QFUNCTION_HELPER CeedScalar VelocityGradientMagnitude(const CeedScalar strain_sframe[3], const CeedScalar vorticity_sframe[3]) {
return sqrt(Dot3(strain_sframe, strain_sframe) + 0.5 * Dot3(vorticity_sframe, vorticity_sframe));
};

// @brief Change the order of basis vectors so that they align with vector and obey right-hand rule
// @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
// The basis vectors are assumed to form the rows of the basis matrix.
CEED_QFUNCTION_HELPER void OrientBasisWithVector(CeedScalar basis[3][3], const CeedScalar vector[3]) {
CeedScalar alignment[3] = {0.}, cross[3];

MatVec3(basis, vector, CEED_NOTRANSPOSE, alignment);

if (alignment[0] < 0) ScaleN(basis[0], -1, 3);
if (alignment[2] < 0) ScaleN(basis[2], -1, 3);

Cross3(basis[2], basis[0], cross);
CeedScalar basis_1_orientation = Dot3(cross, basis[1]);
if (basis_1_orientation < 0) ScaleN(basis[1], -1, 3);
}

// @brief Denormalize outputs using min-max (de-)normalization
CEED_QFUNCTION_HELPER void DenormalizeDDOutputs(CeedScalar output[6], const CeedScalar new_bounds[6][2], const CeedScalar old_bounds[6][2]) {
CeedScalar bounds_ratio;
for (int i = 0; i < 6; i++) {
bounds_ratio = (new_bounds[i][1] - new_bounds[i][0]) / (old_bounds[i][1] - old_bounds[i][0]);
output[i] = bounds_ratio * (output[i] - old_bounds[i][1]) + new_bounds[i][1];
}
}

/**
* @brief Compute model inputs for anisotropic data-driven model
*
* @param[in] grad_velo_aniso Gradient of velocity in physical (anisotropic) coordinates
* @param[in] km_A_ij Anisotropy tensor, in Kelvin-Mandel notation
* @param[in] delta Length used to create anisotropy tensor
* @param[in] viscosity Kinematic viscosity
* @param[out] eigenvectors Eigenvectors of the (anisotropic) velocity gradient
* @param[out] inputs Data-driven model inputs
* @param[out] grad_velo_magnitude Frobenius norm of the velocity gradient
*/
CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropicInputs(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) {
CeedScalar strain_sframe[3] = {0.}, vorticity_sframe[3] = {0.};
CeedScalar A_ij[3][3] = {{0.}}, grad_velo_iso[3][3] = {{0.}};

// -- Transform physical, anisotropic velocity gradient to isotropic
KMUnpack(km_A_ij, A_ij);
MatMat3(grad_velo_aniso, A_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, grad_velo_iso);

{ // -- Get Eigenframe
CeedScalar kmstrain_iso[6], strain_iso[3][3];
CeedInt work_vector[3] = {0};
KMStrainRate(grad_velo_iso, kmstrain_iso);
KMUnpack(kmstrain_iso, strain_iso);
Diagonalize3(strain_iso, strain_sframe, eigenvectors, work_vector, SORT_DECREASING_EVALS, true, 5);
}

{ // -- Get vorticity in S-frame
CeedScalar rotation_iso[3][3];
RotationRate(grad_velo_iso, rotation_iso);
CeedScalar vorticity_iso[3] = {-2 * rotation_iso[1][2], 2 * rotation_iso[0][2], -2 * rotation_iso[0][1]};
OrientBasisWithVector(eigenvectors, vorticity_iso);
MatVec3(eigenvectors, vorticity_iso, CEED_NOTRANSPOSE, vorticity_sframe);
}

// -- Calculate DD model inputs
*grad_velo_magnitude = VelocityGradientMagnitude(strain_sframe, vorticity_sframe);
inputs[0] = strain_sframe[0];
inputs[1] = strain_sframe[1];
inputs[2] = strain_sframe[2];
inputs[3] = vorticity_sframe[0];
inputs[4] = vorticity_sframe[1];
inputs[5] = viscosity / Square(delta);
ScaleN(inputs, 1 / (*grad_velo_magnitude + CEED_EPSILON), 6);
}

/**
* @brief Compute the physical SGS stresses from the neural-network output
*
* @param[in,out] outputs Outputs from the neural-network
* @param[in] delta Length used to create anisotropy tensor
* @param[in] eigenvectors Eigenvectors of the (anisotropic) velocity gradient
* @param[in] new_bounds Bounds used for min-max de-normalization
* @param[in] grad_velo_magnitude Magnitude of the velocity gradient
* @param[out] kmsgs_stress Physical SGS stresses in Kelvin-Mandel notation
*/
CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropicOutputs(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]) {
CeedScalar old_bounds[6][2] = {{0}};
for (int j = 0; j < 6; j++) old_bounds[j][1] = 1;
DenormalizeDDOutputs(outputs, new_bounds, old_bounds);

// Re-dimensionalize sgs_stress
ScaleN(outputs, Square(delta) * Square(grad_velo_magnitude), 6);

CeedScalar sgs_stress[3][3] = {{0.}};
{ // Rotate SGS Stress back to physical frame, SGS_physical = E^T SGS_sframe E
CeedScalar Evec_sgs[3][3] = {{0.}};
const CeedScalar sgs_sframe[3][3] = {
{outputs[0], outputs[3], outputs[4]},
{outputs[3], outputs[1], outputs[5]},
{outputs[4], outputs[5], outputs[2]},
};
MatMat3(eigenvectors, sgs_sframe, CEED_TRANSPOSE, CEED_NOTRANSPOSE, Evec_sgs);
MatMat3(Evec_sgs, eigenvectors, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, sgs_stress);
}

KMPack(sgs_stress, kmsgs_stress);
}

#endif // sgs_dd_utils_h