Entry point.
91 """!
92 @brief Entry point.
93 @param[in] argv Command-line style argument list supplied to the function.
94 @return Process exit status.
95 """
96 parser = argparse.ArgumentParser(description=__doc__,
97 formatter_class=argparse.RawDescriptionHelpFormatter)
98 parser.add_argument("--checkpoint", required=True,
99 help="Committed checkpoint bundle directory (output/checkpoints/step_NNN).")
100 parser.add_argument("--grid", required=True, help="Canonical PICGRID path for the run.")
101 parser.add_argument("--window", default="stationary", help="Statistics window name.")
102 parser.add_argument("--block", type=int, default=0, help="Block index.")
103 parser.add_argument("--stream-axis", default="Zeta", choices=("Xi", "Eta", "Zeta"),
104 help="Streamwise (driven, homogeneous) axis token.")
105 parser.add_argument("--viscosity", type=float, required=True,
106 help="Kinematic viscosity in solver units.")
107 source = parser.add_mutually_exclusive_group()
108 source.add_argument("--body-force", type=float,
109 help="Converged driving body force f; u_tau = sqrt(f*A/P).")
110 source.add_argument("--u-tau", type=float, help="Friction velocity supplied directly.")
111 parser.add_argument("--no-fold", action="store_true",
112 help="Keep the raw cross-section instead of folding it onto one octant.")
113 parser.add_argument("--output-prefix", required=True, help="Prefix for the four output CSVs.")
114 args = parser.parse_args(argv)
115
116 channel = load_channel_tool()
117 pairs = channel.M2_PAIRS
118 stream_kji = channel.AXIS_TO_KJI[args.stream_axis]
119 stream_c = channel.AXIS_TO_COMPONENT[args.stream_axis]
120 stats = channel.homogeneous_statistics(args.checkpoint, args.window, args.grid, args.block,
121 (stream_kji,))
122 numpy, nodes, window = stats["numpy"], stats["nodes"], stats["window"]
123 mean, covariance = stats["mean"], stats["covariance"]
124
125
126
127 cross_kji = [axis for axis in (0, 1, 2) if axis != stream_kji]
128 kji_to_component = {0: 2, 1: 1, 2: 0}
129 comp = [kji_to_component[axis] for axis in cross_kji]
130 names = "xyz"
131
132 def axis_nodes(kji, component):
133 line = {0: nodes[:, 0, 0, :], 1: nodes[0, :, 0, :], 2: nodes[0, 0, :, :]}[kji]
134 return line[:, component]
135
136 edges = [axis_nodes(kji, c) for kji, c in zip(cross_kji, comp)]
137
138
139 span = [abs(e[-1] - e[0]) for e in edges]
140 for kji, c, e, length in zip(cross_kji, comp, edges, span):
141 full = nodes[..., c]
142 reference = {0: e[:, None, None], 1: e[None, :, None], 2: e[None, None, :]}[kji]
143 if numpy.max(numpy.abs(full - reference)) > 1.0e-9 * length:
144 raise SystemExit("the cross-section is not rectilinear; this reduction assumes it is.")
145 centres = [0.5 * (e[:-1] + e[1:]) for e in edges]
146 widths = [numpy.abs(numpy.diff(e)) for e in edges]
147 n0, n1 = mean.shape[0], mean.shape[1]
148 if (centres[0].size, centres[1].size) != (n0, n1):
149 raise SystemExit(f"cross-section cell counts {(centres[0].size, centres[1].size)} "
150 f"do not match the statistics {(n0, n1)}.")
151
152 def symmetric(c, e, length):
153 return numpy.max(numpy.abs((c - e[0]) - (e[-1] - c[::-1]))) <= 1.0e-9 * length
154
155 folded = "raw (--no-fold)"
156 if not args.no_fold:
157 if not all(symmetric(c, e, length) for c, e, length in zip(centres, edges, span)):
158 raise SystemExit("the cross-section grid is not mirror-symmetric; rerun with --no-fold.")
159 identity = (0, 1, 2)
160 images = []
161 for flip0 in (False, True):
162 for flip1 in (False, True):
163 m, v = mean, covariance
164 sign = [1.0, 1.0, 1.0]
165 if flip0:
166 m, v = m[::-1, :], v[::-1, :]
167 sign[comp[0]] = -1.0
168 if flip1:
169 m, v = m[:, ::-1], v[:, ::-1]
170 sign[comp[1]] = -1.0
171 images.append(transform(numpy, m, v, pairs, identity, sign))
172 diagonal = n0 == n1 and numpy.max(
173 numpy.abs((edges[0] - edges[0][0]) - (edges[1] - edges[1][0]))) <= 1.0e-9 * span[0]
174 if diagonal:
175 swap =
list(identity)
176 swap[comp[0]], swap[comp[1]] = comp[1], comp[0]
177 images += [transform(numpy, numpy.swapaxes(m, 0, 1), numpy.swapaxes(v, 0, 1),
178 pairs, swap, (1.0, 1.0, 1.0)) for m, v in images]
179 mean = sum(m for m, _ in images) / len(images)
180 covariance = sum(v for _, v in images) / len(images)
181 folded = f"folded over {len(images)} symmetry images"
182
183 nu = args.viscosity
184 U = mean[..., stream_c]
185 area = numpy.outer(widths[0], widths[1])
186 bulk = float(numpy.sum(U * area) / numpy.sum(area))
187
188
189 first = [abs(centres[0][0] - edges[0][0]), abs(edges[0][-1] - centres[0][-1]),
190 abs(centres[1][0] - edges[1][0]), abs(edges[1][-1] - centres[1][-1])]
191 walls = [
192 (f"{names[comp[0]]}_min", centres[1], widths[1], nu * U[0, :] / first[0]),
193 (f"{names[comp[0]]}_max", centres[1], widths[1], nu * U[-1, :] / first[1]),
194 (f"{names[comp[1]]}_min", centres[0], widths[0], nu * U[:, 0] / first[2]),
195 (f"{names[comp[1]]}_max", centres[0], widths[0], nu * U[:, -1] / first[3]),
196 ]
197 perimeter = sum(float(numpy.sum(w)) for _, _, w, _ in walls)
198 tau_mean = sum(float(numpy.sum(t * w)) for _, _, w, t in walls) / perimeter
199
200 if args.body_force is not None:
201 cross_area = span[0] * span[1]
202 u_tau = math.sqrt(args.body_force * cross_area / perimeter)
203 source_note = f"sqrt(f*A/P) with f={args.body_force:.8e}"
204 elif args.u_tau is not None:
205 u_tau = args.u_tau
206 source_note = "supplied with --u-tau"
207 else:
208 u_tau = math.sqrt(tau_mean)
209 source_note = "sqrt of the perimeter-mean nu*U/d from the first cells (needs a resolved sublayer)"
210 if not u_tau > 0.0:
211 raise SystemExit("computed a non-positive friction velocity; check the inputs.")
212
213 secondary = numpy.hypot(mean[..., comp[0]], mean[..., comp[1]])
214 peak = numpy.unravel_index(numpy.argmax(secondary), secondary.shape)
215 half_width = 0.5 * min(span)
216 header = (f"# u_tau = {u_tau:.10e} ({source_note})\n"
217 f"# Re_tau = u_tau*a/nu = {u_tau * half_width / nu:.6f} with half-width a = {half_width:.10e}\n"
218 f"# U_b = {bulk:.10e}, max secondary speed / U_b = {float(secondary[peak]) / bulk:.6e} "
219 f"at ({centres[0][peak[0]]:.6e}, {centres[1][peak[1]]:.6e})\n"
220 f"# nu = {nu:.10e}, window = {args.window} ({window.get('state')}, "
221 f"{window.get('sample_count')} samples), {folded}\n")
222 prefix = args.output_prefix
223 pair_names = [names[a] + names[b] for a, b in pairs]
224
225 with open(prefix + "_cross_section.csv", "w", encoding="utf-8") as handle:
226 handle.write(header)
227 handle.write(f"{names[comp[0]]},{names[comp[1]]},u_x,u_y,u_z,"
228 + ",".join("cov_" + p for p in pair_names) + "\n")
229 for a in range(n0):
230 for b in range(n1):
231 handle.write(f"{centres[0][a]:.10e},{centres[1][b]:.10e},"
232 + ",".join(f"{value:.10e}" for value in mean[a, b])
233 + "," + ",".join(f"{value:.10e}" for value in covariance[a, b]) + "\n")
234
235 def cov(a_index, b_index, c0, c1):
236 return covariance[a_index, b_index, pairs.index(tuple(sorted((c0, c1))))]
237
238
239 middle = [n1 // 2] if n1 % 2 else [n1 // 2 - 1, n1 // 2]
240 orientation = 1.0 if edges[0][-1] > edges[0][0] else -1.0
241 with open(prefix + "_wall_bisector.csv", "w", encoding="utf-8") as handle:
242 handle.write(header)
243 handle.write("d,d_plus,U_plus,u_rms_plus,v_rms_plus,w_rms_plus,minus_uv_plus,V_over_Ub\n")
244 for a in range((n0 + 1) // 2):
245 pick = lambda f: sum(f(a, b) for b in middle) / len(middle)
246 d = abs(centres[0][a] - edges[0][0])
247 values = (
248 pick(lambda i, j: U[i, j]) / u_tau,
249 math.sqrt(max(pick(lambda i, j: cov(i, j, stream_c, stream_c)), 0.0)) / u_tau,
250 math.sqrt(max(pick(lambda i, j: cov(i, j, comp[0], comp[0])), 0.0)) / u_tau,
251 math.sqrt(max(pick(lambda i, j: cov(i, j, comp[1], comp[1])), 0.0)) / u_tau,
252 -orientation * pick(lambda i, j: cov(i, j, stream_c, comp[0])) / u_tau ** 2,
253 orientation * pick(lambda i, j: mean[i, j, comp[0]]) / bulk,
254 )
255 handle.write(f"{d:.10e},{d * u_tau / nu:.10e}," + ",".join(f"{v:.10e}" for v in values) + "\n")
256
257 if n0 == n1:
258
259
260 sign0 = 1.0 if edges[0][-1] > edges[0][0] else -1.0
261 sign1 = 1.0 if edges[1][-1] > edges[1][0] else -1.0
262 with open(prefix + "_corner_bisector.csv", "w", encoding="utf-8") as handle:
263 handle.write(header)
264 handle.write("d,d_plus,U_plus,along_diagonal_over_Ub,k_plus\n")
265 for a in range((n0 + 1) // 2):
266 d = math.hypot(centres[0][a] - edges[0][0], centres[1][a] - edges[1][0])
267 along = (sign0 * mean[a, a, comp[0]] + sign1 * mean[a, a, comp[1]]) / math.sqrt(2.0)
268 k = 0.5 * sum(cov(a, a, c, c) for c in range(3))
269 handle.write(f"{d:.10e},{d * u_tau / nu:.10e},{U[a, a] / u_tau:.10e},"
270 f"{along / bulk:.10e},{k / u_tau ** 2:.10e}\n")
271
272 with open(prefix + "_wall_shear.csv", "w", encoding="utf-8") as handle:
273 handle.write(header)
274 handle.write("wall,s,tau_w,tau_w_over_mean\n")
275 for name, coordinate, _width, tau in walls:
276 for s_value, t_value in zip(coordinate, tau):
277 handle.write(f"{name},{s_value:.10e},{t_value:.10e},{t_value / tau_mean:.10e}\n")
278
279 print(header, end="")
280 print(f"wrote {prefix}_cross_section.csv, {prefix}_wall_bisector.csv, "
281 + (f"{prefix}_corner_bisector.csv, " if n0 == n1 else "")
282 + f"{prefix}_wall_shear.csv")
283 return 0
284
285
int main(int argc, char **argv)
Entry point for the postprocessor executable.
Head of a generic C-style linked list.
1.9.8