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 bool grid_based_width; 35 CeedScalar width_scaling[3]; 36 CeedScalar kernel_scaling; 37 struct NewtonianIdealGasContext_ gas; 38 }; 39 40 CEED_QFUNCTION_HELPER int DifferentialFilter_RHS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 41 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 42 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 43 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 44 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 45 46 DifferentialFilterContext context = (DifferentialFilterContext)ctx; 47 NewtonianIdealGasContext gas = &context->gas; 48 49 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 50 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 51 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 52 const CeedScalar wdetJ = q_data[0][i]; 53 const State s = StateFromQ(gas, qi, x_i, state_var); 54 55 v[DIFF_FILTER_PRESSURE][i] = wdetJ * s.Y.pressure; 56 v[DIFF_FILTER_VELOCITY_X][i] = wdetJ * s.Y.velocity[0]; 57 v[DIFF_FILTER_VELOCITY_Y][i] = wdetJ * s.Y.velocity[1]; 58 v[DIFF_FILTER_VELOCITY_Z][i] = wdetJ * s.Y.velocity[2]; 59 v[DIFF_FILTER_TEMPERATURE][i] = wdetJ * s.Y.temperature; 60 v[DIFF_FILTER_VELOCITY_SQUARED_XX][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[0]; 61 v[DIFF_FILTER_VELOCITY_SQUARED_YY][i] = wdetJ * s.Y.velocity[1] * s.Y.velocity[1]; 62 v[DIFF_FILTER_VELOCITY_SQUARED_ZZ][i] = wdetJ * s.Y.velocity[2] * s.Y.velocity[2]; 63 v[DIFF_FILTER_VELOCITY_SQUARED_YZ][i] = wdetJ * s.Y.velocity[1] * s.Y.velocity[2]; 64 v[DIFF_FILTER_VELOCITY_SQUARED_XZ][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[2]; 65 v[DIFF_FILTER_VELOCITY_SQUARED_XY][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[1]; 66 } 67 return 0; 68 } 69 70 CEED_QFUNCTION(DifferentialFilter_RHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 71 return DifferentialFilter_RHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 72 } 73 74 CEED_QFUNCTION(DifferentialFilter_RHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 75 return DifferentialFilter_RHS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 76 } 77 78 CEED_QFUNCTION_HELPER int DifferentialFilter_LHS_N(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, const CeedInt N) { 79 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 80 const CeedScalar(*Grad_q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 81 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 82 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 83 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 84 CeedScalar(*Grad_v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 85 86 DifferentialFilterContext context = (DifferentialFilterContext)ctx; 87 88 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 89 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 90 const CeedScalar wdetJ = q_data[0][i]; 91 const CeedScalar dXdx[3][3] = { 92 {q_data[1][i], q_data[2][i], q_data[3][i]}, 93 {q_data[4][i], q_data[5][i], q_data[6][i]}, 94 {q_data[7][i], q_data[8][i], q_data[9][i]} 95 }; 96 97 CeedScalar Delta_ij[3][3] = {{0.}}; 98 if (context->grid_based_width) { 99 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]}; 100 const CeedScalar delta = A_ij_delta[6][i]; 101 ScaleN(km_A_ij, delta, 6); // Dimensionalize the normalized anisotropy tensor 102 KMUnpack(km_A_ij, Delta_ij); 103 } else { 104 Delta_ij[0][0] = Delta_ij[1][1] = Delta_ij[2][2] = 1; 105 } 106 107 CeedScalar scaling_matrix[3][3] = {{0.}}; 108 scaling_matrix[0][0] = context->width_scaling[0]; 109 scaling_matrix[1][1] = context->width_scaling[1]; 110 scaling_matrix[2][2] = context->width_scaling[2]; 111 112 CeedScalar scaled_Delta_ij[3][3] = {{0.}}; 113 MatMat3(scaling_matrix, Delta_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, scaled_Delta_ij); 114 CopyMat3(scaled_Delta_ij, Delta_ij); 115 116 CeedScalar alpha_ij[3][3] = {{0.}}; 117 MatMat3(Delta_ij, Delta_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, alpha_ij); 118 ScaleN((CeedScalar *)alpha_ij, context->kernel_scaling, 9); 119 120 v[j][i] = wdetJ * q[j][i]; 121 CeedScalar dq[3], dq_dXdx[3] = {0.}, dq_dXdx_a[3] = {0.}; 122 for (int k = 0; k < 3; k++) { 123 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]; 124 } 125 MatVec3(dXdx, dq, CEED_NOTRANSPOSE, dq_dXdx); 126 MatVec3(alpha_ij, dq_dXdx, CEED_NOTRANSPOSE, dq_dXdx_a); 127 for (int k = 0; k < 3; k++) { 128 Grad_v[k * N + j][i] = wdetJ * dq_dXdx_a[k]; 129 } 130 } 131 } 132 return 0; 133 } 134 135 CEED_QFUNCTION(DifferentialFilter_LHS_1)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 136 return DifferentialFilter_LHS_N(ctx, Q, in, out, 1); 137 } 138 139 CEED_QFUNCTION(DifferentialFilter_LHS_5)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 140 return DifferentialFilter_LHS_N(ctx, Q, in, out, 5); 141 } 142 143 CEED_QFUNCTION(DifferentialFilter_LHS_11)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 144 return DifferentialFilter_LHS_N(ctx, Q, in, out, 11); 145 } 146