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 /// Implementation of differential filtering 10 11 #include <ceed.h> 12 13 #include "newtonian_state.h" 14 #include "newtonian_types.h" 15 #include "utils.h" 16 17 enum DifferentialFilterComponent { 18 DIFF_FILTER_PRESSURE, 19 DIFF_FILTER_VELOCITY_X, 20 DIFF_FILTER_VELOCITY_Y, 21 DIFF_FILTER_VELOCITY_Z, 22 DIFF_FILTER_TEMPERATURE, 23 DIFF_FILTER_VELOCITY_SQUARED_XX, 24 DIFF_FILTER_VELOCITY_SQUARED_YY, 25 DIFF_FILTER_VELOCITY_SQUARED_ZZ, 26 DIFF_FILTER_VELOCITY_SQUARED_YZ, 27 DIFF_FILTER_VELOCITY_SQUARED_XZ, 28 DIFF_FILTER_VELOCITY_SQUARED_XY, 29 DIFF_FILTER_NUM_COMPONENTS, 30 }; 31 32 typedef struct DifferentialFilterContext_ *DifferentialFilterContext; 33 struct DifferentialFilterContext_ { 34 CeedScalar delta; 35 CeedScalar filter_scaling; 36 struct NewtonianIdealGasContext_ gas; 37 }; 38 39 CEED_QFUNCTION_HELPER int DifferentialFilter_RHS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 40 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 41 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 42 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 43 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 44 45 DifferentialFilterContext context = (DifferentialFilterContext)ctx; 46 NewtonianIdealGasContext gas = &context->gas; 47 48 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 49 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 50 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 51 const CeedScalar wdetJ = q_data[0][i]; 52 const State s = StateFromQ(gas, qi, x_i, state_var); 53 54 v[DIFF_FILTER_PRESSURE][i] = wdetJ * s.Y.pressure; 55 v[DIFF_FILTER_VELOCITY_X][i] = wdetJ * s.Y.velocity[0]; 56 v[DIFF_FILTER_VELOCITY_Y][i] = wdetJ * s.Y.velocity[1]; 57 v[DIFF_FILTER_VELOCITY_Z][i] = wdetJ * s.Y.velocity[2]; 58 v[DIFF_FILTER_TEMPERATURE][i] = wdetJ * s.Y.temperature; 59 v[DIFF_FILTER_VELOCITY_SQUARED_XX][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[0]; 60 v[DIFF_FILTER_VELOCITY_SQUARED_YY][i] = wdetJ * s.Y.velocity[1] * s.Y.velocity[1]; 61 v[DIFF_FILTER_VELOCITY_SQUARED_ZZ][i] = wdetJ * s.Y.velocity[2] * s.Y.velocity[2]; 62 v[DIFF_FILTER_VELOCITY_SQUARED_YZ][i] = wdetJ * s.Y.velocity[1] * s.Y.velocity[2]; 63 v[DIFF_FILTER_VELOCITY_SQUARED_XZ][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[2]; 64 v[DIFF_FILTER_VELOCITY_SQUARED_XY][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[1]; 65 } 66 return 0; 67 } 68 69 CEED_QFUNCTION(DifferentialFilter_RHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 70 return DifferentialFilter_RHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 71 } 72 73 CEED_QFUNCTION(DifferentialFilter_RHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 74 return DifferentialFilter_RHS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 75 } 76 77 CEED_QFUNCTION_HELPER int DifferentialFilter_LHS_N(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, const CeedInt N) { 78 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 79 const CeedScalar(*Grad_q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 80 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 81 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 82 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 83 CeedScalar(*Grad_v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 84 85 DifferentialFilterContext context = (DifferentialFilterContext)ctx; 86 87 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 88 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 89 const CeedScalar wdetJ = q_data[0][i]; 90 const CeedScalar dXdx[3][3] = { 91 {q_data[1][i], q_data[2][i], q_data[3][i]}, 92 {q_data[4][i], q_data[5][i], q_data[6][i]}, 93 {q_data[7][i], q_data[8][i], q_data[9][i]} 94 }; 95 96 CeedScalar Delta_ij[3][3] = {{0.}}; 97 if (context->delta > 0) { 98 Delta_ij[0][0] = Delta_ij[1][1] = Delta_ij[2][2] = context->delta; 99 } else { 100 CeedScalar km_A_ij[6] = {A_ij_delta[0][i], A_ij_delta[1][i], A_ij_delta[2][i], A_ij_delta[3][i], A_ij_delta[4][i], A_ij_delta[5][i]}; 101 const CeedScalar delta = A_ij_delta[6][i]; 102 ScaleN(km_A_ij, delta, 6); // Dimensionalize the normalized anisotropy tensor 103 KMUnpack(km_A_ij, Delta_ij); 104 } 105 106 CeedScalar alpha_ij[3][3] = {{0.}}; 107 MatMat3(Delta_ij, Delta_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, alpha_ij); 108 ScaleN((CeedScalar *)alpha_ij, context->filter_scaling, 9); 109 110 v[j][i] = wdetJ * q[j][i]; 111 CeedScalar dq[3], dq_dXdx[3] = {0.}, dq_dXdx_a[3] = {0.}; 112 for (int k = 0; k < 3; k++) { 113 dq[k] = Grad_q[0 * N + j][i] * dXdx[0][k] + Grad_q[1 * N + j][i] * dXdx[1][k] + Grad_q[2 * N + j][i] * dXdx[2][k]; 114 } 115 MatVec3(dXdx, dq, CEED_NOTRANSPOSE, dq_dXdx); 116 MatVec3(alpha_ij, dq_dXdx, CEED_NOTRANSPOSE, dq_dXdx_a); 117 for (int k = 0; k < 3; k++) { 118 Grad_v[k * N + j][i] = wdetJ * dq_dXdx_a[k]; 119 } 120 } 121 } 122 return 0; 123 } 124 125 CEED_QFUNCTION(DifferentialFilter_LHS_1)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 126 return DifferentialFilter_LHS_N(ctx, Q, in, out, 1); 127 } 128 129 CEED_QFUNCTION(DifferentialFilter_LHS_5)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 130 return DifferentialFilter_LHS_N(ctx, Q, in, out, 5); 131 } 132 133 CEED_QFUNCTION(DifferentialFilter_LHS_11)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 134 return DifferentialFilter_LHS_N(ctx, Q, in, out, 11); 135 } 136