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 DifferentialFilterStateComponent { 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_STATE_NUM, 24 }; 25 26 enum DifferentialFilterVelocitySquared { 27 DIFF_FILTER_VELOCITY_SQUARED_XX, 28 DIFF_FILTER_VELOCITY_SQUARED_YY, 29 DIFF_FILTER_VELOCITY_SQUARED_ZZ, 30 DIFF_FILTER_VELOCITY_SQUARED_YZ, 31 DIFF_FILTER_VELOCITY_SQUARED_XZ, 32 DIFF_FILTER_VELOCITY_SQUARED_XY, 33 DIFF_FILTER_VELOCITY_SQUARED_NUM, 34 }; 35 36 enum DifferentialFilterDampingFunction { DIFF_FILTER_DAMP_NONE, DIFF_FILTER_DAMP_VAN_DRIEST, DIFF_FILTER_DAMP_MMS }; 37 static const char *const DifferentialFilterDampingFunctions[] = { 38 "none", "van_driest", "mms", "DifferentialFilterDampingFunction", "DIFF_FILTER_DAMP_", NULL}; 39 40 typedef struct DifferentialFilterContext_ *DifferentialFilterContext; 41 struct DifferentialFilterContext_ { 42 bool grid_based_width; 43 CeedScalar width_scaling[3]; 44 CeedScalar kernel_scaling; 45 CeedScalar friction_length; 46 enum DifferentialFilterDampingFunction damping_function; 47 CeedScalar damping_constant; 48 struct NewtonianIdealGasContext_ gas; 49 }; 50 51 CEED_QFUNCTION_HELPER int DifferentialFilter_RHS(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, StateVariable state_var) { 52 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 53 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 54 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 55 CeedScalar(*v0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 56 CeedScalar(*v1)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 57 58 DifferentialFilterContext context = (DifferentialFilterContext)ctx; 59 NewtonianIdealGasContext gas = &context->gas; 60 61 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 62 const CeedScalar qi[5] = {q[0][i], q[1][i], q[2][i], q[3][i], q[4][i]}; 63 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 64 const CeedScalar wdetJ = q_data[0][i]; 65 const State s = StateFromQ(gas, qi, x_i, state_var); 66 67 v0[DIFF_FILTER_PRESSURE][i] = wdetJ * s.Y.pressure; 68 v0[DIFF_FILTER_VELOCITY_X][i] = wdetJ * s.Y.velocity[0]; 69 v0[DIFF_FILTER_VELOCITY_Y][i] = wdetJ * s.Y.velocity[1]; 70 v0[DIFF_FILTER_VELOCITY_Z][i] = wdetJ * s.Y.velocity[2]; 71 v0[DIFF_FILTER_TEMPERATURE][i] = wdetJ * s.Y.temperature; 72 v1[DIFF_FILTER_VELOCITY_SQUARED_XX][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[0]; 73 v1[DIFF_FILTER_VELOCITY_SQUARED_YY][i] = wdetJ * s.Y.velocity[1] * s.Y.velocity[1]; 74 v1[DIFF_FILTER_VELOCITY_SQUARED_ZZ][i] = wdetJ * s.Y.velocity[2] * s.Y.velocity[2]; 75 v1[DIFF_FILTER_VELOCITY_SQUARED_YZ][i] = wdetJ * s.Y.velocity[1] * s.Y.velocity[2]; 76 v1[DIFF_FILTER_VELOCITY_SQUARED_XZ][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[2]; 77 v1[DIFF_FILTER_VELOCITY_SQUARED_XY][i] = wdetJ * s.Y.velocity[0] * s.Y.velocity[1]; 78 } 79 return 0; 80 } 81 82 CEED_QFUNCTION(DifferentialFilter_RHS_Conserv)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 83 return DifferentialFilter_RHS(ctx, Q, in, out, STATEVAR_CONSERVATIVE); 84 } 85 86 CEED_QFUNCTION(DifferentialFilter_RHS_Prim)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 87 return DifferentialFilter_RHS(ctx, Q, in, out, STATEVAR_PRIMITIVE); 88 } 89 90 CEED_QFUNCTION_HELPER CeedScalar VanDriestWallDamping(const CeedScalar wall_dist_plus, const CeedScalar A_plus) { 91 return -expm1(-wall_dist_plus / A_plus); 92 } 93 94 CEED_QFUNCTION_HELPER int DifferentialFilter_LHS_N(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out, const CeedInt N) { 95 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 96 const CeedScalar(*Grad_q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 97 const CeedScalar(*A_ij_delta)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2]; 98 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[3]; 99 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[4]; 100 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 101 CeedScalar(*Grad_v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[1]; 102 103 DifferentialFilterContext context = (DifferentialFilterContext)ctx; 104 105 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 106 CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { 107 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 108 const CeedScalar wdetJ = q_data[0][i]; 109 const CeedScalar dXdx[3][3] = { 110 {q_data[1][i], q_data[2][i], q_data[3][i]}, 111 {q_data[4][i], q_data[5][i], q_data[6][i]}, 112 {q_data[7][i], q_data[8][i], q_data[9][i]} 113 }; 114 115 CeedScalar Delta_ij[3][3] = {{0.}}; 116 if (context->grid_based_width) { 117 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]}; 118 const CeedScalar delta = A_ij_delta[6][i]; 119 ScaleN(km_A_ij, delta, 6); // Dimensionalize the normalized anisotropy tensor 120 KMUnpack(km_A_ij, Delta_ij); 121 } else { 122 Delta_ij[0][0] = Delta_ij[1][1] = Delta_ij[2][2] = 1; 123 } 124 125 CeedScalar scaling_matrix[3][3] = {{0}}; 126 if (context->damping_function == DIFF_FILTER_DAMP_VAN_DRIEST) { 127 const CeedScalar damping_coeff = VanDriestWallDamping(x_i[1] / context->friction_length, context->damping_constant); 128 scaling_matrix[0][0] = Max(1, damping_coeff * context->width_scaling[0]); 129 scaling_matrix[1][1] = damping_coeff * context->width_scaling[1]; 130 scaling_matrix[2][2] = Max(1, damping_coeff * context->width_scaling[2]); 131 } else if (context->damping_function == DIFF_FILTER_DAMP_NONE) { 132 scaling_matrix[0][0] = context->width_scaling[0]; 133 scaling_matrix[1][1] = context->width_scaling[1]; 134 scaling_matrix[2][2] = context->width_scaling[2]; 135 } else if (context->damping_function == DIFF_FILTER_DAMP_MMS) { 136 const CeedScalar damping_coeff = tanh(60 * x_i[1]); 137 scaling_matrix[0][0] = 1; 138 scaling_matrix[1][1] = damping_coeff; 139 scaling_matrix[2][2] = 1; 140 } 141 142 CeedScalar scaled_Delta_ij[3][3] = {{0.}}; 143 MatMat3(scaling_matrix, Delta_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, scaled_Delta_ij); 144 CopyMat3(scaled_Delta_ij, Delta_ij); 145 146 CeedScalar alpha_ij[3][3] = {{0.}}; 147 MatMat3(Delta_ij, Delta_ij, CEED_NOTRANSPOSE, CEED_NOTRANSPOSE, alpha_ij); 148 ScaleN((CeedScalar *)alpha_ij, context->kernel_scaling, 9); 149 150 v[j][i] = wdetJ * q[j][i]; 151 CeedScalar dq[3], dq_dXdx[3] = {0.}, dq_dXdx_a[3] = {0.}; 152 for (int k = 0; k < 3; k++) { 153 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]; 154 } 155 MatVec3(dXdx, dq, CEED_NOTRANSPOSE, dq_dXdx); 156 MatVec3(alpha_ij, dq_dXdx, CEED_NOTRANSPOSE, dq_dXdx_a); 157 for (int k = 0; k < 3; k++) { 158 Grad_v[k * N + j][i] = wdetJ * dq_dXdx_a[k]; 159 } 160 } 161 } 162 return 0; 163 } 164 165 CEED_QFUNCTION(DifferentialFilter_LHS_1)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 166 return DifferentialFilter_LHS_N(ctx, Q, in, out, 1); 167 } 168 169 CEED_QFUNCTION(DifferentialFilter_LHS_5)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 170 return DifferentialFilter_LHS_N(ctx, Q, in, out, 5); 171 } 172 173 CEED_QFUNCTION(DifferentialFilter_LHS_6)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 174 return DifferentialFilter_LHS_N(ctx, Q, in, out, 6); 175 } 176 177 CEED_QFUNCTION(DifferentialFilter_LHS_11)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 178 return DifferentialFilter_LHS_N(ctx, Q, in, out, 11); 179 } 180 181 CEED_QFUNCTION_HELPER CeedScalar MMS_Solution(const CeedScalar x_i[3], const CeedScalar omega) { 182 return sin(6 * omega * x_i[0]) + sin(6 * omega * x_i[1]); 183 } 184 185 CEED_QFUNCTION(DifferentialFilter_MMS_RHS)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 186 const CeedScalar(*q)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 187 const CeedScalar(*q_data)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[1]; 188 CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 189 190 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 191 const CeedScalar wdetJ = q_data[0][i]; 192 v[0][i] = wdetJ * q[0][i]; 193 } 194 return 0; 195 } 196 197 // @brief Generate initial condition such that the solution of the differential filtering is given by MMS_Solution() above 198 // 199 // This requires a *very* specific grid, as the anisotropic filtering is grid dependent. 200 // It's for a 75x75x1 grid on a [0,0.5]x3 domain. 201 // The grid is evenly distributed in x, but distributed based on the analytical mesh distribution \Delta_y = .001 + .01*tanh(6*y). 202 // The MMS test can optionally include a wall damping function (must also be enabled for the differential filtering itself). 203 // It can be run via: 204 // ./navierstokes -options_file tests-output/blasius_test.yaml -diff_filter_enable -diff_filter_view cgns:filtered_solution.cgns -ts_max_steps 0 205 // -diff_filter_mms -diff_filter_kernel_scaling 1 -diff_filter_wall_damping_function mms -dm_plex_box_upper 0.5,0.5,0.5 -dm_plex_box_faces 75,75,1 206 // -platemesh_y_node_locs_path tests-output/diff_filter_mms_y_spacing.dat -platemesh_top_angle 0 -diff_filter_grid_based_width 207 CEED_QFUNCTION(DifferentialFilter_MMS_IC)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { 208 const CeedScalar(*x)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; 209 CeedScalar(*q0)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; 210 211 CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { 212 const CeedScalar x_i[3] = {x[0][i], x[1][i], x[2][i]}; 213 214 const CeedScalar aniso_scale_factor = 1; // Must match the one passed in by -diff_filter_aniso_scale 215 const CeedScalar omega = 2 * M_PI; 216 const CeedScalar omega6 = 6 * omega; 217 const CeedScalar phi_bar = MMS_Solution(x_i, omega); 218 const CeedScalar dx = 0.5 / 75; 219 const CeedScalar dy_analytic = .001 + .01 * tanh(6 * x_i[1]); 220 const CeedScalar dy = dy_analytic; 221 const CeedScalar d_dy_dy = 0.06 * Square(1 / cosh(6 * x_i[1])); // Change of \Delta_y w.r.t. y 222 CeedScalar alpha[2] = {Square(dx) * aniso_scale_factor, Square(dy) * aniso_scale_factor}; 223 bool damping = true; 224 CeedScalar dalpha1dy; 225 if (damping) { 226 CeedScalar damping_coeff = tanh(60 * x_i[1]); 227 CeedScalar d_damping_coeff = 60 / Square(cosh(60 * x_i[1])); 228 dalpha1dy = aniso_scale_factor * 2 * (damping_coeff * dy) * (dy * d_damping_coeff + damping_coeff * d_dy_dy); 229 alpha[1] *= Square(damping_coeff); 230 } else { 231 dalpha1dy = aniso_scale_factor * 2 * dy * d_dy_dy; 232 } 233 234 CeedScalar phi = phi_bar + alpha[0] * Square(omega6) * sin(6 * omega * x_i[0]) + alpha[1] * Square(omega6) * sin(omega6 * x_i[1]); 235 phi -= dalpha1dy * omega6 * cos(omega6 * x_i[1]); 236 237 for (CeedInt j = 0; j < 5; j++) q0[j][i] = phi; 238 } 239 return 0; 240 } 241