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