Skip to content

Commit

Permalink
fluids: Extract out the DD output processing
Browse files Browse the repository at this point in the history
  • Loading branch information
jrwrigh committed Aug 23, 2023
1 parent 1abde36 commit b32ff2a
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 32 deletions.
34 changes: 2 additions & 32 deletions examples/fluids/qfunctions/sgs_dd_model.h
Original file line number Diff line number Diff line change
Expand Up @@ -38,15 +38,6 @@ struct SGS_DD_ModelContext_ {
CeedScalar data[1];
};

// @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];
}
}

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 @@ -72,32 +63,11 @@ 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 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];

ComputeSGS_DDAnisotropicInputs(grad_velo_aniso, km_A_ij, delta, viscosity, eigenvectors, inputs, &grad_velo_magnitude);

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;
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_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
44 changes: 44 additions & 0 deletions examples/fluids/qfunctions/sgs_dd_utils.h
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,15 @@ CEED_QFUNCTION_HELPER void OrientBasisWithVector(CeedScalar basis[3][3], const C
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
*
Expand Down Expand Up @@ -98,4 +107,39 @@ CEED_QFUNCTION_HELPER void ComputeSGS_DDAnisotropicInputs(const CeedScalar grad_
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

0 comments on commit b32ff2a

Please sign in to comment.